gels.c 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232
  1. /*
  2. C
  3. C ..................................................................
  4. C
  5. C SUBROUTINE GELS
  6. C
  7. C PURPOSE
  8. C TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH
  9. C SYMMETRIC COEFFICIENT MATRIX UPPER TRIANGULAR PART OF WHICH
  10. C IS ASSUMED TO BE STORED COLUMNWISE.
  11. C
  12. C USAGE
  13. C CALL GELS(R,A,M,N,EPS,IER,AUX)
  14. C
  15. C DESCRIPTION OF PARAMETERS
  16. C R - M BY N RIGHT HAND SIDE MATRIX. (DESTROYED)
  17. C ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
  18. C A - UPPER TRIANGULAR PART OF THE SYMMETRIC
  19. C M BY M COEFFICIENT MATRIX. (DESTROYED)
  20. C M - THE NUMBER OF EQUATIONS IN THE SYSTEM.
  21. C N - THE NUMBER OF RIGHT HAND SIDE VECTORS.
  22. C EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
  23. C TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
  24. C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  25. C IER=0 - NO ERROR,
  26. C IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
  27. C PIVOT ELEMENT AT ANY ELIMINATION STEP
  28. C EQUAL TO 0,
  29. C IER=K - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  30. C CANCE INDICATED AT ELIMINATION STEP K+1,
  31. C WHERE PIVOT ELEMENT WAS LESS THAN OR
  32. C EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
  33. C ABSOLUTELY GREATEST MAIN DIAGONAL
  34. C ELEMENT OF MATRIX A.
  35. C AUX - AN AUXILIARY STORAGE ARRAY WITH DIMENSION M-1.
  36. C
  37. C REMARKS
  38. C UPPER TRIANGULAR PART OF MATRIX A IS ASSUMED TO BE STORED
  39. C COLUMNWISE IN M*(M+1)/2 SUCCESSIVE STORAGE LOCATIONS, RIGHT
  40. C HAND SIDE MATRIX R COLUMNWISE IN N*M SUCCESSIVE STORAGE
  41. C LOCATIONS. ON RETURN SOLUTION MATRIX R IS STORED COLUMNWISE
  42. C TOO.
  43. C THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
  44. C GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
  45. C ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
  46. C INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
  47. C SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
  48. C INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
  49. C GIVEN IN CASE M=1.
  50. C ERROR PARAMETER IER=-1 DOES NOT NECESSARILY MEAN THAT
  51. C MATRIX A IS SINGULAR, AS ONLY MAIN DIAGONAL ELEMENTS
  52. C ARE USED AS PIVOT ELEMENTS. POSSIBLY SUBROUTINE GELG (WHICH
  53. C WORKS WITH TOTAL PIVOTING) WOULD BE ABLE TO FIND A SOLUTION.
  54. C
  55. C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  56. C NONE
  57. C
  58. C METHOD
  59. C SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
  60. C PIVOTING IN MAIN DIAGONAL, IN ORDER TO PRESERVE
  61. C SYMMETRY IN REMAINING COEFFICIENT MATRICES.
  62. C
  63. C ..................................................................
  64. C
  65. */
  66. #include <math.h>
  67. #ifdef ANSIPROT
  68. extern double fabs ( double );
  69. #else
  70. double fabs();
  71. #endif
  72. gels( A, R, M, EPS, AUX )
  73. double A[],R[];
  74. int M;
  75. double EPS;
  76. double AUX[];
  77. {
  78. int I, J, K, L, IER;
  79. int II, LL, LLD, LR, LT, LST, LLST, LEND;
  80. double tb, piv, tol, pivi;
  81. if( M <= 0 )
  82. {
  83. fatal:
  84. IER = -1;
  85. goto done;
  86. }
  87. /* SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT */
  88. /* Diagonal elements are at A(i,i) = 1, 3, 6, 10, ...
  89. * A(i,j) = A( i(i-1)/2 + j )
  90. */
  91. IER = 0;
  92. piv = 0.0;
  93. L = 0;
  94. for( K=1; K<=M; K++ )
  95. {
  96. L += K;
  97. tb = fabs( A[L-1] );
  98. if( tb > piv )
  99. {
  100. piv = tb;
  101. I = L;
  102. J = K;
  103. }
  104. }
  105. tol = EPS * piv;
  106. /*
  107. C MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.
  108. C PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
  109. */
  110. /* START ELIMINATION LOOP */
  111. LST = 0;
  112. LEND = M - 1;
  113. for( K=1; K<=M; K++ )
  114. {
  115. /* TEST ON USEFULNESS OF SYMMETRIC ALGORITHM */
  116. if( piv <= 0.0 )
  117. goto fatal;
  118. if( IER == 0 )
  119. {
  120. if( piv <= tol )
  121. {
  122. IER = K - 1;
  123. }
  124. }
  125. LT = J - K;
  126. LST += K;
  127. /* PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R */
  128. pivi = 1.0 / A[I-1];
  129. L = K;
  130. LL = L + LT;
  131. tb = pivi * R[LL-1];
  132. R[LL-1] = R[L-1];
  133. R[L-1] = tb;
  134. /* IS ELIMINATION TERMINATED */
  135. if( K >= M )
  136. break;
  137. /*
  138. C ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
  139. C ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.
  140. */
  141. LR = LST + (LT*(K+J-1))/2;
  142. LL = LR;
  143. L=LST;
  144. for( II=K; II<=LEND; II++ )
  145. {
  146. L += II;
  147. LL += 1;
  148. if( L == LR )
  149. {
  150. A[LL-1] = A[LST-1];
  151. tb = A[L-1];
  152. goto lab13;
  153. }
  154. if( L > LR )
  155. LL = L + LT;
  156. tb = A[LL-1];
  157. A[LL-1] = A[L-1];
  158. lab13:
  159. AUX[II-1] = tb;
  160. A[L-1] = pivi * tb;
  161. }
  162. /* SAVE COLUMN INTERCHANGE INFORMATION */
  163. A[LST-1] = LT;
  164. /* ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT */
  165. piv = 0.0;
  166. LLST = LST;
  167. LT = 0;
  168. for( II=K; II<=LEND; II++ )
  169. {
  170. pivi = -AUX[II-1];
  171. LL = LLST;
  172. LT += 1;
  173. for( LLD=II; LLD<=LEND; LLD++ )
  174. {
  175. LL += LLD;
  176. L = LL + LT;
  177. A[L-1] += pivi * A[LL-1];
  178. }
  179. LLST += II;
  180. LR = LLST + LT;
  181. tb =fabs( A[LR-1] );
  182. if( tb > piv )
  183. {
  184. piv = tb;
  185. I = LR;
  186. J = II + 1;
  187. }
  188. LR = K;
  189. LL = LR + LT;
  190. R[LL-1] += pivi * R[LR-1];
  191. }
  192. }
  193. /* END OF ELIMINATION LOOP */
  194. /* BACK SUBSTITUTION AND BACK INTERCHANGE */
  195. if( LEND <= 0 )
  196. {
  197. if( LEND < 0 )
  198. goto fatal;
  199. goto done;
  200. }
  201. II = M;
  202. for( I=2; I<=M; I++ )
  203. {
  204. LST -= II;
  205. II -= 1;
  206. L = A[LST-1] + 0.5;
  207. J = II;
  208. tb = R[J-1];
  209. LL = J;
  210. K = LST;
  211. for( LT=II; LT<=LEND; LT++ )
  212. {
  213. LL += 1;
  214. K += LT;
  215. tb -= A[K-1] * R[LL-1];
  216. }
  217. K = J + L;
  218. R[J-1] = R[K-1];
  219. R[K-1] = tb;
  220. }
  221. done:
  222. return( IER );
  223. }