libm/double/gels.c

   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
  73 gels( A, R, M, EPS, AUX )
  74 double A[],R[];
  75 int M;
  76 double EPS;
  77 double AUX[];
  78 {
  79 int I, J, K, L, IER;
  80 int II, LL, LLD, LR, LT, LST, LLST, LEND;
  81 double tb, piv, tol, pivi;
  82
  83 if( M <= 0 )
  84         {
  85 fatal:
  86         IER = -1;
  87         goto done;
  88         }
  89 /* SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT */
  90
  91 /*  Diagonal elements are at A(i,i) = 1, 3, 6, 10, ...
  92  *  A(i,j) = A( i(i-1)/2 + j )
  93  */
  94 IER = 0;
  95 piv = 0.0;
  96 L = 0;
  97 for( K=1; K<=M; K++ )
  98         {
  99         L += K;
 100         tb = fabs( A[L-1] );
 101         if( tb > piv )
 102                 {
 103                 piv = tb;
 104                 I = L;
 105                 J = K;
 106                 }
 107         }
 108 tol = EPS * piv;
 109
 110 /*
 111 C     MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.
 112 C     PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
 113 */
 114
 115 /*     START ELIMINATION LOOP */
 116 LST = 0;
 117 LEND = M - 1;
 118 for( K=1; K<=M; K++ )
 119         {
 120 /*     TEST ON USEFULNESS OF SYMMETRIC ALGORITHM */
 121         if( piv <= 0.0 )
 122                 goto fatal;
 123         if( IER == 0 )
 124                 {
 125                 if( piv <= tol )
 126                         {
 127                         IER = K - 1;
 128                         }
 129                 }
 130         LT = J - K;
 131         LST += K;
 132
 133 /*  PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R */
 134         pivi = 1.0 / A[I-1];
 135         L = K;
 136         LL = L + LT;
 137         tb = pivi * R[LL-1];
 138         R[LL-1] = R[L-1];
 139         R[L-1] = tb;
 140 /* IS ELIMINATION TERMINATED */
 141         if( K >= M )
 142                 break;
 143 /*
 144 C     ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
 145 C     ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.
 146 */
 147         LR = LST + (LT*(K+J-1))/2;
 148         LL = LR;
 149         L=LST;
 150         for( II=K; II<=LEND; II++ )
 151                 {
 152                 L += II;
 153                 LL += 1;
 154                 if( L == LR )
 155                         {
 156                         A[LL-1] = A[LST-1];
 157                         tb = A[L-1];
 158                         goto lab13;
 159                         }
 160                 if( L > LR )
 161                         LL = L + LT;
 162
 163                 tb = A[LL-1];
 164                 A[LL-1] = A[L-1];
 165 lab13:
 166                 AUX[II-1] = tb;
 167                 A[L-1] = pivi * tb;
 168                 }
 169 /* SAVE COLUMN INTERCHANGE INFORMATION */
 170         A[LST-1] = LT;
 171 /* ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT */
 172         piv = 0.0;
 173         LLST = LST;
 174         LT = 0;
 175         for( II=K; II<=LEND; II++ )
 176                 {
 177                 pivi = -AUX[II-1];
 178                 LL = LLST;
 179                 LT += 1;
 180                 for( LLD=II; LLD<=LEND; LLD++ )
 181                         {
 182                         LL += LLD;
 183                         L = LL + LT;
 184                         A[L-1] += pivi * A[LL-1];
 185                         }
 186                 LLST += II;
 187                 LR = LLST + LT;
 188                 tb =fabs( A[LR-1] );
 189                 if( tb > piv )
 190                         {
 191                         piv = tb;
 192                         I = LR;
 193                         J = II + 1;
 194                         }
 195                 LR = K;
 196                 LL = LR + LT;
 197                 R[LL-1] += pivi * R[LR-1];
 198                 }
 199         }
 200 /* END OF ELIMINATION LOOP */
 201
 202 /* BACK SUBSTITUTION AND BACK INTERCHANGE */
 203
 204 if( LEND <= 0 )
 205         {
 206         if( LEND < 0 )
 207                 goto fatal;
 208         goto done;
 209         }
 210 II = M;
 211 for( I=2; I<=M; I++ )
 212         {
 213         LST -= II;
 214         II -= 1;
 215         L = A[LST-1] + 0.5;
 216         J = II;
 217         tb = R[J-1];
 218         LL = J;
 219         K = LST;
 220         for( LT=II; LT<=LEND; LT++ )
 221                 {
 222                 LL += 1;
 223                 K += LT;
 224                 tb -= A[K-1] * R[LL-1];
 225                 }
 226         K = J + L;
 227         R[J-1] = R[K-1];
 228         R[K-1] = tb;
 229         }
 230 done:
 231 return( IER );
 232 }