#include <stdio.h>
#include <math.h>
int
simq(
double
A[],
double
B[],
double
X[],
int
n,
int
flag,
int
IPS[] )
{
int
i, j, ij, ip, ipj, ipk, ipn;
int
idxpiv = 0, iback;
int
k, kp, kp1, kpk, kpn;
int
nip, nkp, nm1 = n-1;
double
em, q, rownrm, big, size, pivot, sum;
if
( flag < 0 )
goto
solve;
ij=0;
for
( i=0; i<n; i++ )
{
IPS[i] = i;
rownrm = 0.0;
for
( j=0; j<n; j++ )
{
q =
fabs
( A[ij] );
if
( rownrm < q )
rownrm = q;
++ij;
}
if
( rownrm == 0.0 )
{
puts
(
"SIMQ ROWNRM=0"
);
return
(1);
}
X[i] = 1.0/rownrm;
}
for
( k=0; k<nm1; k++ )
{
big= 0.0;
for
( i=k; i<n; i++ )
{
ip = IPS[i];
ipk = n*ip + k;
size =
fabs
( A[ipk] ) * X[ip];
if
( size > big )
{
big = size;
idxpiv = i;
}
}
if
( big == 0.0 )
{
puts
(
"SIMQ BIG=0"
);
return
(2);
}
if
( idxpiv != k )
{
j = IPS[k];
IPS[k] = IPS[idxpiv];
IPS[idxpiv] = j;
}
kp = IPS[k];
kpk = n*kp + k;
pivot = A[kpk];
kp1 = k+1;
for
( i=kp1; i<n; i++ )
{
ip = IPS[i];
ipk = n*ip + k;
em = -A[ipk]/pivot;
A[ipk] = -em;
nip = n*ip;
nkp = n*kp;
for
( j=kp1; j<n; j++ )
{
ipj = nip + j;
A[ipj] = A[ipj] + em * A[nkp + j];
}
}
}
kpn = n * IPS[n-1] + n - 1;
if
( A[kpn] == 0.0 )
{
puts
(
"SIMQ A[kpn]=0"
);
return
(3);
}
solve:
ip = IPS[0];
X[0] = B[ip];
for
( i=1; i<n; i++ )
{
ip = IPS[i];
ipj = n * ip;
sum = 0.0;
for
( j=0; j<i; j++ )
{
sum += A[ipj] * X[j];
++ipj;
}
X[i] = B[ip] - sum;
}
ipn = n * IPS[n-1] + n - 1;
X[n-1] = X[n-1]/A[ipn];
for
( iback=1; iback<n; iback++ )
{
i = nm1 - iback;
ip = IPS[i];
nip = n*ip;
sum = 0.0;
for
( j=i+1; j<n; j++ )
sum += A[nip+j] * X[j];
X[i] = (X[i] - sum)/A[nip+i];
}
return
(0);
}