summaryrefslogtreecommitdiff
path: root/libm/double/simq.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/simq.c')
-rw-r--r--libm/double/simq.c180
1 files changed, 0 insertions, 180 deletions
diff --git a/libm/double/simq.c b/libm/double/simq.c
deleted file mode 100644
index 96d63e521..000000000
--- a/libm/double/simq.c
+++ /dev/null
@@ -1,180 +0,0 @@
-/* simq.c
- *
- * Solution of simultaneous linear equations AX = B
- * by Gaussian elimination with partial pivoting
- *
- *
- *
- * SYNOPSIS:
- *
- * double A[n*n], B[n], X[n];
- * int n, flag;
- * int IPS[];
- * int simq();
- *
- * ercode = simq( A, B, X, n, flag, IPS );
- *
- *
- *
- * DESCRIPTION:
- *
- * B, X, IPS are vectors of length n.
- * A is an n x n matrix (i.e., a vector of length n*n),
- * stored row-wise: that is, A(i,j) = A[ij],
- * where ij = i*n + j, which is the transpose of the normal
- * column-wise storage.
- *
- * The contents of matrix A are destroyed.
- *
- * Set flag=0 to solve.
- * Set flag=-1 to do a new back substitution for different B vector
- * using the same A matrix previously reduced when flag=0.
- *
- * The routine returns nonzero on error; messages are printed.
- *
- *
- * ACCURACY:
- *
- * Depends on the conditioning (range of eigenvalues) of matrix A.
- *
- *
- * REFERENCE:
- *
- * Computer Solution of Linear Algebraic Systems,
- * by George E. Forsythe and Cleve B. Moler; Prentice-Hall, 1967.
- *
- */
-
-/* simq 2 */
-
-#include <stdio.h>
-#define fabs(x) ((x) < 0 ? -(x) : (x))
-
-int simq( A, B, X, n, flag, IPS )
-double A[], B[], X[];
-int n, flag;
-int IPS[];
-{
-int i, j, ij, ip, ipj, ipk, ipn;
-int idxpiv, iback;
-int k, kp, kp1, kpk, kpn;
-int nip, nkp, nm1;
-double em, q, rownrm, big, size, pivot, sum;
-
-nm1 = n-1;
-if( flag < 0 )
- goto solve;
-
-/* Initialize IPS and X */
-
-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 )
- {
- printf("SIMQ ROWNRM=0");
- return(1);
- }
- X[i] = 1.0/rownrm;
- }
-
-/* simq 3 */
-/* Gaussian elimination with partial pivoting */
-
-for( k=0; k<nm1; k++ )
- {
- big= 0.0;
- idxpiv = 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 )
- {
- printf( "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; /* last element of IPS[n] th row */
-if( A[kpn] == 0.0 )
- {
- printf( "SIMQ A[kpn]=0");
- return(3);
- }
-
-/* simq 4 */
-/* back substitution */
-
-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 goes (n-1),...,1 */
- 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);
-}