C 
C Program MATMUL to compute product M3=M1*M2 of two matrices M1 and M2.
C Matrices M1 and M2 may be transposed before the multiplication.
C
C Version: 6.10
C Date: 2007, June 13
C
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Filenames of the header files of the matrices:
C     MATIN1='string' .. Name of the header file of the input matrix M1.
C             No default, MATIN1 must be specified and cannot be blank.
C     MATIN2='string' .. Name of the header file of the input matrix M2.
C             Default: MATIN2=' ' ... no input matrix M2, matrix M1 will
C             be written on output. This option may be used to transpose
C             matrix M1, or to change the form of the matrix M1.
C     MATOUT='string' . Name of the header file of the output matrix M3.
C             No default, MATOUT must be specified and cannot be blank.
C     For general description of the files with matrices refer to file
C     forms.htm.
C Form of the output file with the matrix M3:
C     SYMMETRY='string' ... Symmetry of the matrix. For SYMMETRY='diag',
C             only the diagonal components of the matrix are calculated
C             and written to the output file, other components are not
C             calculated.  Simmilarly for other symmetries.
C             Possible values of SYMMETRY are:
C               SYMMETRY='diag' ... diagonal matrix
C               SYMMETRY='sym'  ... symmetric matrix
C               SYMMETRY='skew' ... antisymmetric matrix
C               SYMMETRY=' '    ... general matrix
C             Default:
C               For matrix M2 not specified (MATIN2=' '):
C                 SYMMETRY equal to symmetry of matrix M1
C               For both input matrices M1 and M2 diagonal:
C                 SYMMETRY='diag' (output matrix M3 will be diagonal)
C               In other cases:
C                 SYMMETRY=' ' (output matrix M3 will be general)
C     MATFORM='string' ... Form of the output file with the matrix.
C             Possible values of MATFORM are:
C               MATFORM='formatted'   ... formatted file
C               MATFORM='unformatted' ... unformatted file
C             Default: MATFORM='formatted'
C     SPARSE=integer ... Identifies whether the matrix should be sparse.
C             Possible values of SPARSE are:
C               SPARSE=1  ... sparse matrix
C               SPARSE=0  ... automatic selection of the sparseness:
C                 matrix with half or more zero elements will be sparse
C               SPARSE=-1 ... normal (not sparse) matrix
C             Default: SPARSE=0 (automatic selection of the sparseness)
C     For detailed description of the forms of the files with matrices
C     refer to file forms.htm.
C Switches identifying transposition of the input matrices:
C     MATT1=integer ... Identifies whether the matrix M1 is to be
C             transposed before the multiplication:
C             MATT1=0 ... no transposition
C             MATT1=1 ... matrix M1 is to be transposed
C             Default: MATT1=0 (no transposition of M1)
C     MATT2=integer ... Same as MATT1, but for matrix M2.
C             Default: MATT2=0 (no transposition of M2)
C Optional parameter specifying the form of the output formatted
C matrix data files:
C     NUMLIN=positive integer ... Number of reals or integers in one
C             line of the output file. See the description in file
C             mat.for.
C
C For detailed description of storage of matrices in the memory
C refer to file mat.for.
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3I,RSEP3T,RMATH,RMATD,WMATH,WMATD,
     *TMATR,SMATRE,SMATR,GSMAT,GSMATR,SGMATR,ISYM,NELMAT,MSHIFT,
     *MMILOC,MMIRVE,OMULS,OMULNS
      INTEGER ISYM,NELMAT
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T,SSEP ...
C             File sep.for.
C     RMATH,RMATD,WMATH,WMATD  ...
C     TMATR,SMATRE,SMATR,GSMAT,GSMATR,SGMATR,ISYM,NELMAT,MSHIFT  ...
C             File forms.for.
C     MMILOC,MMIRVE,OMULS,OMULNS ... This file.
C Subroutines and external functions referred by above subroutines:
C     LENGTH ... File length.for.
C     WSEP3I,WSEP3T ... File sep.for.
C     RARRAI,RARRAY,WARRAI,WARRAY ...
C             File forms.for.
C     INDEXI ... File indexi.for.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
      CHARACTER*80 FILSEP,FILE1,FILE2,FILE3,FILED1,FILED2,FILED3
      CHARACTER*80 SYM1,SYM2,SYM3,SYM3I,SYMC,FORM1,FORM2,FORM3
      INTEGER MATT1,MATT2,M1,M2,M2B,M3,NEL1,NEL2,NELC,OA,OA1,OAN,NA,
     *OB,OB1,OBN,NB,OC,NC,MC,NCTMP,NCCOL,OTMP,MTMP,
     *NSPAR1,NSPAR2,ISPAR3,NSPARC,ISYM1,ISYM2,ISYM3,ISYMC,
     *ISTORA,ISTORB,ISTORC,LU1,
     *I1,II,JJ,KK,II1,II2,II3,IIA,JJ1,JJ2,JJ3,JJA,IIKK,OCTMP1
      REAL AA,BB,CIK,RSPAR1,RSPAR2,RSPAR3
      CHARACTER*72 TXTERR
      PARAMETER (LU1=1)
C
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      WRITE(*,'(A)') '+MATMUL: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
C
C     Reading all the data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU1,FILSEP)
      ELSE
C       MATMUL-01
        CALL ERROR('MATMUL-01: SEP file not given.')
      ENDIF
C
C     Reading the names of matrices header files:
      CALL RSEP3T('MATIN1',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
C       MATMUL-02
        CALL ERROR('MATMUL-02: File MATIN1 not given.')
      ENDIF
      CALL RSEP3T('MATIN2',FILE2,' ')
      CALL RSEP3T('MATOUT',FILE3,' ')
      FILED3=' '
      IF (FILE3.EQ.' ') THEN
C       MATMUL-03
        CALL ERROR('MATMUL-03: File MATOUT not given.')
      ENDIF
C     Reading the header files of the input matrices:
      CALL RMATH(LU1,FILE1,FILED1,M1,M2,NSPAR1,SYM1,FORM1,NEL1)
      ISYM1=ISYM(SYM1)
      IF (FILE2.NE.' ') THEN
        CALL RMATH(LU1,FILE2,FILED2,M2B,M3,NSPAR2,SYM2,FORM2,NEL2)
        ISYM2=ISYM(SYM2)
      ENDIF
C     Reading the properties of the output file:
      IF (FILE2.EQ.' ') THEN
        SYM3I=SYM1
      ELSEIF ((SYM1.EQ.'diag').AND.(SYM2.EQ.'diag')) THEN
        SYM3I='diag'
      ELSE
        SYM3I=' '
      ENDIF
      CALL RSEP3T('SYMMETRY',SYM3,SYM3I)
      CALL LOWER(SYM3)
      ISYM3=ISYM(SYM3)
      CALL RSEP3T('MATFORM',FORM3,'formatted')
      CALL LOWER(FORM3)
      CALL RSEP3I('SPARSE',ISPAR3,0)
C     Reading the transposition switches:
      CALL RSEP3I('MATT1',MATT1,0)
      CALL RSEP3I('MATT2',MATT2,0)
C
C     Reading input matrix A:
      IF (NSPAR1.GE.0) THEN
        NA=2*NEL1
        OA=1
        OTMP=NA+1
        MTMP=MRAM-NA
      ELSE
        NA=NEL1
        OA=MRAM-NA+1
        OTMP=1
        MTMP=OA-1
      ENDIF
      IF (NA+1.GT.MRAM) THEN
C       MATMUL-04
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATMUL-04: Array RAM too small,',NA+1-MRAM,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
      CALL RMATD(LU1,FILED1,NSPAR1,FORM1,NEL1,IRAM(OA),RAM(OA))
      ISTORA=0
C     Checking whether the non-sparse matrix is to be changed to the
C     sparse matrix, changing the matrix:
      IF (NSPAR1.LT.0) THEN
        RSPAR1=0.5
        CALL GSMAT(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM,
     *                  OA,NA,ISTORA,RSPAR1)
        IF (NSPAR1.GE.0) THEN
          OTMP=NA+1
          MTMP=MRAM-NA
        ENDIF
      ELSE
        RSPAR1=FLOAT(NSPAR1)/FLOAT(NSPAR1+NEL1)
      ENDIF
C     Transposing input matrix A:
      IF (MATT1.EQ.1) THEN
        CALL TMATR(M1,M2,SYM1,NSPAR1,NEL1,IRAM(OA),RAM(OA),
     *             IRAM(OTMP),RAM(OTMP),MTMP)
      ENDIF
C     Extending sparse matrix A:
      IF ((NSPAR1.GE.0).AND.((SYM1.EQ.'sym').OR.(SYM1.EQ.'skew'))) THEN
        CALL SMATRE
     *       (M1,M2,SYM1,NSPAR1,NEL1,IRAM(OA),RAM(OA),MRAM,NA,ISTORA)
        ISYM1=ISYM(SYM1)
      ENDIF
C     Reorganizing sparse matrix A in the memory:
      IF (NSPAR1.GE.0) THEN
        CALL SMATR(M1,M2,NSPAR1,NEL1,IRAM,RAM,MRAM,
     *             OA,NA,ISTORA)
      ENDIF
C
      IF (FILE2.EQ.' ') THEN
C       Matrix A is to be written as output matrix C:
        M3=M2
        IF (SYM3.EQ.SYM1) THEN
C         No change of matrix symmetry:
          OC=OA
          NC=NA
          NSPARC=NSPAR1
          NELC=NEL1
          SYMC=SYM1
          ISYMC=ISYM1
          ISTORC=ISTORA
C         Copying matrix A to the position of C:
          IF (NSPARC.GE.0) THEN
            DO 2, I1=1,M3+1
              IRAM(I1)=IRAM(OC-1+I1)-OC+1
   2        CONTINUE
            DO 3, I1=1,NC-M3-1-1,2
              IRAM(1+M3+I1)=IRAM(OC+M3+I1)
              RAM(1+M3+I1+1)=RAM(OC+M3+I1+1)
   3        CONTINUE
          ELSE
            DO 4, I1=1,NC
              RAM(I1)=RAM(OC-1+I1)
   4        CONTINUE
          ENDIF
          OC=1
        ELSE
C         Changing matrix symmetry:
          IF ((ISYM3.LE.3).AND.(M1.NE.M3)) THEN
C           MATMUL-05
            CALL ERROR('MATMUL-05: Wrong symmetry of output matrix.')
C           If the output matrix should have symmetry 'diag', 'sym' or
C           'skew', it must have the same number of rows and columns.
          ENDIF
          OC=1
          IF     (ISYM3.EQ.2) THEN
C           Symmetric matrix
            NC=M1*(M1+1)/2
          ELSEIF (ISYM3.EQ.3) THEN
C           Antisymmetric matrix
            NC=M1*(M1-1)/2
          ELSEIF (ISYM3.EQ.1) THEN
C           Diagonal matrix
            NC=M1
          ELSE
C           General matrix
            NC=M1*M3
          ENDIF
          MC=OA-1
          NSPARC=-1
          NELC=NC
          SYMC=SYM3
          ISYMC=ISYM3
          ISTORC=0
C
          NCTMP=0
          OCTMP1=0
          DO 15, KK=1,M3
C           Preparing sufficient number of the storage locations:
C           Number NCCOL of storage locations for column KK:
            IF     (ISYMC.EQ.1) THEN
              NCCOL=1
            ELSEIF (ISYMC.EQ.2) THEN
              NCCOL=KK
            ELSEIF (ISYMC.EQ.3) THEN
              NCCOL=KK-1
            ELSE
              NCCOL=M1
            ENDIF
            IF (NCTMP+NCCOL.GT.MC) THEN
C             Erasing unneeded columns 1 to KK-1 of A:
              IF (NSPAR1.LT.0) THEN
C               OA stays unchanged, MC is redefined:
                IF     (ISYM1.EQ.1) THEN
                  MC=OA+KK-1-1
                ELSEIF (ISYM1.EQ.2) THEN
                  MC=OA+(KK-1)*KK/2-1
                ELSEIF (ISYM1.EQ.3) THEN
                  MC=OA+(KK-2)*(KK-1)/2-1
                ELSE
                  MC=OA+M1*(KK-1)-1
                ENDIF
              ELSE
C               OA must be changed, addresses of the ends of the columns
C               must be moved:
                OAN=IRAM(OA+KK-1)-M2
                DO 12, I1=M2,0,-1
                  IF (I1.GE.KK) THEN
                    IRAM(OAN+I1)=IRAM(OA+I1)
                  ELSE
                    IRAM(OAN+I1)=IRAM(OA+KK-1)
                  ENDIF
  12            CONTINUE
                OA=OAN
                MC=OA-1
              ENDIF
              IF (NCTMP+NCCOL.GT.MC) THEN
C               Changing the written part of C to sparse matrix:
                CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OA-1,
     *                      OC,NCTMP,ISTORC,KK-1,OCTMP1)
              ENDIF
              IF (NCTMP+NCCOL.GT.MC) THEN
C               MATMUL-06
                WRITE(TXTERR,'(A,I9,A)')
     *         'MATMUL-06: Array RAM too small,',NCTMP+NCCOL-MC,
     *         ' units missing.'
                CALL ERROR(TXTERR)
              ENDIF
            ENDIF
C           Initiating the values in column KK of C:
            DO 11, I1=NCTMP+1,NCTMP+NCCOL
              RAM(I1)=0.
  11        CONTINUE
C           Preparing loop over column KK of A:
            CALL MMILOC(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,KK,II1,II2,II3)
C           Loop over nonzero elements of column KK of A:
            DO 13, IIA=II1,II2,II3
C             Computing index II of row and value AA of element of A:
              CALL MMIRVE(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,KK,IIA,II1,II3,
     *                                     II,AA)
              IF     ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN
                IIKK=OCTMP1+II
              ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN
                IIKK=OCTMP1+(KK-1)*KK/2+II
              ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN
                IIKK=OCTMP1+(KK-2)*(KK-1)/2+II
              ELSEIF (ISYMC.EQ.4)                  THEN
                IIKK=OCTMP1+(KK-1)*M1+II
              ELSE
                GOTO 13
              ENDIF
C             Storing the element:
              RAM(IIKK)=AA
  13        CONTINUE
            NCTMP=NCTMP+NCCOL
  15      CONTINUE
          IF (NSPARC.GE.0) THEN
C           Changing the remaining part of C to sparse matrix:
            CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OA-1,
     *                  OC,NCTMP,ISTORC,M3,OCTMP1)
          ENDIF
        ENDIF
        GOTO 100
      ENDIF
C
C     Reading input matrix B:
      IF (((MATT2.EQ.0).AND.(M2B.NE.M2)).OR.
     *    ((MATT2.EQ.1).AND.(M3 .NE.M2))) THEN
C       MATMUL-07
        CALL ERROR('MATMUL-07: Mismatch in dimensions of the matrices.')
      ENDIF
      IF ((ISYM3.LE.3).AND.(((MATT2.EQ.0).AND.(M1.NE.M3)).OR.
     *                      ((MATT2.EQ.1).AND.(M2.NE.M3)))) THEN
C       MATMUL-08
        CALL ERROR('MATMUL-08: Wrong symmetry of output matrix.')
C       If the output matrix should have symmetry 'diag', 'sym' or
C       'skew', it must have the same number of rows and columns.
      ENDIF
      IF (NSPAR2.GE.0) THEN
        NB=2*NEL2
        OB=1
        OTMP=NB+1
        MTMP=MRAM-NA-NB
      ELSE
        NB=NEL2
        OB=MRAM-NA-NB+1
        OTMP=1
        MTMP=OB-1
      ENDIF
      IF (NB+1.GT.MRAM-NA) THEN
C       MATMUL-09
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATMUL-09: Array RAM too small,',NB+1-MRAM+NA,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
      CALL RMATD(LU1,FILED2,NSPAR2,FORM2,NEL2,IRAM(OB),RAM(OB))
      ISTORB=0
C     Checking whether the non-sparse matrix is to be changed to the
C     sparse matrix, changing the matrix:
      IF (NSPAR2.LT.0) THEN
        RSPAR2=0.5
        CALL GSMAT(M2B,M3,SYM2,NSPAR2,NEL2,RAM,IRAM,OA-1,
     *                  OB,NB,ISTORB,RSPAR2)
        IF (NSPAR2.GE.0) THEN
          OTMP=NB+1
          MTMP=MRAM-NA-NB
        ENDIF
      ELSE
        RSPAR2=FLOAT(NSPAR2)/FLOAT(NSPAR2+NEL2)
      ENDIF
C     Transposing input matrix B:
      IF (MATT2.EQ.1) THEN
        CALL TMATR(M2B,M3,SYM2,NSPAR2,NEL2,IRAM(OB),RAM(OB),
     *             IRAM(OTMP),RAM(OTMP),MTMP)
      ENDIF
C     Extending sparse matrix B:
      IF ((NSPAR2.GE.0).AND.((SYM2.EQ.'sym').OR.(SYM2.EQ.'skew'))) THEN
        CALL SMATRE
     *       (M2,M3,SYM2,NSPAR2,NEL2,IRAM(OB),RAM(OB),OA-1,NB,ISTORB)
        ISYM2=ISYM(SYM2)
      ENDIF
C     Reorganizing sparse matrix B in the memory:
      IF (NSPAR2.GE.0) THEN
        CALL SMATR(M2,M3,NSPAR2,NEL2,IRAM,RAM,OA-1,
     *             OB,NB,ISTORB)
      ENDIF
C
C     Optimizations for the speed of calculation:
      IF (((1.-RSPAR1)*(1.-RSPAR2).LE.0.5).AND.
     *    (.NOT.((NSPAR1.LT.0).AND.(ISYM1.EQ.2.OR.ISYM1.EQ.3))).AND.
     *    (.NOT.((NSPAR2.LT.0).AND.(ISYM2.EQ.2.OR.ISYM2.EQ.3)))) THEN
C       Matrices are to be multiplied as sparse matrices:
C       Checking if there is enough of memory for multiplication:
C       Calculating OC for the case that B will be sparse:
        CALL OMULS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,OA-1,
     *             OB,NB,ISTORB,OC)
        OAN=OA-OC+1
C       Shifting OC if A is to be changed to sparse:
        IF (NSPAR1.LT.0) THEN
          OC=OC + NA - (2*NINT((1.-RSPAR1)*FLOAT(NA))+M2+1)
        ENDIF
        IF (OC.GT.0) THEN
C         Matrices may be (and will be) changed to sparse matrices
C         and shifted to optimum position for multiplication:
C         Shifting B to the beginning of RAM:
          CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM,
     *                OB,NB,ISTORB,1)
          IF (NSPAR2.LT.0) THEN
C           Changing B to sparse, form "as in the memory":
            ISTORB=1
            CALL GSMATR(M2,M3,SYM2,NSPAR2,NEL2,RAM,IRAM,OA-1,
     *                  OB,NB,ISTORB)
          ENDIF
          IF (NSPAR1.LT.0) THEN
C           Shifting A behind B:
            CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM,
     *                  OA,NA,ISTORA,NB+1)
C           Changing A to sparse, form "as in the memory":
            ISTORA=1
            CALL GSMATR(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM,
     *                  OA,NA,ISTORA)
          ENDIF
C         Shifting A to the optimum position for multiplication:
          CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM,
     *                OA,NA,ISTORA,OAN)
C         Shifting B to the optimum position for multiplication:
          CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM,
     *                OB,NB,ISTORB,OAN-NB)
        ENDIF
      ELSE
C       Matrices are to be multiplied as non-sparse matrices:
C       Calculating OC for the case that B will be non-sparse:
        CALL OMULNS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,OA-1,
     *              OB,NB,ISTORB,OC)
        OAN=OA-OC+1
C       Shifting OC if A is to be changed to non-sparse:
        IF (NSPAR1.GE.0) THEN
          OC=OC + NA - NELMAT(M1,M2,ISYM1)
        ENDIF
C       Shifting B to the beginning of RAM:
        CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM,
     *              OB,NB,ISTORB,1)
        IF (NSPAR2.GE.0) THEN
C         Changing B to non-sparse:
          CALL SGMATR(M2,M3,SYM2,NSPAR2,NEL2,RAM,IRAM,OA-1,
     *                OB,NB,ISTORB)
        ENDIF
        IF (NSPAR1.GE.0) THEN
C         Shifting A behind B:
          CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM,
     *                OA,NA,ISTORA,NB+1)
C         Changing A to non-sparse:
          CALL SGMATR(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM,
     *                OA,NA,ISTORA)
        ENDIF
C       Shifting A to the optimum position for multiplication:
        IF (OC.LE.0) OAN=MRAM-NA+1
        CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM,
     *              OA,NA,ISTORA,OAN)
C       Shifting B to the optimum position for multiplication:
        CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM,
     *              OB,NB,ISTORB,OAN-NB)
      ENDIF
C
C     Multiplication:
      WRITE(*,'(A)') '+MATMUL: Multiplying...            '
C
*      IF ((SYM2.EQ.'diag').AND.
*     *    (((SYM1.EQ.'DIAG').OR.(SYM1.EQ.' ')).AND.(SYM1.EQ.SYM3)))THEN
*C       Matrix B is diagonal, A is diagonal or general,
*C       thus C will be A multiplied by the diagonal of B:
*        OC=OA
*        NC=NA
*        NSPARC=NSPAR1
*        NELC=NEL1
*        SYMC=SYM1
*        ISYMC=ISYM1
*        ISTORC=ISTORA
*C       Loop over diagonal of B:
*        DO 10, KK=1,M2
*          WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M2
*C         Calculating the value BB of diagonal element of B:
*          IF (NSPAR2.GE.0) THEN
*            IF (IRAM(OB+KK).EQ.IRAM(OB+KK-1)) THEN
*              BB=0.
*            ELSE
*              BB=RAM(IRAM(OB+KK))
*            ENDIF
*          ELSE
*            BB=RAM(OB-1+KK)
*          ENDIF
*C         Preparing loop over column KK of A:
*          IF (NSPAR1.GE.0) THEN
*            II1=IRAM(OA+KK-1)+2
*            II2=IRAM(OA+KK)
*            II3=2
*          ELSE
*            IF (SYM1.EQ.'diag') THEN
*              II1=OA+KK-1
*              II2=II1
*              II3=1
*            ELSEIF (SYM1.EQ.'sym') THEN
*              II1=OA+(KK-1)*KK/2
*              II2=II1+KK-1
*              II3=1
*            ELSEIF (SYM1.EQ.'skew') THEN
*              II1=OA+(KK-1)*(KK-2)/2
*              II2=II1+KK-2
*              II3=1
*            ELSE
*              II1=OA+(KK-1)*M1
*              II2=II1+M1-1
*              II3=1
*            ENDIF
*          ENDIF
*C         Loop over column KK of A:
*          DO 8, II=II1,II2,II3
*C           Multiplying:
*            RAM(II)=RAM(II)*BB
*   8      CONTINUE
*  10    CONTINUE
*C       Shifting matrix C:
*        IF (NSPARC.GE.0) THEN
*          DO 12, I1=1,M3+1
*            IRAM(I1)=IRAM(OC-1+I1)-OC+1
*  12      CONTINUE
*          DO 13, I1=1,NC-M3-1-1,2
*            IRAM(1+M3+I1)=IRAM(OC+M3+I1)
*            RAM(1+M3+I1+1)=RAM(OC+M3+I1+1)
*  13      CONTINUE
*        ELSE
*          DO 14, I1=1,NC
*            RAM(I1)=RAM(OC-1+I1)
*  14      CONTINUE
*        ENDIF
*        OC=1
*      ELSEIF ((SYM1.EQ.'diag').AND.
*     *    (((SYM2.EQ.'DIAG').OR.(SYM2.EQ.' ')).AND.(SYM2.EQ.SYM3)))THEN
*C       Matrix A is diagonal, B is diagonal or general,
*C       thus C will be B multiplied by the diagonal of A:
*        OC=OB
*        NC=NB
*        NELC=NEL2
*        NSPARC=NSPAR2
*        SYMC=SYM2
*        ISYMC=ISYM2
*        ISTORC=ISTORB
*C       Loop over all columns of B:
*        DO 20, KK=1,M3
*          WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M3
*C         Preparing loop over column KK of B:
*          IF (NSPAR2.GE.0) THEN
*            JJ1=IRAM(OB+KK-1)+2
*            JJ2=IRAM(OB+KK)
*            JJ3=2
*          ELSE
*            IF (SYM2.EQ.'sym') THEN
*              JJ1=OB+(KK-1)*KK/2
*              JJ2=JJ1+KK-1
*              JJ3=1
*            ELSEIF (SYM2.EQ.'skew') THEN
*              JJ1=OB+(KK-1)*(KK-2)/2
*              JJ2=JJ1+KK-2
*              JJ3=1
*            ELSE
*              JJ1=OB+(KK-1)*M2
*              JJ2=JJ1+M2-1
*              JJ3=1
*            ENDIF
*          ENDIF
*C         Loop over column KK of B:
*          DO 18, JJ=JJ1,JJ2,JJ3
*C           Calculating index II of row of B:
*            IF (NSPAR2.GE.0) THEN
*              II=IRAM(JJ-1)
*            ELSE
*              II=(JJ-JJ1+JJ3)/JJ3
*            ENDIF
*C           Calculating the value AA of diagonal element of A:
*            IF (NSPAR1.GE.0) THEN
*              IF (IRAM(OA+II).EQ.IRAM(OA+II-1)) THEN
*                AA=0.
*              ELSE
*                AA=RAM(IRAM(OA+II))
*              ENDIF
*            ELSE
*              AA=RAM(OA-1+II)
*            ENDIF
*C           Multiplying:
*            RAM(JJ)=RAM(JJ)*AA
*  18      CONTINUE
*  20    CONTINUE
*C       Shifting matrix C:
*        IF (NSPARC.GE.0) THEN
*          DO 22, I1=1,M3+1
*            IRAM(I1)=IRAM(OC+I1-1)-OC+1
*  22      CONTINUE
*          DO 23, I1=1,NC-M3-1-1,2
*            IRAM(1+M3+I1)=IRAM(OC+M3+I1)
*            RAM(1+M3+I1+1)=RAM(OC+M3+I1+1)
*  23      CONTINUE
*        ELSE
*          DO 24, I1=1,NC
*            RAM(I1)=RAM(OC+I1-1)
*  24      CONTINUE
*        ENDIF
*        OC=1
*      ELSEIF ((SYM1.EQ.' ').AND.(SYM2.EQ.' ').AND.
      IF ((NSPAR1.LT.0).AND.(NSPAR2.LT.0).AND.
ccc  *    ((M1.LT.OB).AND.(M1*M3.LT.OA-M2))) THEN
     *    (OC.GT.0)) THEN
C       Multiplication for non-sparse matrices:
        OC=1
        NC=M1*M3
        NSPARC=-1
        NELC=NC
        SYMC=SYM3
        ISYMC=ISYM3
        ISTORC=0
        IIKK=0
        II1=1
        II2=M1
        OA1=OA-1
        OB1=OB-1
C       Loop over columns of C:
        DO 27, KK=1,M3
          WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M3
C         Loop over lines of column KK of matrix C:
          IF     (ISYMC.EQ.1) THEN
            II1=KK
            II2=KK
          ELSEIF (ISYMC.EQ.2) THEN
            II2=KK
          ELSEIF (ISYMC.EQ.3) THEN
            II2=KK-1
          ENDIF
          DO 26, II=II1,II2
            IIKK=IIKK+1
            CIK=0.
            DO 25, JJ=1,M2
C             Element of the matrix A:
              IF     (ISYM1.EQ.1) THEN
                IF (II.EQ.JJ) THEN
                  AA=RAM(OA1+JJ)
                ELSE
                  GOTO 25
                ENDIF
              ELSEIF (ISYM1.EQ.2) THEN
                IF (II.LE.JJ) THEN
                  AA=RAM(OA1+(JJ-1)*JJ/2+II)
                ELSE
                  AA=RAM(OA1+(II-1)*II/2+JJ)
                ENDIF
              ELSEIF (ISYM1.EQ.3) THEN
                IF (II.EQ.JJ) THEN
                  GOTO 25
                ELSEIF (II.LT.JJ) THEN
                  AA= RAM(OA1+(JJ-1)*(JJ-2)/2+II)
                ELSE
                  AA=-RAM(OA1+(II-1)*(II-2)/2+JJ)
                ENDIF
              ELSE
                AA=RAM(OA1+(JJ-1)*M1+II)
              ENDIF
C             Element of the matrix B:
              IF     (ISYM2.EQ.1) THEN
                IF (JJ.EQ.KK) THEN
                  BB=RAM(OB1+KK)
                ELSE
                  GOTO 25
                ENDIF
              ELSEIF (ISYM2.EQ.2) THEN
                IF (JJ.LE.KK) THEN
                  BB=RAM(OB1+(KK-1)*KK/2+JJ)
                ELSE
                  BB=RAM(OB1+(JJ-1)*JJ/2+KK)
                ENDIF
              ELSEIF (ISYM2.EQ.3) THEN
                IF (JJ.EQ.KK) THEN
                  GOTO 25
                ELSEIF (JJ.LT.KK) THEN
                  BB= RAM(OB1+(KK-1)*(KK-2)/2+JJ)
                ELSE
                  BB=-RAM(OB1+(JJ-1)*(JJ-2)/2+KK)
                ENDIF
              ELSE
                BB=RAM(OB1+(KK-1)*M2+JJ)
              ENDIF
              CIK=CIK+AA*BB
  25        CONTINUE
            RAM(IIKK)=CIK
  26      CONTINUE
  27    CONTINUE
      ELSE
C       Multiplication for the case that at least one matrix is sparse,
C       or OC is negative (no space for whole matrix C):
        OC=1
        IF (SYM3.EQ.'sym') THEN
C         Symmetric matrix
          NC=M1*(M1+1)/2
        ELSEIF (SYM3.EQ.'skew') THEN
C         Antisymmetric matrix
          NC=M1*(M1-1)/2
        ELSEIF (SYM3.EQ.'diag') THEN
C         Diagonal matrix
          NC=M1
        ELSE
C         General matrix
          NC=M1*M3
        ENDIF
        MC=OB-1
        NSPARC=-1
        NELC=NC
        SYMC=SYM3
        ISYMC=ISYM3
        ISTORC=0
C
        NCTMP=0
        OCTMP1=0
        DO 50, KK=1,M3
          WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M3
C         Preparing sufficient number of the storage locations:
C         Number NCCOL of storage locations for column KK:
          IF     (ISYMC.EQ.1) THEN
            NCCOL=1
          ELSEIF (ISYMC.EQ.2) THEN
            NCCOL=KK
          ELSEIF (ISYMC.EQ.3) THEN
            NCCOL=KK-1
          ELSE
            NCCOL=M1
          ENDIF
          IF (NCTMP+NCCOL.GT.MC) THEN
C           Erasing unneeded columns 1 to KK-1 of B:
            IF (NSPAR2.LT.0) THEN
C             OB stays unchanged, MC is redefined:
              IF     (ISYM2.EQ.1) THEN
                MC=OB+KK-1-1
              ELSEIF (ISYM2.EQ.2) THEN
                MC=OB+(KK-1)*KK/2-1
              ELSEIF (ISYM2.EQ.3) THEN
                MC=OB+(KK-2)*(KK-1)/2-1
              ELSE
                MC=OB+M2*(KK-1)-1
              ENDIF
            ELSE
C             OB must be changed, addresses of the ends of the columns
C             must be moved:
              OBN=IRAM(OB+KK-1)-M3
              DO 28, I1=M3,0,-1
                IF (I1.GE.KK) THEN
                  IRAM(OBN+I1)=IRAM(OB+I1)
                ELSE
                  IRAM(OBN+I1)=IRAM(OB+KK-1)
                ENDIF
  28          CONTINUE
              OB=OBN
              MC=OB-1
            ENDIF
            IF (NCTMP+NCCOL.GT.MC) THEN
C             Changing the calculated part of C to sparse matrix:
              CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OB-1,
     *                    OC,NCTMP,ISTORC,KK-1,OCTMP1)
            ENDIF
            IF (NCTMP+NCCOL.GT.MC) THEN
C             MATMUL-10
              WRITE(TXTERR,'(A,I9,A)')
     *       'MATMUL-10: Array RAM too small,',NCTMP+NCCOL-MC,
     *       ' units missing.'
              CALL ERROR(TXTERR)
            ENDIF
          ENDIF
C         Initiating the values in column KK of C:
          DO 29, I1=NCTMP+1,NCTMP+NCCOL
            RAM(I1)=0.
  29      CONTINUE
C         Preparing loop over column KK of B:
          CALL MMILOC(M2,M3,ISYM2,NSPAR2,RAM,IRAM,OB,KK,JJ1,JJ2,JJ3)
C         Loop over nonzero elements of column KK of B:
          DO 40, JJA=JJ1,JJ2,JJ3
C           Computing index JJ of row and value BB of the element of B:
            CALL MMIRVE(M2,M3,ISYM2,NSPAR2,RAM,IRAM,OB,KK,JJA,JJ1,JJ3,
     *                                     JJ,BB)
C           Preparing loop over column JJ of A:
            CALL MMILOC(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,JJ,II1,II2,II3)
C           Loop over nonzero elements of column JJ of A:
            DO 30, IIA=II1,II2,II3
C             Computing index II of row and value AA of element of A:
              CALL MMIRVE(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,JJ,IIA,II1,II3,
     *                                     II,AA)
              IF     ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN
                IIKK=OCTMP1+II
              ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN
                IIKK=OCTMP1+(KK-1)*KK/2+II
              ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN
                IIKK=OCTMP1+(KK-2)*(KK-1)/2+II
              ELSEIF (ISYMC.EQ.4)                  THEN
                IIKK=OCTMP1+(KK-1)*M1+II
              ELSE
                GOTO 30
              ENDIF
C             Multiplying:
              RAM(IIKK)=RAM(IIKK) + AA*BB
  30        CONTINUE
  40      CONTINUE
          NCTMP=NCTMP+NCCOL
  50    CONTINUE
        IF (NSPARC.GE.0) THEN
C         Changing the remaining part of C to sparse matrix:
          CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OB-1,
     *                OC,NCTMP,ISTORC,M3,OCTMP1)
        ENDIF
      ENDIF
C
 100  CONTINUE
C
      IF ((ISPAR3.EQ.1).AND.(NSPARC.LT.0)) THEN
C       Changing non-sparse matrix to sparse matrix:
        CALL GSMATR(M1,M3,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,OC,NC,ISTORC)
      ELSEIF ((ISPAR3.EQ.-1).AND.(NSPARC.GE.0)) THEN
C       Changing sparse matrix to non-sparse matrix:
        CALL SGMATR(M1,M3,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,OC,NC,ISTORC)
      ELSEIF (ISPAR3.EQ.0) THEN
C       Automatic selection of the sparseness:
        IF (NSPARC.LT.0) THEN
C         Matrix is non-sparse, it will be changed to sparse if its
C         sparseness is 0.5 or more:
          RSPAR3=0.5
          CALL GSMAT(M1,M3,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,
     *                  OC,NC,ISTORC,RSPAR3)
        ELSE
C         Matrix is sparse.
          RSPAR3=FLOAT(NSPARC)/FLOAT(NSPARC+NELC)
          IF (RSPAR3.LT.0.5) THEN
C           Matrix is sparse, sparseness is less than 0.5, thus changing
C           the sparse matrix to non-sparse matrix:
            CALL SGMATR(M1,M3,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,
     *                  OC,NC,ISTORC)
          ENDIF
        ENDIF
      ENDIF
C
      IF ((NSPARC.GE.0).AND.(ISTORC.EQ.1)) THEN
C       Reorganizing sparse matrix in the memory:
        CALL SMATR(M1,M3,NSPARC,NELC,IRAM,RAM,MRAM,
     *             OC,NC,ISTORC)
      ENDIF
C
C
C     Writing output matrix C:
      CALL WMATH(LU1,FILE3,FILED3,M1,M3,NSPARC,SYM3,FORM3)
      CALL WMATD(LU1,FILED3,NSPARC,FORM3,NELC,IRAM(OC),RAM(OC))
C
      WRITE(*,'(A)') '+MATMUL: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE MMILOC(M1,M2,ISYM,NSPARS,ARRAY,IARRAY,OARRAY,ICOLUM,
     *                 IELEM1,IELEM2,IELEM3)
      INTEGER M1,M2,ISYM,NSPARS,OARRAY,ICOLUM,IELEM1,IELEM2,IELEM3
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine designed to Initiate Loop Over Column ICOLUM of a matrix.
C
C Input:
C     M1 ...  Number of rows of the matrix.
C     M2 ...  Number of columns of the matrix.
C     ISYM... Index of the symmetry of the matrix.
C     NSPARS . Sparseness of the matrix.
C     ARRAY,IARRAY ... Array with indices and values of the matrix.
C             Sparse matrix must be stored in the form "as in memory".
C     OARRAY..Address of the first storage location of the matrix.
C     ICOLUM..Number of the column under considertation.
C     For detailed description of storage of matrices in the memory
C     refer to file mat.for.
C
C Output:
C     IELEM1..Address of the value of the first element of column ICOLUM
C             in array ARRAY.
C     IELEM2..Address of the value of the last element of column ICOLUM
C             in array ARRAY.
C     IELEM3..Step between two consecutive values of the matrix elements
C             in array ARRAY.
C
C Coded by Petr Bulant
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C.......................................................................
C
C     Preparing loop over column ICOLUM of the matrix:
      IF (NSPARS.GE.0) THEN
        IELEM1=IARRAY(OARRAY+ICOLUM-1)+2
        IELEM2=IARRAY(OARRAY+ICOLUM)
        IELEM3=2
      ELSE
        IF     (ISYM.EQ.1) THEN
          IELEM1=OARRAY+ICOLUM-1
          IELEM2=IELEM1
          IELEM3=1
        ELSEIF (ISYM.EQ.2) THEN
          IELEM1=OARRAY+(ICOLUM-1)*ICOLUM/2
          IELEM2=IELEM1+M1-1
          IELEM3=1
        ELSEIF (ISYM.EQ.3) THEN
          IELEM1=OARRAY+(ICOLUM-1)*(ICOLUM-2)/2
          IELEM2=IELEM1+M1-1
          IELEM3=1
        ELSE
          IELEM1=OARRAY+(ICOLUM-1)*M1
          IELEM2=IELEM1+M1-1
          IELEM3=1
        ENDIF
      ENDIF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE MMIRVE(M1,M2,ISYM,NSPARS,ARRAY,IARRAY,OARRAY,ICOLUM,
     *                 IELEMA,IELEM1,IELEM3,
     *                                     IROW,VELEM)
      INTEGER M1,M2,ISYM,NSPARS,OARRAY,ICOLUM,IELEMA,IELEM1,IELEM3,IROW
      REAL     ARRAY(*),VELEM
      INTEGER IARRAY(*)
C
C Subroutine designed to calculate Index of Row and Value of Element
C of a matrix.
C
C Input:
C     M1 ...  Number of rows of the matrix.
C     M2 ...  Number of columns of the matrix.
C     ISYM... Index of the symmetry of the matrix.
C     NSPARS. Sparseness of the matrix.
C     ARRAY,IARRAY ... Array with indices and values of the matrix.
C             Sparse matrix must be stored in the form "as in memory".
C     OARRAY..Address of the first storage location of the matrix.
C     ICOLUM..Number of the column under considertation.
C     IELEMA..Address of the value of the current element of column
C             ICOLUM in array ARRAY.
C     IELEM1..Address of the value of the first element of column ICOLUM
C             in array ARRAY.
C     IELEM3..Step between two consecutive values of the matrix elements
C             in array ARRAY.
C     For detailed description of storage of matrices in the memory
C     refer to file mat.for.
C
C Output:
C     IROW .. Number of the row corresponding to the matrix element
C             under considertation.
C     VELEM...Value of the considered matrix element.
C
C Coded by Petr Bulant
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C.......................................................................
C
C     Calculating index IROW of row of the matrix:
      IF (NSPARS.GE.0) THEN
        IROW=IARRAY(IELEMA-1)
      ELSE
        IF (ISYM.EQ.1) THEN
          IROW=ICOLUM
        ELSE
          IROW=(IELEMA-IELEM1+IELEM3)/IELEM3
        ENDIF
      ENDIF
C     Calculating the value VELEM of element of the matrix:
      IF (NSPARS.GE.0) THEN
        VELEM=ARRAY(IELEMA)
      ELSE
        IF     (ISYM.EQ.2) THEN
C         'sym'
          IF (IROW.LE.ICOLUM) THEN
            VELEM=ARRAY(IELEMA)
          ELSE
            VELEM=ARRAY(OARRAY-1+(IROW-1)*IROW/2+ICOLUM)
          ENDIF
        ELSEIF (ISYM.EQ.3) THEN
C         'skew'
          IF (IROW.LT.ICOLUM) THEN
            VELEM=ARRAY(IELEMA)
          ELSEIF (IROW.EQ.ICOLUM) THEN
            VELEM=0.
          ELSE
            VELEM=-ARRAY(OARRAY-1+(IROW-1)*(IROW-2)/2+ICOLUM)
          ENDIF
        ELSE
C         'diag' or ' '
          VELEM=ARRAY(IELEMA)
        ENDIF
      ENDIF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OMULS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,
     *                 IARRAY,ARRAY,MARRAY,OARRAY,NARRAY,ISTOR,OC)
      INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MARRAY,OARRAY,NARRAY,
     *        ISTOR,OC
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine to estimate optimum origin OC where to write matrix C,
C assuming that matrix B given to OMUL will be stored as sparse matrix,
C and will be stored so that it will terminate at the same position
C where it terminates now.
C OC is estimated as minimum over KK=1,M2 of End(KK-1)-M2-AllKKC,
C where End(KK-1) is the address of end of row KK-1 of matrix B if
C matrix B were stored as sparse matrix in the form "as in the memory",
C and AllKKC is the number of elements of columns 1 to KK of matrix C.
C
C Input:
C     M1C ... Number of rows of matrix C.
C     ISYMC.. Index of symmetry of matrix C.
C     M1...   Number of rows of matrix B.
C     M2...   Number of columns of matrix B.
C     ISYM .. Index of symmetry of matrix B.
C     NSPAR...Sparseness of matrix B.
C     NELEM...Number of elements of matrix B.
C     IARRAY..Components of matrix B.
C     ARRAY...Components of matrix B.
C     MARRAY..Dimension of array ARRAY.
C     OARRAY..First storage location of matrix B.
C     NARRAY..Number of storage locations for matrix B.
C     ISTOR...Form of matrix B.
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C
C Output:
C     OC  ... Origin where to write matrix C.
C
C Coded by Petr Bulant
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NELMAT
      INTEGER NELMAT
C
C     Local storage location:
      INTEGER II,KK,JJ,NN,IB,II1,II2,OTMP
C
C.......................................................................
C
      IF (NSPAR.GE.0.AND.ISTOR.EQ.0) THEN
C       MATMUL-11
        CALL ERROR('MATMUL-11: Wrong ISTOR.')
C       This option is not yet coded, sparse matrix must be in the form
C       "as in the memory".
      ENDIF
C
      IF (NSPAR.LT.0) THEN
C       Calculating addresses of ends of columns of B if B were stored
C       as sparse matrix ending at the same position where it ends now:
        OTMP=OARRAY-M2-1
        IF (OTMP.LE.0) THEN
C         There is no space for ends of rows, thus OC is put to -1:
          OC=-1
          RETURN
        ENDIF
        IB=OARRAY+NARRAY-1
        JJ=IB
        IARRAY(OTMP+M2)=JJ
        II1=1
        II2=M1
C       Loop over columns of B:
        DO 12, KK=M2,1,-1
C         Loop over lines of B:
          IF     (ISYM.EQ.1) THEN
            II1=KK
            II2=KK
          ELSEIF (ISYM.EQ.2) THEN
            II2=KK
          ELSEIF (ISYM.EQ.3) THEN
            II2=KK-1
          ENDIF
C         Number NN of nonzero elements of column KK:
          NN=0
          DO 10, II=II1,II2
            IF (ARRAY(IB).NE.0.) NN=NN+1
            IB=IB-1
  10      CONTINUE
          JJ=JJ-2*NN
          IARRAY(OTMP+KK-1)=JJ
  12    CONTINUE
      ELSE
        OTMP=OARRAY
      ENDIF
      OC=OTMP
      DO 20, KK=1,M2
        OC=MIN0(OC,IARRAY(OTMP+KK-1)-M2-NELMAT(M1C,KK,ISYMC))
  20  CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OMULNS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,
     *                  IARRAY,ARRAY,MARRAY,OARRAY,NARRAY,ISTOR,OC)
      INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MARRAY,OARRAY,NARRAY,
     *        ISTOR,OC
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine to estimate optimum origin OC where to write matrix C,
C assuming that matrix B given to OMUL will be stored as non-sparse
C matrix, and will be stored so that it will terminate at the same
C position where it terminates now.
C OC is estimated as minimum over KK=1,M2 of End(KK-1)-AllKKC+1,
C where End(KK-1) is the address of end of row KK-1 of matrix B if
C matrix B were stored as non-sparse matrix,
C and AllKKC is the number of elements of columns 1 to KK of matrix C.
C
C Input:
C     M1C ... Number of rows of matrix C.
C     ISYMC.. Index of symmetry of matrix C.
C     M1...   Number of rows of matrix B.
C     M2...   Number of columns of matrix B.
C     ISYM .. Index of symmetry of matrix B.
C     NSPAR...Sparseness of matrix B.
C     NELEM...Number of elements of matrix B.
C     IARRAY..Components of matrix B.
C     ARRAY...Components of matrix B.
C     MARRAY..Dimension of array ARRAY.
C     OARRAY..First storage location of matrix B.
C     NARRAY..Number of storage locations for matrix B.
C     ISTOR...Form of matrix B.
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C
C Output:
C     OC  ... Origin where to write matrix C.
C
C Coded by Petr Bulant
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL NELMAT
      INTEGER NELMAT
C
C     Local storage location:
      INTEGER KK,OTMP
C
C.......................................................................
C
      IF (NSPAR.LT.0) THEN
        OTMP=OARRAY
      ELSE
        OTMP=OARRAY+NARRAY-NELMAT(M1,M2,ISYM)
      ENDIF
      OC=OTMP
      DO 20, KK=1,M2
        OC=MIN0(OC,OTMP+NELMAT(M1,KK-1,ISYM)-NELMAT(M1C,KK,ISYMC))
  20  CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE GSPART(M1,M2,ISYM,NSPAR,NELEM,ARRAY,IARRAY,MARRAY,
     *                  OARRAY,NARRAY,ISTOR,ICOL,OTMP)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,OARRAY,NARRAY,MARRAY,ISTOR,
     *        ICOL,OTMP
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine designed to consecutively change parts of non-sparse matrix
C to the sparse matrix. For the first invocation, it is assumed on input
C that non-sparse part of the matrix for columns 1 to ICOL is stored
C in ARRAY, and this part is changed to sparse matrix, form "as in the
C memory" and written to ARRAY and IARRAY together with M2+1 addresses
C of ends of columns. For second and other invocations, it is assumed
C that IARRAY and ARRAY contain M2+1 addresses of ends of columns of the
C sparse part of the matrix, followed by the sparse part of the matrix
C (columns 1 to ICOL0), and finally by non-sparse part of the matrix
C (columns ICOL0+1 to ICOL). The non-sparse part is converted to the
C sparse and written after the already stored sparse part, and the
C addressses of columns ICOL0+1 to M2 are updated.
C
C Input:
C     M1...   Number of rows of the full matrix.
C     M2...   Number of columns of the full matrix.
C     ISYM ...Index of symmetry of the matrix.
C     ARRAY...Components of the given matrix.
C     IARRAY..Components of the given matrix.
C     MARRAY..Dimension of array ARRAY.
C     OARRAY..Address of the first storage location in array ARRAY
C             used for the matrix.
C     NARRAY..Number of storage locations for the matrix.
C     ICOL ...Index of the last calculated column of the matrix.
C
C Output:
C     NSPAR...Sparseness of the M1*ICOL matrix.
C     NELEM...Number of elements of the M1*ICOL matrix.
C     NARRAY..New number of storage locations for the M1*ICOL matrix.
C     IARRAY..Components of the M1*ICOL matrix.
C     ARRAY...Components of the M1*ICOL matrix.
C     ISTOR...Form of the matrix.
C             ISTOR=1 ... form "as in the memory"
C     OTMP ...Fictitious position in ARRAY just before the first storage
C             location of the matrix if the matrix would stay non-sparse
C             and would terminate at the same location where it
C             terminates after its conversion to the sparse matrix.
C
C Coded by Petr Bulant
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NSPMAT,NELMAT
      INTEGER NSPMAT,NELMAT
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,I20,I11,I12,J1,OAR,NAR,NSPA,NELE,ICOL0
      CHARACTER*72 TXTERR
      SAVE ICOL0
      DATA ICOL0/0/
C
C.......................................................................
C
C     Beginning OAR and number of elements NAR of the non-sparse part
C     of ARRAY:
      IF (ICOL0.EQ.0) THEN
        OAR=OARRAY
        NAR=NARRAY
      ELSE
        OAR=IARRAY(OARRAY+M2)+1
        NAR=NARRAY-(OAR-OARRAY)
      ENDIF
C     Number NSPA of zero and NELE of nonzero elements:
      NSPA=NSPMAT(NAR,ARRAY(OAR))
      NELE=NAR-NSPA
      IF (((ICOL0.EQ.0).AND.(OAR+M2+1+2*NELE.GT.MARRAY)).OR.
     *    ((ICOL0.NE.0).AND.(OAR+     2*NELE.GT.MARRAY))) THEN
C       MATMUL-12
        IF (ICOL0.EQ.0) THEN
          I1=OAR+M2+1+2*NELE - MARRAY
        ELSE
          I1=OAR+     2*NELE - MARRAY
        ENDIF
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATMUL-12: Array RAM too small,',I1,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     Moving the nonzero elements of the matrix to the end of ARRAY,
C     storing the corresponding matrix indices of the elements:
      J1=MARRAY
      I3=OARRAY-1+NARRAY
      DO 20, I2=ICOL,ICOL0+1,-1
        IF     (ISYM.EQ.1) THEN
          I11=I2
          I12=I2
        ELSEIF (ISYM.EQ.2) THEN
          I11=I2
          I12=1
        ELSEIF (ISYM.EQ.3) THEN
          I11=I2-1
          I12=1
        ELSE
          I11=M1
          I12=1
        ENDIF
        DO 18, I1=I11,I12,-1
          IF (ARRAY(I3).NE.0.) THEN
            IF (J1-1.LT.I3) THEN
C             MATMUL-13
              WRITE(TXTERR,'(A,I9,A)')
     *        'MATMUL-13: Array RAM too small,',I3-J1+1,
     *        ' units missing.'
              CALL ERROR(TXTERR)
            ENDIF
            ARRAY(J1)=ARRAY(I3)
            IARRAY(J1-1)=(I2-1)*M1+I1
            J1=J1-2
          ENDIF
          I3=I3-1
  18    CONTINUE
  20  CONTINUE
      IF (I3.NE.OAR-1) THEN
C       MATMUL-14
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATMUL-14: Array RAM too small,',OAR-1-I3,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     J1 now points just before the index of the first nonzero element
C
      ISTOR=1
      IF (ICOL0.EQ.0) OAR=OARRAY+M2+1
C     Rearranging the matrix to the form 'as in the memory':
C     Moving the indices and the values to the new positions:
      DO 42, I3=1,NELE
        ARRAY(OAR-1+2*I3)=ARRAY(J1+2*I3)
        IARRAY(OAR-1+2*I3-1)=IARRAY(J1+2*I3-1)
  42  CONTINUE
C     Calculating end addresses of the columns,
C     changing matrix indices of the elements to the row numbers:
      I20=M2+1
      DO 44, I3=NELE,1,-1
        I1=MOD(IARRAY(OAR-1+2*I3-1),M1)
        IF (I1.EQ.0) I1=M1
        I2=(IARRAY(OAR-1+2*I3-1)-1)/M1+1
        IARRAY(OAR-1+2*I3-1)=I1
        IF (I2.NE.I20) THEN
          DO 43, I4=I20-1,I2,-1
            IARRAY(OARRAY+I4)=OAR-1+2*I3
  43      CONTINUE
          I20=I2
        ENDIF
  44  CONTINUE
      DO 45, I4=I20-1,ICOL0,-1
        IARRAY(OARRAY+I4)=OAR-1
  45  CONTINUE
C
      ICOL0=ICOL
      NARRAY=IARRAY(OARRAY+M2)-OARRAY+1
      NELEM=(NARRAY-M2-1)/2
      I1=NELMAT(M1,ICOL,ISYM)
      NSPAR=I1-NELEM
      OTMP=IARRAY(OARRAY+M2)-I1
C
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'mat.for'
C     mat.for
      INCLUDE 'indexi.for'
C     indexi.for
C
C=======================================================================
C