File Coverage

lib/PDL/MatrixOps/simq.c
Criterion Covered Total %
statement 67 74 90.5
branch 30 34 88.2
condition n/a
subroutine n/a
pod n/a
total 97 108 89.8


line stmt bran cond sub pod time code
1             /* simq.c
2             *
3             * Solution of simultaneous linear equations AX = B
4             * by Gaussian elimination with partial pivoting
5             *
6             *
7             *
8             * SYNOPSIS:
9             *
10             * double A[n*n], B[n], X[n];
11             * int n, flag;
12             * int IPS[];
13             * int simq();
14             *
15             * ercode = simq( A, B, X, n, flag, IPS );
16             *
17             *
18             *
19             * DESCRIPTION:
20             *
21             * B, X, IPS are vectors of length n.
22             * A is an n x n matrix (i.e., a vector of length n*n),
23             * stored row-wise: that is, A(i,j) = A[ij],
24             * where ij = i*n + j, which is the transpose of the normal
25             * column-wise storage.
26             *
27             * The contents of matrix A are destroyed.
28             *
29             * Set flag=0 to solve.
30             * Set flag=-1 to do a new back substitution for different B vector
31             * using the same A matrix previously reduced when flag=0.
32             *
33             * The routine returns nonzero on error; messages are printed.
34             *
35             *
36             * ACCURACY:
37             *
38             * Depends on the conditioning (range of eigenvalues) of matrix A.
39             *
40             *
41             * REFERENCE:
42             *
43             * Computer Solution of Linear Algebraic Systems,
44             * by George E. Forsythe and Cleve B. Moler; Prentice-Hall, 1967.
45             *
46             */
47            
48             /* simq 2 */
49              
50             #include
51             #include
52              
53 4           int simq( double A[], double B[], double X[], int n, int flag, int IPS[] )
54             {
55             int i, j, ij, ip, ipj, ipk, ipn;
56 4           int idxpiv = 0, iback;
57             int k, kp, kp1, kpk, kpn;
58 4           int nip, nkp, nm1 = n-1;
59             double em, q, rownrm, big, size, pivot, sum;
60              
61 4 50         if( flag < 0 )
62 0           goto solve;
63              
64             /* Initialize IPS and X */
65              
66 4           ij=0;
67 15 100         for( i=0; i
68             {
69 11           IPS[i] = i;
70 11           rownrm = 0.0;
71 44 100         for( j=0; j
72             {
73 33           q = fabs( A[ij] );
74 33 100         if( rownrm < q )
75 18           rownrm = q;
76 33           ++ij;
77             }
78 11 50         if( rownrm == 0.0 )
79             {
80 0           puts("SIMQ ROWNRM=0");
81 0           return(1);
82             }
83 11           X[i] = 1.0/rownrm;
84             }
85            
86             /* simq 3 */
87             /* Gaussian elimination with partial pivoting */
88              
89 11 100         for( k=0; k
90             {
91 7           big= 0.0;
92 25 100         for( i=k; i
93             {
94 18           ip = IPS[i];
95 18           ipk = n*ip + k;
96 18           size = fabs( A[ipk] ) * X[ip];
97 18 100         if( size > big )
98             {
99 7           big = size;
100 7           idxpiv = i;
101             }
102             }
103              
104 7 50         if( big == 0.0 )
105             {
106 0           puts( "SIMQ BIG=0" );
107 0           return(2);
108             }
109 7 100         if( idxpiv != k )
110             {
111 2           j = IPS[k];
112 2           IPS[k] = IPS[idxpiv];
113 2           IPS[idxpiv] = j;
114             }
115 7           kp = IPS[k];
116 7           kpk = n*kp + k;
117 7           pivot = A[kpk];
118 7           kp1 = k+1;
119 18 100         for( i=kp1; i
120             {
121 11           ip = IPS[i];
122 11           ipk = n*ip + k;
123 11           em = -A[ipk]/pivot;
124 11           A[ipk] = -em;
125 11           nip = n*ip;
126 11           nkp = n*kp;
127 32 100         for( j=kp1; j
128             {
129 21           ipj = nip + j;
130 21           A[ipj] = A[ipj] + em * A[nkp + j];
131             }
132             }
133             }
134 4           kpn = n * IPS[n-1] + n - 1; /* last element of IPS[n] th row */
135 4 50         if( A[kpn] == 0.0 )
136             {
137 0           puts( "SIMQ A[kpn]=0");
138 0           return(3);
139             }
140            
141             /* simq 4 */
142             /* back substitution */
143              
144 4           solve:
145 4           ip = IPS[0];
146 4           X[0] = B[ip];
147 11 100         for( i=1; i
148             {
149 7           ip = IPS[i];
150 7           ipj = n * ip;
151 7           sum = 0.0;
152 18 100         for( j=0; j
153             {
154 11           sum += A[ipj] * X[j];
155 11           ++ipj;
156             }
157 7           X[i] = B[ip] - sum;
158             }
159              
160 4           ipn = n * IPS[n-1] + n - 1;
161 4           X[n-1] = X[n-1]/A[ipn];
162              
163 11 100         for( iback=1; iback
164             {
165             /* i goes (n-1),...,1 */
166 7           i = nm1 - iback;
167 7           ip = IPS[i];
168 7           nip = n*ip;
169 7           sum = 0.0;
170 18 100         for( j=i+1; j
171 11           sum += A[nip+j] * X[j];
172 7           X[i] = (X[i] - sum)/A[nip+i];
173             }
174 4           return(0);
175             }