/*
C
C     ..................................................................
C
C        SUBROUTINE GELS
C
C        PURPOSE
C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH
C           SYMMETRIC COEFFICIENT MATRIX UPPER TRIANGULAR PART OF WHICH
C           IS ASSUMED TO BE STORED COLUMNWISE.
C
C        USAGE
C           CALL GELS(R,A,M,N,EPS,IER,AUX)
C
C        DESCRIPTION OF PARAMETERS
C           R      - M BY N RIGHT HAND SIDE MATRIX.  (DESTROYED)
C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
C           A      - UPPER TRIANGULAR PART OF THE SYMMETRIC
C                    M BY M COEFFICIENT MATRIX.  (DESTROYED)
C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
C                    IER=0  - NO ERROR,
C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
C                             PIVOT ELEMENT AT ANY ELIMINATION STEP
C                             EQUAL TO 0,
C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
C                             CANCE INDICATED AT ELIMINATION STEP K+1,
C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
C                             ABSOLUTELY GREATEST MAIN DIAGONAL
C                             ELEMENT OF MATRIX A.
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH DIMENSION M-1.
C
C        REMARKS
C           UPPER TRIANGULAR PART OF MATRIX A IS ASSUMED TO BE STORED
C           COLUMNWISE IN M*(M+1)/2 SUCCESSIVE STORAGE LOCATIONS, RIGHT
C           HAND SIDE MATRIX R COLUMNWISE IN N*M SUCCESSIVE STORAGE
C           LOCATIONS. ON RETURN SOLUTION MATRIX R IS STORED COLUMNWISE
C           TOO.
C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
C           GIVEN IN CASE M=1.
C           ERROR PARAMETER IER=-1 DOES NOT NECESSARILY MEAN THAT
C           MATRIX A IS SINGULAR, AS ONLY MAIN DIAGONAL ELEMENTS
C           ARE USED AS PIVOT ELEMENTS. POSSIBLY SUBROUTINE GELG (WHICH
C           WORKS WITH TOTAL PIVOTING) WOULD BE ABLE TO FIND A SOLUTION.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
C           PIVOTING IN MAIN DIAGONAL, IN ORDER TO PRESERVE
C           SYMMETRY IN REMAINING COEFFICIENT MATRICES.
C
C     ..................................................................
C
*/

#include <stdio.h>
#define fabsl(x) ( (x) < 0.0L ? -(x) : (x) )

int gels( A, R, M, EPS, AUX )
long double A[],R[];
int M;
long double EPS;
long double AUX[];
{
int I, J, K, L, IER;
int II, LL, LLD, LR, LT, LST, LLST, LEND;
long double tb, piv, tol, pivi;

IER = 0;
if( M <= 0 )
 {
fatal:
 IER = -1;
 goto done;
 }
/* SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT */

/*  Diagonal elements are at A(i,i) = 0, 2, 5, 9, 14, ...
 *  A(i,j) = A( i(i-1)/2 + j - 1 )
 */
piv = 0.0L;
I = 0;
J = 0;
L = 0;
for( K=1; K<=M; K++ )
 {
 L += K;
 tb = fabsl( A[L-1] );
 if( tb > piv )
  {
  piv = tb;
  I = L;
  J = K;
  }
 }
tol = EPS * piv;

/*
C     MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.
C     PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
*/

/*     START ELIMINATION LOOP */
LST = 0;
LEND = M - 1;
for( K=1; K<=M; K++ )
 {
/*     TEST ON USEFULNESS OF SYMMETRIC ALGORITHM */
 if( piv <= 0.0L )
   {
     printf( "gels: piv <= 0 at K = %d\n", K );
     goto fatal;
   }
 if( IER == 0 )
  {
  if( piv <= tol )
   {
   IER = K;
/*
   goto done;
*/
   }
  }
 LT = J - K;
 LST += K;

/*  PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R */
 pivi = 1.0L / A[I-1];
 L = K;
 LL = L + LT;
 tb = pivi * R[LL-1];
 R[LL-1] = R[L-1];
 R[L-1] = tb;
/* IS ELIMINATION TERMINATED */
 if( K >= M )
  break;
/*
C     ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
C     ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.
*/
 LR = LST + (LT*(K+J-1))/2;
 LL = LR;
 L=LST;
 for( II=K; II<=LEND; II++ )
  {
  L += II;
  LL += 1;
  if( L == LR )
   {
   A[LL-1] = A[LST-1];
   tb = A[L-1];
   goto lab13;
   }
  if( L > LR )
   LL = L + LT;

  tb = A[LL-1];
  A[LL-1] = A[L-1];
lab13:
  AUX[II-1] = tb;
  A[L-1] = pivi * tb;
  }
/* SAVE COLUMN INTERCHANGE INFORMATION */
 A[LST-1] = LT;
/* ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT */
 piv = 0.0L;
 LLST = LST;
 LT = 0;
 for( II=K; II<=LEND; II++ )
  {
  pivi = -AUX[II-1];
  LL = LLST;
  LT += 1;
  for( LLD=II; LLD<=LEND; LLD++ )
   {
   LL += LLD;
   L = LL + LT;
   A[L-1] += pivi * A[LL-1];
   }
  LLST += II;
  LR = LLST + LT;
  tb =fabsl( A[LR-1] );
  if( tb > piv )
   {
   piv = tb;
   I = LR;
   J = II + 1;
   }
  LR = K;
  LL = LR + LT;
  R[LL-1] += pivi * R[LR-1];
  }
 }
/* END OF ELIMINATION LOOP */

/* BACK SUBSTITUTION AND BACK INTERCHANGE */

if( LEND <= 0 )
 {
 printf( "gels: LEND = %d\n", LEND );
 if( LEND < 0 )
  goto fatal;
 goto done;
 }
II = M;
for( I=2; I<=M; I++ )
 {
 LST -= II;
 II -= 1;
 L = A[LST-1] + 0.5L;
 J = II;
 tb = R[J-1];
 LL = J;
 K = LST;
 for( LT=II; LT<=LEND; LT++ )
  {
  LL += 1;
  K += LT;
  tb -= A[K-1] * R[LL-1];
  }
 K = J + L;
 R[J-1] = R[K-1];
 R[K-1] = tb;
 }
done:
if( IER )
 printf( "gels error %d!\n", IER );
return( IER );
}