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

		SUBROUTINE MADMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      	LOGICAL OMUL
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         OMUL=.FALSE.
         DO I=1,MI
         DO K=1,MK
           IF(OMUL)C(I,K)=0.
         DO J=1,MJ
            C(I,K)=C(I,K)+A(I,J)*B(J,K)
         END DO
         END DO
         END DO
		END SUBROUTINE MADMM

      SUBROUTINE MULMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      	LOGICAL OMUL
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         OMUL=.TRUE.
         DO I=1,MI
         DO K=1,MK
           IF(OMUL)C(I,K)=0.
         DO J=1,MJ
           C(I,K)=C(I,K)+A(I,J)*B(K,J)
         END DO
         END DO
         END DO
      END SUBROUTINE MULMT

      SUBROUTINE MADMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      	LOGICAL OMUL
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         OMUL=.FALSE.
         DO I=1,MI
         DO K=1,MK
           IF(OMUL)C(I,K)=0.
         DO J=1,MJ
           C(I,K)=C(I,K)+A(I,J)*B(K,J)
         END DO
         END DO
         END DO
      END SUBROUTINE MADMT

      SUBROUTINE MULTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      	LOGICAL OMUL
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         OMUL=.TRUE.
         DO I=1,MI
         DO K=1,MK
           IF(OMUL)C(I,K)=0.
         DO J=1,MJ
           C(I,K)=C(I,K)+A(J,I)*B(J,K)
			END DO
			END DO
			END DO
      END SUBROUTINE MULTM

      SUBROUTINE MADTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      	LOGICAL OMUL
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         OMUL=.FALSE.
         DO I=1,MI
         DO K=1,MK
           IF(OMUL)C(I,K)=0.
         DO J=1,MJ
           C(I,K)=C(I,K)+A(J,I)*B(J,K)
         END DO
         END DO
         END DO
      END SUBROUTINE MADTM

      SUBROUTINE MULTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      	LOGICAL OMUL
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         OMUL=.TRUE.
         DO I=1,MI
         DO K=1,MK
           IF(OMUL)C(I,K)=0.
         DO J=1,MJ
           C(I,K)=C(I,K)+A(J,I)*B(K,J)
         END DO
         END DO
         END DO
      END SUBROUTINE MULTT

      SUBROUTINE MADTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      	LOGICAL OMUL
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         OMUL=.FALSE.
         DO I=1,MI
         DO K=1,MK
           IF(OMUL)C(I,K)=0.
         DO J=1,MJ
           C(I,K)=C(I,K)+A(J,I)*B(K,J)
			END DO
			END DO
			END DO
      END SUBROUTINE MADTT

      SUBROUTINE MSBMM(A,B,C,MI,MJ,MK,NA,NB,NC)
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         DO I=1,MI
         DO K=1,MK
         DO J=1,MJ
           C(I,K)=C(I,K)-A(I,J)*B(J,K)
         END DO
         END DO
         END DO
      END SUBROUTINE MSBMM

      SUBROUTINE MSBMT(A,B,C,MI,MJ,MK,NA,NB,NC)
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         DO I=1,MI
         DO K=1,MK
         DO J=1,MJ
           C(I,K)=C(I,K)-A(I,J)*B(K,J)
         END DO
         END DO
         END DO
      END SUBROUTINE MSBMT

      SUBROUTINE MSBTM(A,B,C,MI,MJ,MK,NA,NB,NC)
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         DO I=1,MI
         DO K=1,MK
         DO J=1,MJ
           C(I,K)=C(I,K)-A(J,I)*B(J,K)
         END DO
         END DO
         END DO
      END SUBROUTINE MSBTM

      SUBROUTINE MSBTT(A,B,C,MI,MJ,MK,NA,NB,NC)
      	DIMENSION A(NA,*),B(NB,*),C(NC,*)
         DO I=1,MI
         DO K=1,MK
         DO J=1,MJ
            C(I,K)=C(I,K)-A(J,I)*B(K,J)
         END DO
         END DO
         END DO
      END SUBROUTINE MSBTT

end module mulmm_m


!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE COPM etc                                        C
!  Mainly routines to copy matrices                                            C
!------------------------------------------------------------------------------C
      SUBROUTINE COPM(A,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*)
      ENTRY EQMM(A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=A(I,J)
		END DO
		END DO
      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
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE INVL                                            C
!     Invert lower triangular matrix, possibly in place if A and B are same    C
!                                                                              C
!------------------------------------------------------------------------------C
      SUBROUTINE INVL(A,B,MI,NA,NB)
      DIMENSION A(NA,*),B(NB,*)
      DO J=MI,1,-1
        JM=J-1
        JP=J+1
        DO I=1,JM
          B(I,J)=0.
        END DO
        B(J,J)=1./A(J,J)
        DO I=JP,MI
          IM=I-1
          S=0.
            DO K=J,IM
              S=S+A(I,K)*B(K,J)
            END DO
          B(I,J)=-B(I,I)*S
        END DO
      END DO
      END SUBROUTINE INVL

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE INVLV                                           C
!     Solve linear system involving lower triangular (INVLV) or upper          C
!     triangular (INVUV) matrix, right-hand-side vector U and output vector V  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
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE MULMD etc                                       C
!  Multiply matrices with diagonal matrices                                    C
!                                                                              C
!------------------------------------------------------------------------------C
      SUBROUTINE MULMD(A,D,B,MI,MJ,NA,NB)
      DIMENSION A(NA,*),B(NB,*),D(*)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=A(I,J)*D(J)
		END DO
		END DO
      RETURN
      ENTRY MULTD(A,D,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO  J=1,MJ
        B(I,J)=A(J,I)*D(J)
      END DO
      END DO
      RETURN
      ENTRY MULDM(D,A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=D(I)*A(I,J)
		END DO
		END DO
      RETURN
      ENTRY MULDT(D,A,B,MI,MJ,NA,NB)
      DO I=1,MI
      DO J=1,MJ
        B(I,J)=D(I)*A(J,I)
		END DO
		END DO
      RETURN
      END

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

      SUBROUTINE MULDD(A,B,C,M)
      DIMENSION A(*),B(*),C(*)
      DO I=1,M
        C(I)=A(I)*B(I)
      END DO
      RETURN
      ENTRY ADDDD(A,B,C,M)
      DO I=1,M
        C(I)=A(I)+B(I)
      END DO
      RETURN
      ENTRY SUBDD(A,B,C,M)
      DO I=1,M
        C(I)=A(I)-B(I)
      END DO
      RETURN
      ENTRY DIVDD(A,B,C,M)
      DO I=1,M
        C(I)=A(I)/B(I)
      END DO
      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 I=1,M
      DO J=1,M
        A(I,J)=0.
		END DO
        A(I,I)=S
		END DO
      RETURN
      ENTRY EQDM(D,A,M,NA)
      ENTRY COPDM(D,A,M,NA)
      DO I=1,M
      DO J=1,M
        A(I,J)=0.
		END DO
        A(I,I)=D(I)
		END DO
      RETURN
      ENTRY SUBSM(S,A,B,MI,NA,NB)
!  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)
!  ADD SCALAR S TO DIAGONALS OF SQUARE MATRIX A, RESULT IN B
      DO I=1,M
      DO J=1,M
        B(I,J)=A(I,J)
		END DO
        B(I,I)=B(I,I)+S
		END DO
      RETURN
      ENTRY ADDDM(D,A,B,M,NA,NB)
!  ADD DIAGONAL D TO DIAGONALS OF SQUARE MATRIX A, RESULT IN B
      DO I=1,M
      DO J=1,M
        B(I,J)=A(I,J)
		END DO
        B(I,I)=B(I,I)+D(I)
		END DO
      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)
!  GET LOG OF DETERMINANT OF LOWER TRIANGULAR MATRIX
      DET=0.
      DO I=1,MI
        DET=DET+ALOG(A(I,I))
      END DO
      RETURN
      END

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