SUBROUTINE MATINV (ARRAY, NORDER, DET) C C************************************************************************ C * C Inverts the symmetric double precision matrix ARRAY of degree * C NORDER and calculates its determinant DET * C * C************************************************************************ C * C DESCRIPTION OF PARAMETERS * C ARRAY - Input matrix which is replaced by its inverse. * C Note: original matrix is lost. * C NORDER - Degree of matrix (order of determinant). * C DET - Determinant of input matrix. * C * C RESTRICTIONS * C To use this routine for matrices of degree > 20, the * C of the variables IK and JK below must be increased to * C the desired degree. * C * C PROGRAMMER * C Robert Walraven * C Department of Applied Science * C University of California * C Davis, CA 95616 * C (916) 752-0360 * C * C************************************************************************ C INTEGER*2 IK(20), JK(20) DOUBLE PRECISION ARRAY(1), AMAX, SAVE C C======================================================================= C DET = 1. C C------ Main loop begins here ------------------------------------------ C DO 130 K=1, NORDER AMAX = 0D0 KK = (K-1)*NORDER C C------ Find the largest element ARRAY(I,J) in the rest of the matrix -- C 10 CONTINUE C DO 30 I=K, NORDER DO 20 J=K, NORDER JJ = (J-1)*NORDER IF (DABS(AMAX) .GT. DABS(ARRAY(I+JJ))) GO TO 20 AMAX = ARRAY(I+JJ) IK(K) = I JK(K) = J 20 CONTINUE 30 CONTINUE C C------ Interchange the rows and columns to put AMAX in ARRAY(K,K) ----- C IF (AMAX .NE. 0D0) GO TO 40 DET = 0. RETURN C 40 I = IK(K) IF (I .LT. K) GO TO 10 IF (I .EQ. K) GO TO 60 C DO 50 J=1, NORDER JJ = (J-1)*NORDER SAVE = ARRAY(K+JJ) ARRAY(K+JJ) = ARRAY(I+JJ) ARRAY(I+JJ) = -SAVE 50 CONTINUE C 60 J = JK(K) IF (J .LT. K) GO TO 10 IF (J .EQ. K) GO TO 80 C JJ = (J-1)*NORDER DO 70 I=1, NORDER SAVE = ARRAY(I+KK) ARRAY(I+KK) = ARRAY(I+JJ) ARRAY(I+JJ) = -SAVE 70 CONTINUE C 80 CONTINUE C C------ Accumulate the elements of the inverse matrix in ARRAY --------- C DO 90 I=1, NORDER IF (I .NE. K) ARRAY(I+KK) = -ARRAY(I+KK) / AMAX 90 CONTINUE C DO 110 I=1, NORDER DO 100 J=1, NORDER IF (I.EQ.K .OR. J.EQ.K) GO TO 100 JJ = (J-1)*NORDER ARRAY(I+JJ) = ARRAY(I+JJ) 1 + ARRAY(I+KK)*ARRAY(K+JJ) 100 CONTINUE 110 CONTINUE C DO 120 J=1, NORDER JJ = (J-1)*NORDER IF (J .NE. K) ARRAY(K+JJ) = ARRAY(K+JJ) / AMAX 120 CONTINUE C ARRAY(K+KK) = 1./AMAX DET = DET * AMAX 130 CONTINUE C C====== End of main loop =============================================== C C------ Restore the ordering of the matrix ----------------------------- C DO 170 L=1, NORDER K = NORDER - L + 1 KK = (K-1)*NORDER J = IK(K) IF (J .LE. K) GO TO 150 C JJ = (J-1)*NORDER DO 140 I=1, NORDER SAVE = ARRAY(I+KK) ARRAY(I+KK) = -ARRAY(I+JJ) ARRAY(I+JJ) = SAVE 140 CONTINUE C 150 CONTINUE C I = JK(K) IF (I .LE. K) GO TO 170 C DO 160 J=1, NORDER JJ = (J-1)*NORDER SAVE = ARRAY(K+JJ) ARRAY(K+JJ) = -ARRAY(I+JJ) ARRAY(I+JJ) = SAVE 160 CONTINUE C 170 CONTINUE C C----------------------------------------------------------------------- C RETURN END