C
C SUBROUTINE 'MFSD' 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 MFSD C C PURPOSE C FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX C C USAGE C CALL MFSD(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 THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE C SQUARE-ROOT OF THE DETERMINANT OF THE GIVEN MATRIX. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY. C THE GIVEN MATRIX IS REPRESENTED AS PRODUCT OF TWO TRIANGULAR C MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF C THE RETURNED RIGHT HAND FACTOR. C C .................................................................. C SUBROUTINE MFSD(A,N,EPS,IER) C C DIMENSION A(*) DOUBLE PRECISION DPIV,DSUM C C TEST ON WRONG INPUT PARAMETER N IF(N-1) 12,1,1 1 IER=0 C C INITIALIZE DIAGONAL-LOOP KPIV=0 DO 11 K=1,N KPIV=KPIV+K IND=KPIV LEND=K-1 C C CALCULATE TOLERANCE TOL=ABS(EPS*A(KPIV)) C C START FACTORIZATION-LOOP OVER K-TH ROW DO 11 I=K,N DSUM=0.D0 IF(LEND) 2,4,2 C C START INNER LOOP 2 DO 3 L=1,LEND LANF=KPIV-L LIND=IND-L 3 DSUM=DSUM+DBLE(A(LANF)*A(LIND)) C END OF INNER LOOP C C TRANSFORM ELEMENT A(IND) 4 DSUM=DBLE(A(IND))-DSUM IF(I-K) 10,5,10 C C TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE 5 IF(SNGL(DSUM)-TOL) 6,6,9 6 IF(DSUM) 12,12,7 7 IF(IER) 8,8,9 8 IER=K-1 C C COMPUTE PIVOT ELEMENT 9 DPIV=DSQRT(DSUM) A(KPIV)=DPIV DPIV=1.D0/DPIV GO TO 11 C C CALCULATE TERMS IN ROW 10 A(IND)=DSUM*DPIV 11 IND=IND+I C C END OF DIAGONAL-LOOP RETURN 12 IER=-1 RETURN END C C======================================================================= C