C                                               *****************
C                                               *    MAT1.FOR   *
C                                               *  PURSER 1994  *
C                                               *****************
C   General matrix routines
C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE MULMM                                           C
C  This, and subsequent routines, perform basic algebraic operations on real   C
C  matrices. The task performed by each routine and entry is essentially       C
C  encoded in each routine's name; the first three letters describe the        C
C  operation, the remainder defining the type of operand. This is the "code":  C
C                                                                              C
C  OPERATIONS:                                                                 C
C   ADD     add first two operands, return result as third argument            C
C   CHL     perform cholesky decomposition                                     C
C   CON     copy negative of first operand to second argument                  C
C   COP     copy first operand to second argument                              C
C   DET     evaluate log-determinant                                           C
C   DIV     divide                                                             C
C   DOT     dot product                                                        C
C   INV     invert the matrix, or linear system involving the matrix operand   C
C   MAD     multiply first two operands, but then add result to third          C
C   MUL     multiply first two operands, return result as third argument       C
C   MSB     multiply first two operands, but then subtract result from third   C
C   SUB     subtract first two operands, return result as third argument       C
C   WRT     write out                                                          C
C   ZER     set to zero                                                        C
C           A prefix of "C" on some of the later routines indicates complex    C
C           arguments.                                                         C
C                                                                              C
C  OPERANDS:                                                                   C
C   B       banded matrix                                                      C
C   D       diagonal matrix                                                    C
C   L       lower triangular matrix                                            C
C   M       matrix                                                             C
C   S       scalar                                                             C
C   T       transpose of the matrix                                            C
C   U       upper triangular matrix                                            C
C   V       vector                                                             C
C                                                                              C
C   For matrix routines, the arguments, MI, MJ, etc usually denote the         C
C   sequence of algebraic dimensions of the matrices/vectors used in           C
C   left-to-right order. Thus, for general matrix multiplication, the          C
C   first, MI, is the number of active rows of the first factor, MJ the        C
C   number of active columns of the first factor and the number of active      C
C   rows of the second factor, MK the number of active columns of the          C
C   second factor (and of the product). The arguments, NA, NB, etc usually     C
C   denote the fortran first dimensions of each matrix (but not vector)        C
C   arguments of the routine in the order in which they appear. This allows    C
C   for the fortran dimension of the matrices to exceed the algebraic          C
C   or "logical" dimensions. If it is known in advance that all matrices       C
C   are square with algebraic dimensions identical to the fortran dimensions   C
C   a convenient alternative to these general routine is the set of            C
C   "quick format" routines prefixed by "Q" in which the single dimension N    C
C   replaces the MI,MJ,..., NA,NB,... used in the general routines. Quick      C
C   format routines are grouped in module MAT3.FOR.                            C
C                                                                              C
C   A special collection of various "banded matrix" routine is found in the    C
C   later modules of this series (MAT6 MAT9 MAT10)                             C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE MULMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      LOGICAL OMUL
      DIMENSION A(NA,*),B(NB,*),C(NC,*)
      OMUL=.TRUE.
      GOTO 320
      ENTRY MADMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
320   DO 200 I=1,MI
      DO 200 K=1,MK
      IF(OMUL)C(I,K)=0.
      DO 200 J=1,MJ
      C(I,K)=C(I,K)+A(I,J)*B(J,K)
  200 continue
      RETURN
      ENTRY MULMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.TRUE.
      GOTO 321
      ENTRY MADMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
321   DO 201 I=1,MI
      DO 201 K=1,MK
      IF(OMUL)C(I,K)=0.
      DO 201 J=1,MJ
  201 C(I,K)=C(I,K)+A(I,J)*B(K,J)
      RETURN
      ENTRY MULTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.TRUE.
      GOTO 322
      ENTRY MADTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
322   DO 202 I=1,MI
      DO 202 K=1,MK
      IF(OMUL)C(I,K)=0.
      DO 202 J=1,MJ
  202 C(I,K)=C(I,K)+A(J,I)*B(J,K)
      RETURN
      ENTRY MULTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.TRUE.
      GOTO 323
      ENTRY MADTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      OMUL=.FALSE.
323   DO 203 I=1,MI
      DO 203 K=1,MK
      IF(OMUL)C(I,K)=0.
      DO 203 J=1,MJ
  203 C(I,K)=C(I,K)+A(J,I)*B(K,J)
      RETURN
      ENTRY MSBMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO 206 I=1,MI
      DO 206 K=1,MK
      DO 206 J=1,MJ
206   C(I,K)=C(I,K)-A(I,J)*B(J,K)
      RETURN
      ENTRY MSBMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO 207 I=1,MI
      DO 207 K=1,MK
      DO 207 J=1,MJ
207   C(I,K)=C(I,K)-A(I,J)*B(K,J)
      RETURN
      ENTRY MSBTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO 208 I=1,MI
      DO 208 K=1,MK
      DO 208 J=1,MJ
208   C(I,K)=C(I,K)-A(J,I)*B(J,K)
      RETURN
      ENTRY MSBTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      DO 209 I=1,MI
      DO 209 K=1,MK
      DO 209 J=1,MJ
209   C(I,K)=C(I,K)-A(J,I)*B(K,J)
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE ADDMM etc                                       C
C   Add and subtract matrices                                                  C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE ADDMM(A,B,C,MI,MJ,NA,NB,NC)
      DIMENSION A(NA,*),B(NB,*),C(NC,*)
      DO 210 I=1,MI
      DO 210 J=1,MJ
  210 C(I,J)=A(I,J)+B(I,J)
      RETURN
      ENTRY ADDMT(A,B,C,MI,MJ,NA,NB,NC)
      DO 211 I=1,MI
      DO 211 J=1,MJ
  211 C(I,J)=A(I,J)+B(J,I)
      RETURN
      ENTRY ADDTM(A,B,C,MI,MJ,NA,NB,NC)
      DO 212 I=1,MI
      DO 212 J=1,MJ
  212 C(I,J)=A(J,I)+B(I,J)
      RETURN
      ENTRY ADDTT(A,B,C,MI,MJ,NA,NB,NC)
      DO 213 I=1,MI
      DO 213 J=1,MJ
  213 C(I,J)=A(J,I)+B(J,I)
      RETURN
      ENTRY SUBMM(A,B,C,MI,MJ,NA,NB,NC)
      DO 220 I=1,MI
      DO 220 J=1,MJ
  220 C(I,J)=A(I,J)-B(I,J)
      RETURN
      ENTRY SUBMT(A,B,C,MI,MJ,NA,NB,NC)
      DO 221 I=1,MI
      DO 221 J=1,MJ
  221 C(I,J)=A(I,J)-B(J,I)
      RETURN
      ENTRY SUBTM(A,B,C,MI,MJ,NA,NB,NC)
      DO 222 I=1,MI
      DO 222 J=1,MJ
  222 C(I,J)=A(J,I)-B(I,J)
      RETURN
      ENTRY SUBTT(A,B,C,MI,MJ,NA,NB,NC)
      DO 223 I=1,MI
      DO 223 J=1,MJ
  223 C(I,J)=A(J,I)-B(J,I)
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE CHLML    etc                                    C
C  standard Cholesky decomposition of A to put lower triangular factor         C
C  in B. (A is not changed.) ODIAG flags the detection of ill-posed            C
C  input.                                                                      C
C------------------------------------------------------------------------------C
      SUBROUTINE CHLML(A,B,MI,NA,NB,ODIAG)
      LOGICAL ODIAG
      DIMENSION A(NA,*),B(NB,*)
      ODIAG=.FALSE.
      DO J=1,MI
       JM=J-1
       JP=J+1
       S=A(J,J)
       DO K=1,JM
        S=S-B(J,K)**2
       ENDDO
       IF(S.LT.0.)THEN
        WRITE(6,6000)J,S
6000    FORMAT(1X,'CHLML NEG DIAG:',1X,I5,1X,E13.6)
        S=-S
        ODIAG=.TRUE.
       ENDIF
       B(J,J)=SQRT(S)
       BJJI=1./B(J,J)
       DO I=JP,MI
        S=A(I,J)
        DO K=1,JM
         S=S-B(I,K)*B(J,K)
        ENDDO
        B(I,J)=S*BJJI
       ENDDO
       DO I=1,JM
        B(I,J)=0.
       ENDDO
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE COPM etc                                        C
C  Mainly routines to copy matrices                                            C
C------------------------------------------------------------------------------C
      SUBROUTINE COPM(A,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*)
      ENTRY EQMM(A,B,MI,MJ,NA,NB)
      DO 230 I=1,MI
      DO 230 J=1,MJ
  230 B(I,J)=A(I,J)
      RETURN
      ENTRY NEQMM(A,B,MI,MJ,NA,NB)
      ENTRY CONM(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
         B(I,J)=-A(I,J)
      ENDDO
      ENDDO
      RETURN
      ENTRY MULMS(A,SS,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
         B(I,J)=A(I,J)*SS
      ENDDO
      ENDDO
      RETURN
      ENTRY EQTM(A,B,MI,MJ,NA,NB)
      ENTRY COPT(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
       B(I,J)=A(J,I)
      ENDDO
      ENDDO
      RETURN
      ENTRY CONT(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
       B(I,J)=-A(J,I)
      ENDDO
      ENDDO
      RETURN
      ENTRY ZERM(A,MI,MJ,NA)
      DO I=1,MI
      DO J=1,MJ
       A(I,J)=0.
      ENDDO
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE INVL                                            C
C     Invert lower triangular matrix, possibly in place if A and B are same    C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE INVL(A,B,MI,NA,NB)
      DIMENSION A(NA,*),B(NB,*)
      DO 232 J=MI,1,-1
      JM=J-1
      JP=J+1
      DO 233 I=1,JM
233   B(I,J)=0.
      B(J,J)=1./A(J,J)
      DO 234 I=JP,MI
      IM=I-1
      S=0.
      DO 235 K=J,IM
235   S=S+A(I,K)*B(K,J)
234   B(I,J)=-B(I,I)*S
232   CONTINUE
      RETURN
      END


C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE INVLV                                           C
C     Solve linear system involving lower triangular (INVLV) or upper          C
C     triangular (INVUV) matrix, right-hand-side vector U and output vector V  C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE INVLV(A,U,V,MI,NA)
      DIMENSION A(NA,*),U(*),V(*)
      DO I=1,MI
       S=U(I)
       DO J=1,I-1
        S=S-A(I,J)*V(J)
       ENDDO
       V(I)=S/A(I,I)
      ENDDO
      RETURN
      ENTRY INVUV(A,U,V,MI,NA)
      DO J=MI,1,-1
       S=U(J)
       DO I=J+1,MI
        S=S-A(I,J)*V(I)
       ENDDO
       V(J)=S/A(J,J)
      ENDDO
      RETURN
      END

      FUNCTION DOT(A,B,M)
      DIMENSION A(M),B(M)
      DOT=0.
      DO I=1,M
       DOT=DOT+A(I)*B(I)
      ENDDO
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE MULMD etc                                       C
C  Multiply matrices with diagonal matrices                                    C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE MULMD(A,D,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*),D(*)
      DO 200 I=1,MI
      DO 200 J=1,MJ
  200 B(I,J)=A(I,J)*D(J)
      RETURN
      ENTRY MULTD(A,D,B,MI,MJ,NA,NB)
      DO 201 I=1,MI
      DO 201 J=1,MJ
  201 B(I,J)=A(J,I)*D(J)
      RETURN
      ENTRY MULDM(D,A,B,MI,MJ,NA,NB)
      DO 202 I=1,MI
      DO 202 J=1,MJ
  202 B(I,J)=D(I)*A(I,J)
      RETURN
      ENTRY MULDT(D,A,B,MI,MJ,NA,NB)
      DO 203 I=1,MI
      DO 203 J=1,MJ
  203 B(I,J)=D(I)*A(J,I)
      RETURN
      END

      SUBROUTINE MADMD(A,D,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*),D(*)
      DO 200 I=1,MI
      DO 200 J=1,MJ
  200 B(I,J)=B(I,J)+A(I,J)*D(J)
      RETURN
      ENTRY MADTD(A,D,B,MI,MJ,NA,NB)
      DO 201 I=1,MI
      DO 201 J=1,MJ
  201 B(I,J)=B(I,J)+A(J,I)*D(J)
      RETURN
      ENTRY MADDM(D,A,B,MI,MJ,NA,NB)
      DO 202 I=1,MI
      DO 202 J=1,MJ
  202 B(I,J)=B(I,J)+D(I)*A(I,J)
      RETURN
      ENTRY MADDT(D,A,B,MI,MJ,NA,NB)
      DO 203 I=1,MI
      DO 203 J=1,MJ
  203 B(I,J)=B(I,J)+D(I)*A(J,I)
      RETURN
      ENTRY ADDMD(A,D,B,MI,NA,NB)
      DO I=1,MI
       DO J=1,MI
        B(I,J)=A(I,J)
       ENDDO
       B(I,I)=B(I,I)+D(I)
      ENDDO
      RETURN
      END

      SUBROUTINE MULDD(A,B,C,M)
      DIMENSION A(*),B(*),C(*)
      DO 200 I=1,M
  200 C(I)=A(I)*B(I)
      RETURN
      ENTRY ADDDD(A,B,C,M)
      DO 201 I=1,M
  201 C(I)=A(I)+B(I)
      RETURN
      ENTRY SUBDD(A,B,C,M)
      DO 202 I=1,M
  202 C(I)=A(I)-B(I)
      RETURN
      ENTRY DIVDD(A,B,C,M)
      DO 203 I=1,M
  203 C(I)=A(I)/B(I)
      RETURN
      ENTRY EQDD(A,B,M)
      ENTRY COPD(A,B,M)
      DO I=1,M
       B(I)=A(I)
      ENDDO
      RETURN
      END

      SUBROUTINE COPSM(S,A,M,NA)
      DIMENSION A(NA,*),B(NB,*),D(*)
      ENTRY EQSM(S,A,M,NA)
      DO 200 I=1,M
      DO 201 J=1,M
  201 A(I,J)=0.
  200 A(I,I)=S
      RETURN
      ENTRY EQDM(D,A,M,NA)
      ENTRY COPDM(D,A,M,NA)
      DO 202 I=1,M
      DO 203 J=1,M
  203 A(I,J)=0.
  202 A(I,I)=D(I)
      RETURN
      ENTRY SUBSM(S,A,B,MI,NA,NB)
C  SUBTRACT MATRIX A FROM SCALAR MATRIX S TO GET B
      DO J=1,MI
      DO I=1,MI
      B(I,J)=-A(I,J)
      ENDDO
      B(J,J)=B(J,J)+S
      ENDDO
      RETURN
      ENTRY ADDSM(S,A,B,M,NA,NB)
C  ADD SCALAR S TO DIAGONALS OF SQUARE MATRIX A, RESULT IN B
      DO 204 I=1,M
      DO 205 J=1,M
  205 B(I,J)=A(I,J)
  204 B(I,I)=B(I,I)+S
      RETURN
      ENTRY ADDDM(D,A,B,M,NA,NB)
C  ADD DIAGONAL D TO DIAGONALS OF SQUARE MATRIX A, RESULT IN B
      DO 206 I=1,M
      DO 207 J=1,M
  207 B(I,J)=A(I,J)
  206 B(I,I)=B(I,I)+D(I)
      RETURN
      ENTRY EQLD(A,D,M,NA)
      ENTRY COPLD(A,D,M,NA)
      DO I=1,M
      D(I)=A(I,I)
      ENDDO
      RETURN
      ENTRY DETL(A,DET,MI,NA)
C  GET LOG OF DETERMINANT OF LOWER TRIANGULAR MATRIX
      DET=0.
      DO 208 I=1,MI
208   DET=DET+ALOG(A(I,I))
      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
C                   SUBROUTINE WRTM, WRTIM                                     C
C  Write out contents of a real (WRTM) or integer (WRTIM) matrix               C
C                                                                              C
C  --> A or L   matrix                                                         C
C  --> MI       number of active rows of matrix                                C
C  --> MJ       number of active columns of matrix                             C
C  --> JCOLS    number of columns output at a time                             C
C  --> N        first fortran dimension of matrix                              C
C                                                                              C
C------------------------------------------------------------------------------C
      SUBROUTINE WRTM(A,MI,MJ,JCOLS,N)
C  WRITE A REAL ARRAY
      DIMENSION A(N,*),L(N,*)
      DO 200 J1=1,MJ,JCOLS
      J2=MIN(MJ,J1+JCOLS-1)
      WRITE(6,600)(J,J=J1,J2)
      WRITE(6,604)('________',J=J1,J2)
      DO 202 I=1,MI
  202 WRITE(6,64)I,(A(I,J),J=J1,J2)
      WRITE(6,60)
  200 CONTINUE
      RETURN
      ENTRY WRTIM(L,MI,MJ,JCOLS,N)
C  WRITE AN INTEGER ARRAY
      DO 203 J1=1,MJ,JCOLS
      J2=MIN(MJ,J1+JCOLS-1)
      WRITE(6,601)(J,J=J1,J2)
      WRITE(6,604)('________',J=J1,J2)
      DO 204 I=1,MI
204   WRITE(6,61)I,(L(I,J),J=J1,J2)
      WRITE(6,60)
203   CONTINUE
60    FORMAT(/)
61    FORMAT(1X,I3,'³',8(1X,I5))
64    FORMAT(1X,I3,10(1X,E12.5))
600   FORMAT(5X,10(I3,9X))
601   FORMAT(9X,20(I3,3X))
602   FORMAT(1X,12A8)
604   FORMAT(1X,'________',11A8)
      RETURN
      END
