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.70
C Date: 2012, November 7
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             For 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             Default: MATIN2=' '
C     MATOUT='string'... Name of the header file of the output matrix.
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     FORMM='string'... Form of the output file with the matrix.
C             Possible values of FORMM are:
C               FORMM='FORMATTED'... formatted file
C               FORMM='UNFORMATTED'... unformatted file
C             Default: FORMM='FORMATTED'
C     SPARSE=integer... Identifies whether the matrix should be sparse.
C             Possible values of SPARSE are:
C               SPARSE=1... sparse matrix in the CSC format
C               SPARSE=0... automatic selection of the sparseness:
C                 matrix with half or more zero elements will be sparse
C                 in the CSC format, otherwise non-sparse matrix
C               SPARSE=-1... normal (not sparse) matrix
C             Default: SPARSE=0 (automatic selection of the sparseness)
C     For general 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     NUMLINM=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-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3I,RSEP3T,RMATH,RMATD,WMATH,WMATD,
     *TMATR,SMATRE,GSMAT,GSMATR,SGMATR,ISYM,NELMAT,MSHIFT,VELEM,
     *MMILOC,MMIRVE,OMULS,OMULNS,LOWER,UPPER
      INTEGER ISYM,NELMAT
      REAL VELEM
C     ERROR... File error.for.
C     RSEP1,RSEP3I,RSEP3T,SSEP...File sep.for.
C     RMATH,RMATD,WMATH,WMATD,TMATR,SMATRE,GSMAT,GSMATR,SGMATR,
C     ISYM,NELMAT,MSHIFT,VELEM... File mat.for.
C     MMILOC,MMIRVE,OMULS,OMULNS,GSPART... This file.
C     LOWER,UPPER... File length.for.
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... File forms.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,FORM1,FORM2,FORM3
      CHARACTER*80 SPARS1,SPARS2,SPARSC
      INTEGER MATT1,MATT2,M1,M2,M2B,M3,NEL1,NEL2,NELC,IA,IA1,IAN,NA,
     *IB,IB1,IBN,NB,IC,NC,MC,ICNEW,NCTMP,NCCOL,NTMP1,NTMP2,
     *NSPAR1,NSPAR2,ISPAR3,NSPARC,ISYM1,ISYM2,ISYM3,ISYMC,LU1,
     *I1,II,JJ,KK,II1,II2,II3,IIA,JJ1,JJ2,JJ3,JJA,IIKK,ICTMP1,JJKK,IIJJ
      REAL AA,BB,CIK,RSPAR1,RSPAR2,RSPAR3,BBB
      LOGICAL LSPARC
      CHARACTER*72 TXTERR
      PARAMETER (LU1=1)
C
C-----------------------------------------------------------------------
C
      MINRAM=1
      MAXRAM=MRAM
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,SPARS1,NEL1,SYM1,FORM1)
      ISYM1=ISYM(SYM1)
      IF (SPARS1.EQ.' ') THEN
        NSPAR1=-1
      ELSEIF (SPARS1.EQ.'CSC') THEN
        NSPAR1=NELMAT(M1,M2,ISYM1)-NEL1
      ELSE
C       MATMUL-15
        CALL ERROR('MATMUL-15: Wrong value of SPARSE of input matrix')
C       Only SPARSE=' ' or SPARSE='CSC' is allowed.
      ENDIF
      IF (FILE2.NE.' ') THEN
        CALL RMATH(LU1,FILE2,FILED2,M2B,M3,SPARS2,NEL2,SYM2,FORM2)
        ISYM2=ISYM(SYM2)
        IF (SPARS2.EQ.' ') THEN
          NSPAR2=-1
        ELSEIF (SPARS2.EQ.'CSC') THEN
          NSPAR2=NELMAT(M2B,M3,ISYM2)-NEL2
        ELSE
C         MATMUL-16
          CALL ERROR('MATMUL-16: Wrong value of SPARSE of input matrix')
C         Only SPARSE=' ' or SPARSE='CSC' is allowed.
        ENDIF
      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('FORMM',FORM3,'FORMATTED')
      CALL UPPER(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=M2+1+2*NEL1
      ELSE
        NA=NEL1
      ENDIF
      IA=MAXRAM-NA+1
      IF (NA+1.GT.MAXRAM) THEN
C       MATMUL-04
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATMUL-04: Array RAM too small,',NA+1-MAXRAM,' units missing.'
        CALL ERROR(TXTERR)
C       This is the memory needed to read the first matrix M1.
      ENDIF
      CALL RMATD(LU1,FILED1,M2,SPARS1,NEL1,FORM1,IA)
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,ISYM1,NSPAR1,NEL1,MINRAM,MAXRAM,IA,NA,RSPAR1)
      ELSE
        RSPAR1=FLOAT(NSPAR1)/FLOAT(NSPAR1+NEL1)
      ENDIF
C     Transposing input matrix A:
      IF (MATT1.EQ.1) THEN
        CALL TMATR(M1,M2,ISYM1,NSPAR1,NEL1,MINRAM,MAXRAM,IA,NA)
        CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,MINRAM,MAXRAM,IA,NA,
     *              MAXRAM-NA+1)
      ENDIF
C     Extending sparse matrix A:
CCC   IF ((NSPAR1.GE.0).AND.((SYM1.EQ.'sym').OR.(SYM1.EQ.'skew'))) THEN
CCC     CALL SMATRE
CCC  *       (M1,M2,SYM1,NSPAR1,NEL1,IRAM(IA),RAM(IA),MAXRAM,NA,ISTORA)
CCC     ISYM1=ISYM(SYM1)
CCC   ENDIF
C
      IF (FILE2.EQ.' ') THEN
C       Matrix A is to be written as output matrix C:
        M3=M2
        IF (ISYM3.EQ.ISYM1) THEN
C         No change of matrix symmetry:
          IC=IA
          NC=NA
          NSPARC=NSPAR1
          NELC=NEL1
          ISYMC=ISYM1
C         Copying matrix A to the position of C:
          ICNEW=MINRAM
          CALL MSHIFT(M1,M3,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC,ICNEW)
        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
          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
          IC=MINRAM
          MC=IA-1
          NSPARC=-1
          NELC=NC
          ISYMC=ISYM3
C
          NCTMP=0
          ICTMP1=IC-1
          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 (IC-1+NCTMP+NCCOL.GT.MC) THEN
C             Erasing unneeded columns 1 to KK-1 of A:
              IF (NSPAR1.LT.0) THEN
C               IA stays unchanged, MC is redefined:
                IF     (ISYM1.EQ.1) THEN
                  MC=IA+KK-1-1
                ELSEIF (ISYM1.EQ.2) THEN
                  MC=IA+(KK-1)*KK/2-1
                ELSEIF (ISYM1.EQ.3) THEN
                  MC=IA+(KK-2)*(KK-1)/2-1
                ELSE
                  MC=IA+M1*(KK-1)-1
                ENDIF
              ELSE
C               IA must be changed, addresses of the ends of the columns
C               must be moved:
                IAN=IRAM(IA+KK-1)-M2-1
                DO 12, I1=M2,0,-1
                  IF (I1.GE.KK-1) THEN
                    IRAM(IAN+I1)=IRAM(IA+I1)
                  ELSE
                    IRAM(IAN+I1)=IRAM(IA+KK-1)
                  ENDIF
  12            CONTINUE
                IA=IAN
                MC=IA-1
              ENDIF
              IF (IC-1+NCTMP+NCCOL.GT.MC) THEN
C               Changing the written part of C to sparse matrix:
                CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,MINRAM,IA-1,
     *                      IC,NCTMP,KK-1,ICTMP1)
              ENDIF
              IF (IC-1+NCTMP+NCCOL.GT.MC) THEN
C               MATMUL-06
                WRITE(TXTERR,'(A,I9,A)')
     *         'MATMUL-06: Array RAM too small,',IC-1+NCTMP+NCCOL-MC,
     *         ' units missing.'
                CALL ERROR(TXTERR)
              ENDIF
            ENDIF
C           Initiating the values in column KK of C:
            DO 11, I1=IC-1+NCTMP+1,IC-1+NCTMP+NCCOL
              RAM(I1)=0.
  11        CONTINUE
C           Preparing loop over column KK of A:
            CALL MMILOC(M1,M2,ISYM1,NSPAR1,IA,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,IA,KK,IIA,II1,II3,II,AA)
              IF     ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN
                IIKK=ICTMP1+II
              ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN
                IIKK=ICTMP1+(KK-1)*KK/2+II
              ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN
                IIKK=ICTMP1+(KK-2)*(KK-1)/2+II
              ELSEIF (ISYMC.EQ.4)                  THEN
                IIKK=ICTMP1+(KK-1)*M1+II
              ELSE
                GOTO 13
              ENDIF
C             Storing the element:
              RAM(IIKK)=AA
              IF (((ISYM1.EQ.2).OR.(ISYM1.EQ.3)).AND.(II.NE.KK).AND.
     *             (ISYMC.EQ.4)) THEN
C               A is (anti)symmetric, storing also Aki=(-)Aik:
                IIKK=ICTMP1+(II-1)*M1+KK
                IF (ISYM1.EQ.3) AA=-AA
                RAM(IIKK)=AA
              ENDIF
  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,MINRAM,MAXRAM,
     *                  IC,NCTMP,M3,ICTMP1)
          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=M3+1+2*NEL2
      ELSE
        NB=NEL2
      ENDIF
      IB=IA-NB
      IF (NB+1.GT.MAXRAM-MINRAM-NA) THEN
C       MATMUL-09
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATMUL-09: Array RAM too small,',NB+1-(MAXRAM-MINRAM-NA),
     *  ' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
      CALL RMATD(LU1,FILED2,M3,SPARS2,NEL2,FORM2,IB)
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,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,RSPAR2)
      ELSE
        RSPAR2=FLOAT(NSPAR2)/FLOAT(NSPAR2+NEL2)
      ENDIF
C     Transposing input matrix B:
      IF (MATT2.EQ.1) THEN
        CALL TMATR(M2B,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB)
        CALL MSHIFT(M2B,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,IA-NB)
      ENDIF
C     Extending sparse matrix B:
      IF ((NSPAR1.GE.0).AND.((SYM1.EQ.'sym').OR.(SYM1.EQ.'skew')).AND.
     *    ((SYM2.EQ.'sym').OR.(SYM2.EQ.'skew'))) THEN
        CALL SMATRE(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB)
        SYM2=' '
      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 IC for the case that B will be sparse:
        CALL OMULS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,
     *             IB,NB,IC)
        IAN=MINRAM+IA-IC
C       Shifting IC if A is to be changed to sparse:
        IF (NSPAR1.LT.0) THEN
          IC=IC + NA - (2*NINT((1.-RSPAR1)*FLOAT(NA))+M2+1)
        ENDIF
        IF (IC.GE.MINRAM) THEN
C         Matrices may be (and will be) changed to sparse matrices
C         and shifted to optimum position for multiplication:
          IF (NSPAR2.LT.0) THEN
C           Changing B to sparse:
            CALL GSMATR(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB)
          ENDIF
C         Shifting B to the beginning of RAM:
          CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,MINRAM)
          IF (NSPAR1.LT.0) THEN
C           Changing A to sparse:
            CALL GSMATR(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA)
          ENDIF
C         Shifting A to the optimum position for multiplication:
          CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA,IAN)
C         Shifting B to the optimum position for multiplication:
          CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,IA-NB)
        ENDIF
      ELSE
C       Matrices are to be multiplied as non-sparse matrices:
C       Checking the memory needed for conversion of matrices
C       into non-sparse form:
        NTMP1=0
        NTMP2=0
        IF (NSPAR1.GE.0) NTMP1=M1
        IF (NSPAR2.GE.0) NTMP2=M2
        IF ((NELMAT(M1,M2,ISYM1)+NELMAT(M2,M3,ISYM2)+MAX0(NTMP1,NTMP2))
     *      .LE.MAXRAM) THEN
C         There is enough RAM for both matrices non-sparse plus the
C         auxiliary array to change them to non-sparse.
C         Calculating IC for the case that B will be non-sparse:
          CALL OMULNS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,
     *                IB,NB,IC)
          IAN=MINRAM+IA-IC
C         Shifting IC if A is to be changed to non-sparse:
          IF (NSPAR1.GE.0) THEN
            IC=IC + NA - NELMAT(M1,M2,ISYM1)
          ENDIF
C         Shifting B to the beginning of RAM:
          CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,MAXRAM,
     *                IB,NB,MINRAM)
          IF (NSPAR2.GE.0) THEN
C           Changing B to non-sparse:
            CALL SGMATR(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB)
          ENDIF
          IF (NSPAR1.GE.0) THEN
C           Shifting A behind B:
            CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,
     *                  IA,NA,IB+NB)
C           Changing A to non-sparse:
            CALL SGMATR(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA)
          ENDIF
C         Shifting A to the optimum position for multiplication:
          IF (IC.LT.MINRAM) IAN=MAXRAM-NA+1
          CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA,IAN)
C         Shifting B to the optimum position for multiplication:
          CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,IA-NB)
        ENDIF
      ENDIF
C
C     Multiplication:
      WRITE(*,'(A)') '+MATMUL: Multiplying...            '
C
      IF ((NSPAR1.LT.0).AND.(NSPAR2.LT.0).AND.
ccc  *    ((M1.LT.IB).AND.(M1*M3.LT.IA-M2))) THEN
     *    (IC.GE.MINRAM)) THEN
C       Multiplication for non-sparse matrices:
        IC=MINRAM
        NC=M1*M3
        NSPARC=-1
        NELC=NC
        ISYMC=ISYM3
        IIKK=IC-1
        II1=1
        II2=M1
        IA1=IA-1
        IB1=IB-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(IA1+JJ)
                ELSE
                  GOTO 25
                ENDIF
              ELSEIF (ISYM1.EQ.2) THEN
                IF (II.LE.JJ) THEN
                  AA=RAM(IA1+(JJ-1)*JJ/2+II)
                ELSE
                  AA=RAM(IA1+(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(IA1+(JJ-1)*(JJ-2)/2+II)
                ELSE
                  AA=-RAM(IA1+(II-1)*(II-2)/2+JJ)
                ENDIF
              ELSE
                AA=RAM(IA1+(JJ-1)*M1+II)
              ENDIF
C             Element of the matrix B:
              IF     (ISYM2.EQ.1) THEN
                IF (JJ.EQ.KK) THEN
                  BB=RAM(IB1+KK)
                ELSE
                  GOTO 25
                ENDIF
              ELSEIF (ISYM2.EQ.2) THEN
                IF (JJ.LE.KK) THEN
                  BB=RAM(IB1+(KK-1)*KK/2+JJ)
                ELSE
                  BB=RAM(IB1+(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(IB1+(KK-1)*(KK-2)/2+JJ)
                ELSE
                  BB=-RAM(IB1+(JJ-1)*(JJ-2)/2+KK)
                ENDIF
              ELSE
                BB=RAM(IB1+(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 IC is negative (no space for whole matrix C):
        IC=MINRAM
        LSPARC=.TRUE.
        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=IB-1
        NSPARC=-1
        NELC=NC
        ISYMC=ISYM3
C
        NCTMP=0
        ICTMP1=IC-1
        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 (IC-1+NCTMP+NCCOL.GT.MC) THEN
C           Erasing unneeded columns 1 to KK-1 of B:
            IF (NSPAR2.LT.0) THEN
C             IB stays unchanged, MC is redefined:
              IF     (ISYM2.EQ.1) THEN
                MC=IB+KK-1-1
              ELSEIF (ISYM2.EQ.2) THEN
                MC=IB+(KK-1)*KK/2-1
              ELSEIF (ISYM2.EQ.3) THEN
                MC=IB+(KK-2)*(KK-1)/2-1
              ELSE
                MC=IB+M2*(KK-1)-1
              ENDIF
            ELSE
C             IB must be changed, addresses of the ends of the columns
C             must be moved:
              IBN=IRAM(IB-1+KK)-(M3+1)
              DO 28, I1=M3,0,-1
                IF (I1.GE.KK-1) THEN
                  IRAM(IBN+I1)=IRAM(IB+I1)
                ELSE
                  IRAM(IBN+I1)=IRAM(IB-1+KK)
                ENDIF
  28          CONTINUE
              IB=IBN
              MC=IB-1
            ENDIF
            IF ((IC-1+NCTMP+NCCOL.GT.MC).AND.(LSPARC)) THEN
C             Changing the calculated part of C to sparse matrix:
              CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,MINRAM,IB-1,
     *                    IC,NCTMP,KK-1,ICTMP1)
            ENDIF
            IF (IC-1+NCTMP+NCCOL.GT.MC) THEN
C             MATMUL-10
              WRITE(TXTERR,'(A,I9,A)')
     *       'MATMUL-10: Array RAM too small,',IC-1+NCTMP+NCCOL-MC,
     *       ' units missing.'
              CALL ERROR(TXTERR)
            ENDIF
          ENDIF
C         Initiating the values in column KK of C:
          DO 29, I1=IC-1+NCTMP+1,IC-1+NCTMP+NCCOL
            RAM(I1)=0.
  29      CONTINUE
C         Preparing loop over column KK of B:
          CALL MMILOC(M2,M3,ISYM2,NSPAR2,IB,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,IB,KK,JJA,JJ1,JJ3,JJ,BB)
C           Preparing loop over column JJ of A:
            CALL MMILOC(M1,M2,ISYM1,NSPAR1,IA,JJ,II1,II2,II3)
C           Loop over nonzero elements of column JJ of A:
            DO 32, IIA=II1,II2,II3
C             Computing index II of row and value AA of element of A:
              CALL MMIRVE(M1,M2,ISYM1,NSPAR1,IA,JJ,IIA,II1,II3,II,AA)
              IF     ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN
                IIKK=ICTMP1+II
              ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN
                IIKK=ICTMP1+(KK-1)*KK/2+II
              ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN
                IIKK=ICTMP1+(KK-2)*(KK-1)/2+II
              ELSEIF (ISYMC.EQ.4)                  THEN
                IIKK=ICTMP1+(KK-1)*M1+II
              ELSE
                GOTO 30
              ENDIF
C             Multiplying Aij*Bjk:
              RAM(IIKK)=RAM(IIKK) + AA*BB
  30          CONTINUE
              IF (((ISYM1.EQ.2).OR.(ISYM1.EQ.3)).AND.(II.NE.JJ)) THEN
C               A is (anti)symmetric, multiplying also Aji=(-)Aij:
                IF     ((ISYMC.EQ.1).AND.(JJ.EQ.KK)) THEN
                  JJKK=ICTMP1+JJ
                ELSEIF ((ISYMC.EQ.2).AND.(JJ.LE.KK)) THEN
                  JJKK=ICTMP1+(KK-1)*KK/2+JJ
                ELSEIF ((ISYMC.EQ.3).AND.(JJ.LT.KK)) THEN
                  JJKK=ICTMP1+(KK-2)*(KK-1)/2+JJ
                ELSEIF (ISYMC.EQ.4)                  THEN
                  JJKK=ICTMP1+(KK-1)*M1+JJ
                ELSE
                  GOTO 32
                ENDIF
                IF (ISYM1.EQ.3) AA=-AA
                BBB=VELEM(M2,M3,ISYM2,NSPAR2,IB,II,KK)
C               Multiplying Aji*Bik:
                RAM(JJKK)=RAM(JJKK) + AA*BBB
              ENDIF
  32        CONTINUE
            IF (((ISYM2.EQ.2).OR.(ISYM2.EQ.3)).AND.(JJ.NE.KK)) THEN
C             B is (anti)symmetric, multiplying also Bkj=(-)Bjk:
              IF (ISYM2.EQ.3) BB=-BB
C             In this case matrix C may not be changed to sparse:
              LSPARC=.FALSE.
C             Preparing loop over column KK of A:
              CALL MMILOC(M1,M2,ISYM1,NSPAR1,IA,KK,II1,II2,II3)
C             Loop over nonzero elements stored for column KK of A:
              DO 34, IIA=II1,II2,II3
C               Computing index II of row and value AA of element of A:
                CALL MMIRVE(M1,M2,ISYM1,NSPAR1,IA,KK,IIA,II1,II3,II,AA)
                IF     ((ISYMC.EQ.1).AND.(II.EQ.JJ)) THEN
                  IIJJ=ICTMP1+II
                ELSEIF ((ISYMC.EQ.2).AND.(II.LE.JJ)) THEN
                  IIJJ=ICTMP1+(JJ-1)*JJ/2+II
                ELSEIF ((ISYMC.EQ.3).AND.(II.LT.JJ)) THEN
                  IIJJ=ICTMP1+(JJ-2)*(JJ-1)/2+II
                ELSEIF (ISYMC.EQ.4)                  THEN
                  IIJJ=ICTMP1+(JJ-1)*M1+II
                ELSE
                  GOTO 34
                ENDIF
C               Multiplying Aik*Bkj:
                RAM(IIJJ)=RAM(IIJJ) + AA*BB
  34          CONTINUE
C             If A is (anti)symmetric dense, multiplying the rest
C             of column KK of A under the diagonal.
C             (Here B is (anti)symmetric and thus A cannot be
C                                            (anti)symmetric sparse.)
              IF ((ISYM1.EQ.2).OR.(ISYM1.EQ.3)) THEN
                DO 36, II=KK+1,M1
                  IF     ((ISYMC.EQ.1).AND.(II.EQ.JJ)) THEN
                    IIJJ=ICTMP1+II
                  ELSEIF ((ISYMC.EQ.2).AND.(II.LE.JJ)) THEN
                    IIJJ=ICTMP1+(JJ-1)*JJ/2+II
                  ELSEIF ((ISYMC.EQ.3).AND.(II.LT.JJ)) THEN
                    IIJJ=ICTMP1+(JJ-2)*(JJ-1)/2+II
                  ELSEIF (ISYMC.EQ.4)                  THEN
                    IIJJ=ICTMP1+(JJ-1)*M1+II
                  ELSE
                    GOTO 36
                  ENDIF
                  IF (ISYM1.EQ.2) THEN
C                   Aik=Aki
                    AA=RAM(IA-1+(II-1)*II/2+KK)
C                   Multiplying Aik*Bkj:
                    RAM(IIJJ)=RAM(IIJJ) + AA*BB
                  ELSE
C                   Aik=-Aki
                    AA=-RAM(IA-1+(II-1)*(II-2)/2+KK)
C                   Multiplying Aik*Bkj:
                    RAM(IIJJ)=RAM(IIJJ) + AA*BB
                  ENDIF
  36            CONTINUE
              ENDIF
            ENDIF
  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,MINRAM,MAXRAM,
     *                IC,NCTMP,M3,ICTMP1)
        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,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC)
      ELSEIF ((ISPAR3.EQ.-1).AND.(NSPARC.GE.0)) THEN
C       Changing sparse matrix to non-sparse matrix:
        CALL SGMATR(M1,M3,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC)
      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,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC,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,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC)
          ENDIF
        ENDIF
      ENDIF
C
C     Writing output matrix C:
      SPARSC=' '
      IF (NSPARC.GE.0) SPARSC='CSC'
      CALL WMATH(LU1,FILE3,FILED3,M1,M3,SPARSC,NELC,SYM3,FORM3)
      CALL WMATD(LU1,FILED3,M1,M3,SPARSC,NELC,FORM3,IC)
C
      WRITE(*,'(A)') '+MATMUL: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE MMILOC(M1,M2,ISYM,NSPARS,IA,ICOLUM,
     *                 IELEM1,IELEM2,IELEM3)
      INTEGER M1,M2,ISYM,NSPARS,IA,ICOLUM,IELEM1,IELEM2,IELEM3
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     IA...   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
C             ICOLUM in array ARRAY.
C     IELEM2... Address of the value of the last element of column
C             ICOLUM in array ARRAY.
C     IELEM3... Step between two consecutive values of the matrix
C             elements in array ARRAY.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix.
C.......................................................................
C
C     Preparing loop over column ICOLUM of the matrix:
      IF (NSPARS.GE.0) THEN
        IELEM1=IRAM(IA-1+ICOLUM)+1
        IELEM2=IRAM(IA+ICOLUM)-1
        IELEM3=2
      ELSE
        IF     (ISYM.EQ.1) THEN
          IELEM1=IA+ICOLUM-1
          IELEM2=IELEM1
          IELEM3=1
        ELSEIF (ISYM.EQ.2) THEN
          IELEM1=IA+(ICOLUM-1)*ICOLUM/2
          IELEM2=IELEM1+M1-1
          IELEM3=1
        ELSEIF (ISYM.EQ.3) THEN
          IELEM1=IA+(ICOLUM-1)*(ICOLUM-2)/2
          IELEM2=IELEM1+M1-1
          IELEM3=1
        ELSE
          IELEM1=IA+(ICOLUM-1)*M1
          IELEM2=IELEM1+M1-1
          IELEM3=1
        ENDIF
      ENDIF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE MMIRVE(M1,M2,ISYM,NSPARS,IA,ICOLUM,
     *                 IELEMA,IELEM1,IELEM3,
     *                                     IROW,VELEM)
      INTEGER M1,M2,ISYM,NSPARS,IA,ICOLUM,IELEMA,IELEM1,IELEM3,IROW
      REAL VELEM
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     IA...   Address of the first storage location of the matrix.
C     ICOLUM... Number of the column under consideration.
C     IELEMA... Address of the value of the current element of column
C             ICOLUM in array RAM.
C     IELEM1... Address of the value of the first element of column
C             ICOLUM in array RAM.
C     IELEM3... Step between two consecutive values of the matrix
C             elements in array RAM.
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
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix.
C.......................................................................
C
C     Calculating index IROW of row of the matrix:
      IF (NSPARS.GE.0) THEN
        IROW=IRAM(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=RAM(IELEMA)
      ELSE
        IF     (ISYM.EQ.2) THEN
C         'sym'
          IF (IROW.LE.ICOLUM) THEN
            VELEM=RAM(IELEMA)
          ELSE
            VELEM=RAM(IA-1+(IROW-1)*IROW/2+ICOLUM)
          ENDIF
        ELSEIF (ISYM.EQ.3) THEN
C         'skew'
          IF (IROW.LT.ICOLUM) THEN
            VELEM=RAM(IELEMA)
          ELSEIF (IROW.EQ.ICOLUM) THEN
            VELEM=0.
          ELSE
            VELEM=-RAM(IA-1+(IROW-1)*(IROW-2)/2+ICOLUM)
          ENDIF
        ELSE
C         'diag' or ' '
          VELEM=RAM(IELEMA)
        ENDIF
      ENDIF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OMULS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,
     *                 MIA,MAA,IA,NA,IC)
      INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IC
C
C Subroutine to estimate optimum origin IC 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 IC 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     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.
C     IA...   Address of the first storage location of matrix B.
C     NA...   Number of storage locations of matrix B.
C     For detailed description of storage of matrices in the memory
C     refer to file mat.for.
C
C Output:
C     IC...   Address of the origin where to write matrix C.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NELMAT
      INTEGER NELMAT
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER II,KK,JJ,NN,JB,II1,II2,ITMP
C
C.......................................................................
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:
        ITMP=IA-M2-1
        IF (ITMP.LE.MIA) THEN
C         There is no space for ends of rows, thus IC is put to -1:
          IC=-1
          RETURN
        ENDIF
        JB=IA+NA-1
        JJ=JB
        IRAM(ITMP+M2)=JJ+1
        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 (RAM(JB).NE.0.) NN=NN+1
            JB=JB-1
  10      CONTINUE
          JJ=JJ-2*NN
          IRAM(ITMP+KK-1)=JJ+1
  12    CONTINUE
      ELSE
        ITMP=IA
      ENDIF
      IC=ITMP
      DO 20, KK=1,M2
        IC=MIN0(IC,IRAM(ITMP-1+KK)-M2-1-NELMAT(M1C,KK,ISYMC))
  20  CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OMULNS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,
     *                  MIA,MAA,IA,NA,IC)
      INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IC
C
C Subroutine to estimate optimum origin IC where to write matrix C,
C assuming that matrix B given to OMULNS 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 IC 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     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.
C     IA...   Address of the first storage location of matrix B.
C     NA...   Number of storage locations of matrix B.
C     For detailed description of storage of matrices in the memory
C     refer to file mat.for.
C
C Output:
C     IC...   Origin where to write matrix C.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL NELMAT
      INTEGER NELMAT
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER KK,ITMP
C
C.......................................................................
C
      IF (NSPAR.LT.0) THEN
        ITMP=IA
      ELSE
        ITMP=IA+NA-NELMAT(M1,M2,ISYM)
      ENDIF
      IC=ITMP
      DO 20, KK=1,M2
        IC=MIN0(IC,ITMP+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,MIA,MAA,IA,NA,ICOL,ITMP)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,ICOL,ITMP
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 RAM, and this part is changed to sparse matrix and written
C to RAM and IRAM together with M2+1 addresses of beginnings
C of columns. For second and other invocations, it is assumed that
C IRAM and RAM contain M2+1 addresses of beginnings 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+1 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             SYM='diag' ... ISYM=1
C             SYM='sym'  ... ISYM=2
C             SYM='skew' ... ISYM=3
C             SYM=' '    ... ISYM=4
C     NSPAR... Sparseness of the matrix.
C     NELEM... Number of elements of the matrix stored in array RAM.
C     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.  Entire arrays (I)RAM
C             from MIA to MAA may be used for temporary storage.
C     IA...   Address of the first storage location in array RAM
C             used for the matrix.
C     NA...   Number of storage locations for the input 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     NA...   New number of storage locations for the M1*ICOL matrix.
C     ITMP... Fictitious position in RAM 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     http://sw3d.cz/staff/bulant.htm
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NSPMAT,NELMAT
      INTEGER NSPMAT,NELMAT
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,I11,I12,IAR,IAS,NAR,NSPA,NELE,ICOL0,ISHIFT
      CHARACTER*72 TXTERR
      SAVE ICOL0
      DATA ICOL0/0/
C
C.......................................................................
C
C     Beginning IAR and number of elements NAR of the non-sparse part
C     of RAM:
      IF (ICOL0.EQ.0) THEN
        IAR=IA
        NAR=NA
      ELSE
        IAR=IRAM(IA+M2)
        NAR=NA-(IAR-IA)
      ENDIF
C     Number NSPA of zero and NELE of nonzero elements:
      NSPA=NSPMAT(NAR,RAM(IAR))
      NELE=NAR-NSPA
      IF (((ICOL0.EQ.0).AND.(IAR+M2+1+2*NELE.GT.MAA)).OR.
     *    ((ICOL0.NE.0).AND.(IAR+     2*NELE.GT.MAA))) THEN
C       MATMUL-12
        IF (ICOL0.EQ.0) THEN
          I1=IAR+M2+1+2*NELE - MAA
        ELSE
          I1=IAR+     2*NELE - MAA
        ENDIF
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATMUL-12: Array RAM too small,',I1,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     Moving the dense part of the matrix to the end of available RAM:
      DO 10, I1=0,NAR-1
        RAM(MAA-I1)=RAM(IAR+NAR-1-I1)
  10  CONTINUE
      IAR=MAA-(NAR-1)
C     Beginning IAS where to write the new sparse part of the matrix,
C     shifting index ISHIFT for case of ICOL0.eq.0:
      IF (ICOL0.EQ.0) THEN
        IAS=IA+ICOL+1
        ISHIFT=M2+1-(ICOL+1)
        IRAM(IA)=IAS+ISHIFT
      ELSE
        IAS=IRAM(IA+M2)
        ISHIFT=0
      ENDIF
C     Moving the nonzero elements of the matrix,
C     storing their row indices, recording column indices:
      I3=IAR
      I4=IAS
C     I3 points to the dense array,
C     I4 points to the new part of the sparse array
      DO 20, I2=ICOL0+1,ICOL
        IF     (ISYM.EQ.1) THEN
          I12=I2
          I11=I2
        ELSEIF (ISYM.EQ.2) THEN
          I12=I2
          I11=1
        ELSEIF (ISYM.EQ.3) THEN
          I12=I2-1
          I11=1
        ELSE
          I12=M1
          I11=1
        ENDIF
        DO 18, I1=I11,I12
          IF (RAM(I3).NE.0.) THEN
            IF (I4+1.LE.I3) THEN
C             MATMUL-13
              CALL ERROR('MATMUL-13: Array RAM too small')
C             The output matrix is being formed in the memory,
C             but there was not enough space in RAM for it.
C             Thus the already formed part of the matrix is being
C             converted to the sparse matrix. There is enough space for
C             the sparse matrix, but the conversion is not possible.
C             May be that only several more memory locations will solve
C             the problem, try to enlarge MRAM, if possible.
            ENDIF
            IRAM(I4)=I1
            RAM(I4+1)=RAM(I3)
            I4=I4+2
          ENDIF
          I3=I3+1
  18    CONTINUE
        IRAM(IA-1+I2+1)=I4+ISHIFT
  20  CONTINUE
C
C     I4 now points just after the last nonzero element
      IF (ISHIFT.NE.0) THEN
C       Shifting the sparse part of the matrix:
        IF (I4+ISHIFT.GT.MAA) THEN
C         MATMUL-14
          WRITE(TXTERR,'(A,I9,A)')  'MATMUL-14: Array RAM too small,',
     *    I4+ISHIFT-MAA,' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        DO 30, I1=I4-2,IAS,-2
          RAM(I1+1+ISHIFT)=RAM(I1+1)
          IRAM(I1+ISHIFT)=IRAM(I1)
  30    CONTINUE
      ENDIF
C
C     Recording pointers for columns ICOL+2 to M2+1:
      DO 32, I1=IA-1+ICOL+2,IA-1+M2+1
        IRAM(I1)=IRAM(IA-1+ICOL+1)
  32  CONTINUE
C
C
      ICOL0=ICOL
      NA=IRAM(IA+M2)-IA
      NELEM=(NA-M2-1)/2
      I1=NELMAT(M1,ICOL,ISYM)
      NSPAR=I1-NELEM
      ITMP=IRAM(IA+M2)-1-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
C
C=======================================================================
C