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