C
C SUBROUTINE 'SINV' OF THE IBM SCIENTIFIC SUBROUTINE PACKAGE. C C NOTE: TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS C (1) HAVE BEEN CHANGED TO (*). C C .................................................................. C C SUBROUTINE SINV C C PURPOSE C INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX C C USAGE C CALL SINV(A,N,EPS,IER) C C DESCRIPTION OF PARAMETERS C A - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC C POSITIVE DEFINITE N BY N COEFFICIENT MATRIX. C ON RETURN A CONTAINS THE RESULTANT UPPER C TRIANGULAR MATRIX. C N - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX. 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 WRONG INPUT PARAME- C TER N OR BECAUSE SOME RADICAND IS NON- C POSITIVE (MATRIX A IS NOT POSITIVE C DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI- C FICANCE) C IER=K - WARNING WHICH INDICATES LOSS OF SIGNIFI- C CANCE. THE RADICAND FORMED AT FACTORIZA- C TION STEP K+1 WAS STILL POSITIVE BUT NO C LONGER GREATER THAN ABS(EPS*A(K+1,K+1)). C C REMARKS C THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE C STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS. C IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU- C LAR MATRIX IS STORED COLUMNWISE TOO. C THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL C CALCULATED RADICANDS ARE POSITIVE. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C MFSD C C METHOD C SOLUTION IS DONE USING THE FACTORIZATION BY SUBROUTINE MFSD. C C .................................................................. C SUBROUTINE SINV(A,N,EPS,IER) C C DIMENSION A(*) DOUBLE PRECISION DIN,WORK C C FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE MFSD C A = TRANSPOSE(T) * T CALL MFSD(A,N,EPS,IER) IF(IER) 9,1,1 C C INVERT UPPER TRIANGULAR MATRIX T C PREPARE INVERSION-LOOP 1 IPIV=N*(N+1)/2 IND=IPIV C C INITIALIZE INVERSION-LOOP DO 6 I=1,N DIN=1.D0/DBLE(A(IPIV)) A(IPIV)=DIN MIN=N KEND=I-1 LANF=N-KEND IF(KEND) 5,5,2 2 J=IND C C INITIALIZE ROW-LOOP DO 4 K=1,KEND WORK=0.D0 MIN=MIN-1 LHOR=IPIV LVER=J C C START INNER LOOP DO 3 L=LANF,MIN LVER=LVER+1 LHOR=LHOR+L 3 WORK=WORK+DBLE(A(LVER)*A(LHOR)) C END OF INNER LOOP C A(J)=-WORK*DIN 4 J=J-MIN C END OF ROW-LOOP C 5 IPIV=IPIV-MIN 6 IND=IND-1 C END OF INVERSION-LOOP C C CALCULATE INVERSE(A) BY MEANS OF INVERSE(T) C INVERSE(A) = INVERSE(T) * TRANSPOSE(INVERSE(T)) C INITIALIZE MULTIPLICATION-LOOP DO 8 I=1,N IPIV=IPIV+I J=IPIV C C INITIALIZE ROW-LOOP DO 8 K=I,N WORK=0.D0 LHOR=J C C START INNER LOOP DO 7 L=K,N LVER=LHOR+K-I WORK=WORK+DBLE(A(LHOR)*A(LVER)) 7 LHOR=LHOR+L C END OF INNER LOOP C A(J)=WORK 8 J=J+K C END OF ROW- AND MULTIPLICATION-LOOP C 9 RETURN END C C======================================================================= C