SUBROUTINE SPIN(Q,N,MM,DET,IDEXP) C C SUBROUTINE SPIN IS A MATRIX INVERSION ROUTINE FOR SYMMETRIC C POSITIVE-DEFINITE MATRICES. THE MATRIX INVERTED IS THE UPPER C N BY N PORTION OF THE MATRIX Q WHICH IS DIMENSIONED MM BY MM C IN THE CALLING ROUTINE. C C INPUT: C Q - THE MATRIX DIMENSIONED MM BY MM WHICH CONTAINS THE C MATRIX TO BE INVERTED. C C N - THE DIMENSION OF THE ACTUAL PART( UPPER LEFT CORNER) C OF Q WHICH IS TO BE INVERTED. ( N MAY BE EQUAL BUT MUST C NOT BE LARGER THAN MM) . C MM- DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE. C C C OUTPUT: C C Q - THE UPPER LEFT N BY N PORTION CONTAINS THE INVERSE OF C THE INPUT UPPER LEFT N BY N PORTION. C C DET - THE NON-EXPONENT PORTION OF THE DETERMINANT OF THE C INPUT N BY N (UPPER LEFT PORTION OF Q) MATRIX. SEE C IDEXP BELOW. C C IDEXP - THE EXPONENT (OF 10) PART OF THE DETERMINANT DESCRIBED C UNDER DET ABOVE. THUS THE DETERMINANT IS RETURNED IN C TWO PARTS CORRESPONDING TO C DETERMINANT = DET * 10 ** IDEXP . C THIS IS DONE TO AVOID UNDER OR OVERFLOW IN THE C COMPUTATION OF THE DETERMINANT. TO PRINT THE DETERM- C INANT THE USER SHOULD PRINT BOTH NUMBERS AS FOLLOWS; C (FOR EXAMPLE) C PRINT 10,DET,IDEXP C 10 FORMAT(' ','DETERMINANT= ',F7.4,'D',I4) C C R.R. STEEVES C SEPT., 1979 C C REAL*8 Q(MM,MM),DET,SUM,DSQRT,DABS,RPART,APART DET=0.D0 DO 4 J=1,N DO 4 I=1,J IF(I.EQ.1) GO TO 2 20 M=I-1 SUM=0.0D0 DO 1 K=1,M 1 SUM=SUM+Q(K,I)*Q(K,J) Q(I,J)=Q(I,J)-SUM 2 IF(I.EQ.J) GO TO 3 Q(I,J)=Q(I,J)/Q(I,I) GO TO 4 3 CONTINUE DET=DET+DLOG10(Q(I,I)) Q(I,I)=DSQRT(Q(I,I)) 4 CONTINUE IDEXP=DET RPART=DET-IDEXP APART=DABS(RPART) IF(APART.LT.1.D-20) DET=1.D0 IF(APART.LT.1.D-20) GO TO 10 DET=10**RPART 10 CONTINUE DO 7 J=1,N DO 7 I=1,J IF(I.LT.J) GO TO 5 Q(J,J)=1.0D0/Q(J,J) GO TO 7 5 SUM=0.0D0 M=J-1 DO 6 K=I,M 6 SUM=SUM-Q(I,K)*Q(K,J) Q(I,J)=SUM/Q(J,J) 7 CONTINUE DO 9 J=1,N DO 9 I=1,J SUM=0.0D0 DO 8 K=J,N 8 SUM=SUM+Q(I,K)*Q(J,K) Q(I,J)=SUM IF(I.EQ.J) GO TO 9 Q(J,I)=SUM 9 CONTINUE RETURN END