SUBROUTINE MATINV(ARRAY,N,DET) C C MATRIX INVERSION ROUTINE--BEVINGTON,P.302 C C ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED C N IS THE ORDER, DET IS THE CALCULATED DETERMINANT C THE INVERTED MATRIX REPLACES THE ORIGINAL ONE C c This routine has been optimized for virtual memory c PARAMETER NMAX=50 ! MAX NUMBER OF TERMS DIMENSION ARRAY(N,N) DIMENSION IK(NMAX),JK(NMAX) IF(N .gt. NMAX) THEN CALL MSG('*** ERROR *** -Array dimension larger than 50 in MATINV') CALL EXIT ENDIF DET = 1. d TYPE *,'Original matrix' d DO 1003 J = 1,N d1003 TYPE 1011,(ARRAY(I,J),I=1,N) DO 100 K = 1,N C C FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL C PARTIAL PIVOTING C AMAX = 0.0 21 DO 30 J = K,N DO 30 I = K,N IF(ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30 24 AMAX = ARRAY(I,J) IK(K) = I JK(K) = J 30 CONTINUE d TYPE 1010,K,AMAX,IK(K),JK(K) d1010 FORMAT(I3,'AMAX='1p,G15.6,'(',I2,',',I2,')') IF(AMAX .eq. 0.) THEN ! Matrix element too small ? DET = 0. RETURN ENDIF I = IK(K) IF(I.eq.K) go to 51 ! Current col ? DO 50 J = 1,N ! Reorder columns SAVE = ARRAY(K,J) ARRAY(K,J) = ARRAY(I,J) 50 ARRAY(I,J) = -SAVE 51 J = JK(K) IF(J.eq.K) go to 61 ! Current row ? DO 60 I = 1,N ! Reorder rows SAVE = ARRAY(I,K) ARRAY(I,K) = ARRAY(I,J) 60 ARRAY(I,J) = -SAVE 61 CONTINUE d TYPE *,'Reordered matrix' d DO 1002 J = 1,N d1002 TYPE 1011,(ARRAY(I,J),I=1,N) d1011 FORMAT(1x,1p,8G15.6) C C DIAGONALIZE MATRIX C DO 80 J = 1,N IF(J .eq. K) go to 80 ARRAY(K,J) = -ARRAY(K,J)/AMAX SAVE = ARRAY(K,J) DO 80 I = 1,N IF(I.ne.K)ARRAY(I,J) = ARRAY(I,J)+SAVE*ARRAY(I,K) 80 CONTINUE ARRAY(K,K) = 1. DO 90 I = 1,N ARRAY(I,K) = ARRAY(I,K)/AMAX 90 CONTINUE d TYPE *,'Diagonalized matrix' d DO 1004 J = 1,N d1004 TYPE 1011,(ARRAY(I,J),I=1,N) 100 DET = DET*AMAX C C RESTORE ORDERING OF MATRIX C 101 DO 130 L = 1,N K = N-L+1 J = IK(K) IF(J.le.K)go to 111 ! No reorder ? DO 110 I = 1,N SAVE = ARRAY(I,K) ARRAY(I,K) = -ARRAY(I,J) 110 ARRAY(I,J) = SAVE 111 I = JK(K) IF(I.le.K)go to 130 ! No reorder ? DO 120 J = 1,N SAVE = ARRAY(K,J) ARRAY(K,J) = -ARRAY(I,J) 120 ARRAY(I,J) = SAVE 130 CONTINUE RETURN END SUBROUTINE DMATNV(ARRAY,N,DET) C C MATRIX INVERSION ROUTINE--BEVINGTON,P.302 C C ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED C N IS THE ORDER, DET IS THE CALCULATED DETERMINANT C THE INVERTED MATRIX REPLACES THE ORIGINAL ONE C c This routine has been optimized for virtual memory c PARAMETER NMAX=50 ! MAX NUMBER OF TERMS IMPLICIT REAL*8 (A-H, O-Z) DIMENSION ARRAY(N,N) DIMENSION IK(NMAX),JK(NMAX) IF(N .gt. NMAX) THEN CALL MSG('*** ERROR *** -Array dimension larger than 50 in DMATNV') CALL EXIT ENDIF DET = 1. d TYPE *,'Original matrix' d DO 1003 J = 1,N d1003 TYPE 1011,(ARRAY(I,J),I=1,N) DO 100 K = 1,N C C FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL C PARTIAL PIVOTING C AMAX = 0.0 21 DO 30 J = K,N DO 30 I = K,N IF(ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30 24 AMAX = ARRAY(I,J) IK(K) = I JK(K) = J 30 CONTINUE d TYPE 1010,K,AMAX,IK(K),JK(K) d1010 FORMAT(I3,'AMAX='1p,G15.6,'(',I2,',',I2,')') IF(AMAX .eq. 0.) THEN ! Matrix element too small ? DET = 0. RETURN ENDIF I = IK(K) IF(I.eq.K) go to 51 ! Current col ? DO 50 J = 1,N ! Reorder columns SAVE = ARRAY(K,J) ARRAY(K,J) = ARRAY(I,J) 50 ARRAY(I,J) = -SAVE 51 J = JK(K) IF(J.eq.K) go to 61 ! Current row ? DO 60 I = 1,N ! Reorder rows SAVE = ARRAY(I,K) ARRAY(I,K) = ARRAY(I,J) 60 ARRAY(I,J) = -SAVE 61 CONTINUE d TYPE *,'Reordered matrix' d DO 1002 J = 1,N d1002 TYPE 1011,(ARRAY(I,J),I=1,N) d1011 FORMAT(1x,1p,8G15.6) C C DIAGONALIZE MATRIX C DO 80 J = 1,N IF(J .eq. K) go to 80 ARRAY(K,J) = -ARRAY(K,J)/AMAX SAVE = ARRAY(K,J) DO 80 I = 1,N IF(I.ne.K)ARRAY(I,J) = ARRAY(I,J)+SAVE*ARRAY(I,K) 80 CONTINUE ARRAY(K,K) = 1. DO 90 I = 1,N ARRAY(I,K) = ARRAY(I,K)/AMAX 90 CONTINUE d TYPE *,'Diagonalized matrix' d DO 1004 J = 1,N d1004 TYPE 1011,(ARRAY(I,J),I=1,N) 100 DET = DET*AMAX C C RESTORE ORDERING OF MATRIX C 101 DO 130 L = 1,N K = N-L+1 J = IK(K) IF(J.le.K)go to 111 ! No reorder ? DO 110 I = 1,N SAVE = ARRAY(I,K) ARRAY(I,K) = -ARRAY(I,J) 110 ARRAY(I,J) = SAVE 111 I = JK(K) IF(I.le.K)go to 130 ! No reorder ? DO 120 J = 1,N SAVE = ARRAY(K,J) ARRAY(K,J) = -ARRAY(I,J) 120 ARRAY(I,J) = SAVE 130 CONTINUE RETURN END SUBROUTINE HMATNV(ARRAY,N,DET) C C MATRIX INVERSION ROUTINE--BEVINGTON,P.302 C C ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED C N IS THE ORDER, DET IS THE CALCULATED DETERMINANT C THE INVERTED MATRIX REPLACES THE ORIGINAL ONE C c This routine has been optimized for virtual memory c PARAMETER NMAX=50 ! MAX NUMBER OF TERMS IMPLICIT REAL*16 (A-H, O-Z) DIMENSION ARRAY(N,N) DIMENSION IK(NMAX),JK(NMAX) IF(N .gt. NMAX) THEN CALL MSG('*** ERROR *** -Array dimension larger than 50 in DMATNV') CALL EXIT ENDIF DET = 1. d TYPE *,'Original matrix' d DO 1003 J = 1,N d1003 TYPE 1011,(ARRAY(I,J),I=1,N) DO 100 K = 1,N C C FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL C PARTIAL PIVOTING C AMAX = 0.0 21 DO 30 J = K,N DO 30 I = K,N IF(ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30 24 AMAX = ARRAY(I,J) IK(K) = I JK(K) = J 30 CONTINUE d TYPE 1010,K,AMAX,IK(K),JK(K) d1010 FORMAT(I3,'AMAX='1p,G15.6,'(',I2,',',I2,')') IF(AMAX .eq. 0.) THEN ! Matrix element too small ? DET = 0. RETURN ENDIF I = IK(K) IF(I.eq.K) go to 51 ! Current col ? DO 50 J = 1,N ! Reorder columns SAVE = ARRAY(K,J) ARRAY(K,J) = ARRAY(I,J) 50 ARRAY(I,J) = -SAVE 51 J = JK(K) IF(J.eq.K) go to 61 ! Current row ? DO 60 I = 1,N ! Reorder rows SAVE = ARRAY(I,K) ARRAY(I,K) = ARRAY(I,J) 60 ARRAY(I,J) = -SAVE 61 CONTINUE d TYPE *,'Reordered matrix' d DO 1002 J = 1,N d1002 TYPE 1011,(ARRAY(I,J),I=1,N) d1011 FORMAT(1x,1p,8G15.6) C C DIAGONALIZE MATRIX C DO 80 J = 1,N IF(J .eq. K) go to 80 ARRAY(K,J) = -ARRAY(K,J)/AMAX SAVE = ARRAY(K,J) DO 80 I = 1,N IF(I.ne.K)ARRAY(I,J) = ARRAY(I,J)+SAVE*ARRAY(I,K) 80 CONTINUE ARRAY(K,K) = 1. DO 90 I = 1,N ARRAY(I,K) = ARRAY(I,K)/AMAX 90 CONTINUE d TYPE *,'Diagonalized matrix' d DO 1004 J = 1,N d1004 TYPE 1011,(ARRAY(I,J),I=1,N) 100 DET = DET*AMAX C C RESTORE ORDERING OF MATRIX C 101 DO 130 L = 1,N K = N-L+1 J = IK(K) IF(J.le.K)go to 111 ! No reorder ? DO 110 I = 1,N SAVE = ARRAY(I,K) ARRAY(I,K) = -ARRAY(I,J) 110 ARRAY(I,J) = SAVE 111 I = JK(K) IF(I.le.K)go to 130 ! No reorder ? DO 120 J = 1,N SAVE = ARRAY(K,J) ARRAY(K,J) = -ARRAY(I,J) 120 ARRAY(I,J) = SAVE 130 CONTINUE RETURN END