C------------------------------------------------------------- ************ C MXM C ************ SUBROUTINE MXM(A,NAR,B,NAC,C,NBC) C C Routine on VAX to emulate CRAY SCILIB routine, MXM p. 4-22 C LIBRARY manual, implemented in VAX single precision. C INTEGER NAR,NAC,NBC,I,J,K DOUBLE PRECISION A(NAR,1),B(NAC,1),C(NAR,1) DO 30 J=1,NBC DO 20 I = 1,NAR C(I,J) = 0.0 DO 10 K = 1,NAC C(I,J) = C(I,J) + A(I,K)*B(K,J) 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN END C------------------------------------------------------------- ************ C MXV C ************ SUBROUTINE MXV(A,NAR,B,NBR,C) C C Computes a matrix, A, times a vector, B, and assumes a skip C distance between elements of the matrix to be 1. C INTEGER NAR,NBR,I,K DOUBLE PRECISION A(NAR,1),B(1),C(1) DO 20 I = 1,NAR C(I) = 0.0 DO 10 K = 1,NBR 10 C(I) = C(I) + A(I,K)*B(K) 20 CONTINUE RETURN END C------------------------------------------------------------- ************ C MXMA C ************ SUBROUTINE MXMA(A,NA,IAD,B,NB,IBD,C,NC,ICD,NAR,NAC,NBC) C C Computes a matrix, A, times a matrix, B, and allows for C arbitrary spacing of matrix elements. C Martin J. McBride. 8/21/85. C General Electric CRD, Information System Operation. C INTEGER NA,IAD,NB,IBD,NC,ICD,NAR,NAC,NBC INTEGER AC,AR,BC,BR,CC,CR,I,J,K DOUBLE PRECISION A(1,1),B(1,1),C(1,1) BC = 1 CC = 1 DO 30 J = 1,NBC CR = 1 AR = 1 DO 20 I = 1,NAR C(CR,CC) = 0.0 AC = 1 BR = 1 DO 10 K = 1,NAC C(CR,CC) = C(CR,CC) + A(AR,AC)*B(BR,BC) AC = AC + IAD BR = BR + NB 10 CONTINUE CR = CR + NC AR = AR + NA 20 CONTINUE CC = CC + ICD BC = BC + IBD 30 CONTINUE RETURN END C------------------------------------------------------------- ************ C MXVA C ************ SUBROUTINE MXVA(A,NA,IAD,B,NB,C,NC,NAR,NBR) C C Computes a matrix, A, times a vector, B, and allows for C arbitrary spacing of matrix elements. C Martin J. McBride. 8/21/85. C General Electric CRD, Information System Operation. C INTEGER NA,IAD,NB,NC,NAR,NBR INTEGER AC,AR,IB,IC,I,K DOUBLE PRECISION A(1,1),B(1),C(1) IC = 1 AR = 1 DO 20 I = 1,NAR C(IC) = 0.0 AC = 1 IB = 1 DO 10 K = 1,NBR C(IC) = C(IC) + A(AR,AC)*B(IB) AC = AC + IAD IB = IB + NB 10 CONTINUE IC = IC + NC AR = AR + NA 20 CONTINUE RETURN END C------------------------------------------------------------- ************ C MINV C ************ SUBROUTINE MINV(AB,N,ND,SCRATCH,DET,EPS,M,MODE) C C A subroutine that calculates the determinant and inverse of C a matrix, as well as solving systems of linear equations. C Martin J. McBride. 11/25/85. C General Electric CRD, Information System Operation. C INTEGER N,ND,M,MODE,OUTER,ROW,COL,I,SCOL,SROW,PIVCNT DOUBLE PRECISION AB(ND,1),SCRATCH(1),DET,EPS,MULT,COLNUM,TEMP C Initialize scratch space, with 1 to N holding the diagonal of the identity C matrix used to compute the inverse and N+1 to 2N holding the positions of C the first N columns of the matrix (for use when pivot occurs). DO 5 I = 1,N 5 SCRATCH(I) = 1.0 COLNUM = 1.0 DO 6 I = N+1,2*N SCRATCH(I) = COLNUM COLNUM = COLNUM + 1.0 6 CONTINUE C Make left, square matrix an upper triangular matrix. DET = 0.0 PIVCNT = 0 DO 10 OUTER = 1,N-1 IF (ABS(AB(OUTER,OUTER)) .LE. EPS) THEN CALL PIVOT(AB,N,ND,OUTER,SCRATCH,EPS) IF (AB(OUTER,OUTER) .EQ. 0.0) THEN PRINT* PRINT*,'*************************************' PRINT*,' MINV called with singular matrix.' PRINT*,'*************************************' PRINT* STOP ENDIF PIVCNT = PIVCNT + 1 ENDIF DO 20 ROW = OUTER+1,N MULT = AB(ROW,OUTER)/AB(OUTER,OUTER) DO 30 COL = OUTER,N+M 30 AB(ROW,COL) = AB(ROW,COL) - AB(OUTER,COL)*MULT DO 25 SCOL = 1,OUTER-1 25 AB(ROW,SCOL) = AB(ROW,SCOL) - AB(OUTER,SCOL)*MULT AB(ROW,OUTER) = AB(ROW,OUTER) - SCRATCH(OUTER)*MULT 20 CONTINUE 10 CONTINUE C Compute determinant. DET = AB(1,1) DO 40 I = 2,N 40 DET = DET*AB(I,I) DET = (-1.0)**PIVCNT * DET C Return if inverse is not to be found and there are no systems of equations C to solve. IF (MODE .EQ. 0 .AND. M .EQ. 0) RETURN C Place ones in diagonal of square matrix A. DO 80 ROW = 1,N DIV = AB(ROW,ROW) DO 90 COL = 1,N+M AB(ROW,COL) = AB(ROW,COL)/DIV 90 CONTINUE SCRATCH(ROW) = SCRATCH(ROW)/DIV 80 CONTINUE C Reduce upper triangle to zeros to give matrix A = I. DO 50 OUTER = 2,N DO 60 ROW = OUTER-1,1,-1 MULT = AB(ROW,OUTER)/AB(OUTER,OUTER) DO 70 COL = OUTER,N+M 70 AB(ROW,COL) = AB(ROW,COL) - AB(OUTER,COL)*MULT DO 65 SCOL = 1,ROW-1 65 AB(ROW,SCOL) = AB(ROW,SCOL) - AB(OUTER,SCOL)*MULT SCRATCH(ROW) = SCRATCH(ROW) - AB(OUTER,ROW)*MULT DO 63 SCOL = ROW+1,OUTER-1 63 AB(ROW,SCOL) = AB(ROW,SCOL) - AB(OUTER,SCOL)*MULT AB(ROW,OUTER) = AB(ROW,OUTER) - SCRATCH(OUTER)*MULT 60 CONTINUE 50 CONTINUE C Move diagonals of inverse to matrix AB. DO 85 I = 1,N 85 AB(I,I) = SCRATCH(I) C If pivot was made, switch rows corresponding to the columns that were C pivoted. IF (PIVCNT .EQ. 0) RETURN ROW = 1 DO 95 I = 1,N-1 SROW = INT(SCRATCH(ROW+N)) IF (SROW .NE. ROW) THEN DO 92 COL = 1,N+M TEMP = AB(ROW,COL) AB(ROW,COL) = AB(SROW,COL) AB(SROW,COL) = TEMP 92 CONTINUE TEMP = SCRATCH(ROW+N) SCRATCH(ROW+N) = SCRATCH(SROW+N) SCRATCH(SROW+N) = TEMP ELSE ROW = ROW + 1 ENDIF 95 CONTINUE RETURN END C-------------------------------------------------------------- SUBROUTINE PIVOT(AB,N,ND,OUTER,SCRATCH,EPS) C C This subroutine switches two columns of a matrix to get C a nonzero entry in the diagonal. C Martin J. McBride. 12/04/85. C General Electric CRD, Information System Operation. C INTEGER N,ND,COL,OUTER,I DOUBLE PRECISION AB(ND,1),SCRATCH(1),TEMP,EPS C Get first column with non-zero element in row OUTER. COL = OUTER + 1 10 IF (COL .GT. N) GO TO 90 IF (ABS(AB(OUTER,COL)) .GT. EPS) GO TO 20 COL = COL + 1 GO TO 10 C Switch column OUTER with column COL, which has non-zero element in C row OUTER. 20 DO 30 I = 1,N TEMP = AB(I,OUTER) AB(I,OUTER) = AB(I,COL) AB(I,COL) = TEMP 30 CONTINUE TEMP = SCRATCH(N+OUTER) SCRATCH(N+OUTER) = SCRATCH(N+COL) SCRATCH(N+COL) = TEMP 90 CONTINUE RETURN END