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