Example 21. General Out-of-core Matrix Multiply

This example demonstrates the use of the FMS general matrix multiply routine FMSMM. This example is designed to test the FMSMM subroutine.
C       E X A M P L E   21
C
C	Program name:
	CHARACTER*10 MYNAME
	PARAMETER (MYNAME='EXAMPLE_21')
C
C	FMS Parameters:
	INTEGER IACCUM
C
C	FMS matrix and vector file attributes:
C
C	FMS memory management requires the following arrays:
	POINTER (CMD_PTR, CMD)
	POINTER (RMD_PTR, RMD)
	POINTER (IMD_PTR, IMD)
	COMPLEX*16 CMD(0:1)
	REAL*8     RMD(0:1)
	INTEGER    IMD(0:1)
C
C	Input data:
	LOGICAL ASK
	INTEGER ASK_I
C
C	Profile vector for a full matrix:
	INTEGER LOWEQ(1)
C
C	Local variables:
	INTEGER     M, N, K
	INTEGER     IDTYPE, I, LX, LDISK, NV
	INTEGER     MBA, MBB, MBC, IORATE
	REAL*8      ERROR
	REAL*8      RALPHA(1)
	REAL*8      SIZE, TEMP
	REAL*8      TCPU_A1, TWALL_A1, TIO_A1
	REAL*8      TCPU_A2, TWALL_A2, TIO_A2
	REAL*8      TCPU_B1, TWALL_B1, TIO_B1
	REAL*8      TCPU_B2, TWALL_B2, TIO_B2
	COMPLEX*16  CALPHA(1)
	CHARACTER*1 TRANSA
	LOGICAL     IFILL
	DATA LOWEQ/-1/
	DATA LU0(1)  /0/
	DATA RALPHA(1)/1.0D0/
	DATA CALPHA(1)/(1.0D0,1.0D0)/
C
C	Common block to communicate with fill routines:
	COMMON /MYDATA/ERROR, IOPT
C
C (1)   Initialize FMS:
	CALL FMSINI
	CALL FMSPSH(MYNAME)
	CALL FMSIST ('MDATAU', 1)
	CALL FMSIST ('LOGTIM',    3)
	CALL FMSIST ('IPRMV' , 1026)
	CALL FMSIST ('IPRVV' , 1026)
    1	CONTINUE
	WRITE (6,*) 'The FMS modules are numbered as follows:'
	WRITE (6,*) '   1 = Real Symmetric'
	WRITE (6,*) '   2 = Real Nonsymmetric'
	WRITE (6,*) '   3 = Complex Hermitian'
	WRITE (6,*) '   4 = Complex Symmetric'
	WRITE (6,*) '   5 = Complex Nonsymmetric'
	MOD = ASK_I('Enter the FMS module number (1 to 5)')
	IF( (MOD.LT.1) .OR. (MOD.GT.5) ) THEN
	   PRINT *,'Illegal FMS Module entered.  Try again.'
	   GO TO 1
	END IF
    2	CONTINUE
	WRITE (6,*) 'The following cases are available:'
	WRITE (6,*) '   1 = {C} = {D} + s*[A]{B}'
	WRITE (6,*) '   3 = {C} = {D} + s*{A}{B}'
	WRITE (6,*) '   4 = {C} = {D} + s*{A}T{B}'
	WRITE (6,*) '   5 = [C] = [D] + s*{A}T{B}'
	ICASE = ASK_I('Enter the case number (1 to 5)')
	IF( (ICASE.LT.1) .OR. (ICASE.GT.5) ) THEN
	   PRINT *,'Illegal case entered.  Try again.'
	   GO TO 2
	END IF
	M = ASK_I('Enter the number of rows    in C')
	IF(ICASE .NE. 5) THEN
	   N = ASK_I('Enter the number of columns in C')
	ELSE
	   N = M
	END IF
	IF( (ICASE .NE. 1)   .AND.
     1	    (ICASE .NE. 2) ) THEN
	   K = ASK_I('Enter the number of rows    in {B}')
	ELSE
	   K = M
	END IF
	IF(ICASE .EQ. 1) THEN
	   WRITE (6,2001) MOD, M, N, M, M, M, N
	ELSE IF (ICASE .EQ. 3) THEN
	   WRITE (6,2003) MOD, M, N, M, K, K, N
	ELSE IF (ICASE .EQ. 4) THEN
	   WRITE (6,2004) MOD, M, N, K, M, K, N
	ELSE
	   WRITE (6,2005) MOD, M, N, K, M, K, N
	END IF
	IFILL = ASK('Do you want to fill')
	WRITE (6,*) 'You may now alter any FMS parameter.'
	WRITE (6,*) 'When you are finished, type the letters RETURN'
	CALL FMSSET
	CALL FMSIGT ('MEMPTR', IMD_PTR)
	CALL FMSIGT ('MEMPTR', RMD_PTR)
	CALL FMSIGT ('MEMPTR', CMD_PTR)
	IF(MOD .LE. 2) THEN
	   IDTYPE = 1
	ELSE
	   IDTYPE = 2
	END IF
        IF(ICASE .LT. 4) THEN
           TRANSA = 'N'
        ELSE
           TRANSA = 'T'
        END IF
C (2)   Open FMS files:
C
C	LUA:
	IF(ICASE .EQ. 1) THEN
C	   This is a matrix file:
	   WRITE(6,*) 'Opening matrix file [A]:'
	   IF(MOD.EQ.1) CALL RSDI (LOWEQ, M, 'LUA', LUA)
	   IF(MOD.EQ.2) CALL RNDI (LOWEQ, M, 'LUA', LUA)
	   IF(MOD.EQ.3) CALL CHDI (LOWEQ, M, 'LUA', LUA)
	   IF(MOD.EQ.4) CALL CSDI (LOWEQ, M, 'LUA', LUA)
	   IF(MOD.EQ.5) CALL CNDI (LOWEQ, M, 'LUA', LUA)
	   SIZE = DFLOAT(LUA(5)) * DFLOAT(LUA(10))
	   IF(LUA(11) .EQ. 2) SIZE = 2.0*SIZE
	   MBA = (SIZE + 131071.)/131072.
	   WRITE(6,*) '   Storage for [A] (Mbytes)=', MBA
	ELSE
C	   This is a vector file:
	   WRITE(6,*) 'Opening vector file {A}:'
	   IF(TRANSA .EQ. 'N') THEN
C	      {A} is not transposed:
	      CALL FMSOV (M, IDTYPE, K, 'LUA', LUA)
	   ELSE
C	      {A} is transposed:
	      CALL FMSOV (K, IDTYPE, M, 'LUA', LUA)
	   END IF
	   SIZE = DFLOAT(LUA(4)) * DFLOAT(LUA(6))
	   MBA = (SIZE + 131071.)/131072.
	   WRITE(6,*) '   Storage for {A} (Mbytes)=', MBA
	END IF
C
C	LUB:
	WRITE(6,*) 'Opening vector file {B}:'
	CALL FMSOV (K, IDTYPE, N, 'LUB', LUB)
	SIZE = DFLOAT(LUB(4)) * DFLOAT(LUB(6))
	MBB = (SIZE + 131071.)/131072.
	WRITE(6,*) '   Storage for {B} (Mbytes)=', MBB
C
C	LUC:
	IF(ICASE .EQ. 5) THEN
C	   This is a matrix file:
	   WRITE(6,*) 'Opening matrix file [C]:'
	   IF(MOD.EQ.1) CALL RSDI (LOWEQ, M, 'LUC', LUC)
	   IF(MOD.EQ.2) CALL RNDI (LOWEQ, M, 'LUC', LUC)
	   IF(MOD.EQ.3) CALL CHDI (LOWEQ, M, 'LUC', LUC)
	   IF(MOD.EQ.4) CALL CSDI (LOWEQ, M, 'LUC', LUC)
	   IF(MOD.EQ.5) CALL CSDN (LOWEQ, M, 'LUC', LUC)
	   SIZE = DFLOAT(LUC(5)) * DFLOAT(LUC(10))
	   IF(LUC(11) .EQ. 2) SIZE = 2.0*SIZE
	   MBC = (SIZE + 131071.)/131072.
	   WRITE(6,*) '   Storage for [C] (Mbytes)=', MBA
	ELSE
C	   This is a vector file:
	   WRITE(6,*) 'Opening vector file {C}:'
	   CALL FMSOV (M, IDTYPE, N, 'LUC', LUC)
	   SIZE = DFLOAT(LUC(4)) * DFLOAT(LUC(6))
	   MBC = (SIZE + 131071.)/131072.
	   WRITE(6,*) '   Storage for {C} (Mbytes)=', MBC
	END IF
C
C (3)   Write data to FMS files:
	IF(IFILL) THEN
C
C	{A}
	CALL FMSTIM (TCPU_A1, TWALL_A1, TIO_A1)
	IOPT = 0
	IF(ICASE .EQ. 1) THEN
C	   This is a matrix file:
	   WRITE(6,*) 'Writing data to matrix file [A]:'
	   IF(MOD.EQ.1) CALL RSDA (LU0, 0, LUA)
	   IF(MOD.EQ.2) CALL RNDA (LU0, 0, LUA)
	   IF(MOD.EQ.3) CALL CHDA (LU0, 0, LUA)
	   IF(MOD.EQ.4) CALL CSDA (LU0, 0, LUA)
	   IF(MOD.EQ.5) CALL CNDA (LU0, 0, LUA)
	ELSE
C	   This is a vector file:
	   WRITE(6,*) 'Writing data to vector file {A}:'
	   CALL VFILL (LUA)
	END IF
	CALL FMSTIM (TCPU_A2, TWALL_A2, TIO_A2)
	TEMP = TWALL_A2 - TWALL_A1
	IF(TEMP .LT. .01) TEMP=.01
	TEMP = DFLOAT(MBA)/TEMP
	IORATE = TEMP
	WRITE(6,*) '   Transfer rate  (Mb/Sec.)=', IORATE
C
C	{B}
	WRITE(6,*) 'Writing data to vector file {B}:'
	CALL FMSTIM (TCPU_B1, TWALL_B1, TIO_B1)
	CALL VFILL (LUB)
	END IF
	CALL FMSTIM (TCPU_B2, TWALL_B2, TIO_B2)
	TEMP = TWALL_B2 - TWALL_B1
	IF(TEMP .LT. .01) TEMP=.01
	TEMP = DFLOAT(MBB)/TEMP
	IORATE = TEMP
	WRITE(6,*) '   Transfer rate  (Mb/Sec.)=', IORATE
C
C (4)   Perform matrix algebra:
	WRITE(6,*) 'Calling FMSMM'
C	Compute {C} = {A}{B}:
	CALL FMSMM (LUC, 0, 0, TRANSA, LUA, LUB)
C	Compute {C} = {C} - {A}{B}:
	CALL FMSMM (LUC, LUC, 0, TRANSA, LUA, LUB)
C
C (5)   Read data from FMS files:
C       Check the answer; {C} should be {0}:
	IOPT  = 1
	ERROR = 0.0D0
	IF(ICASE .EQ. 5) THEN
C	   [C]is a matrix file:
	   WRITE(6,*) 'Checking matrix file [C]:'
	   IF(MOD.EQ.1) CALL RSDAF
     1	      (LUC, RALPHA, 1, LU0, 0, LU0, LU0, LU0, LU0, 0)
	   IF(MOD.EQ.2) CALL RNDAF (LU0, 0, LUC)
     1	      (LUC, RALPHA, 1, LU0, 0, LU0, LU0, LU0, LU0, 0)
	   IF(MOD.EQ.3) CALL CHDAF (LU0, 0, LUC)
     1	      (LUC, RALPHA, 1, LU0, 0, LU0, LU0, LU0, LU0, 0)
	   IF(MOD.EQ.4) CALL CSDAF (LU0, 0, LUC)
     1	      (LUC, CALPHA, 1, LU0, 0, LU0, LU0, LU0, LU0, 0)
	   IF(MOD.EQ.5) CALL CNDAF (LU0, 0, LUC)
     1	      (LUC, CALPHA, 1, LU0, 0, LU0, LU0, LU0, LU0, 0)
	ELSE
C	   {C} is a vector file:
	   WRITE(6,*) 'Checking vector file {C}:'
	   CALL VCHECK (LUC, ERROR)
	END IF
	WRITE(6,*) 'MAXIMUM ERROR =', ERROR
C
C (6)   Close FMS files:
C
C	{A}
	IF(ICASE .EQ. 1) THEN
	   CALL FMSCM (LUA)
	ELSE
	   CALL FMSCV (LUA)
	END IF
C
C	{B}
	CALL FMSCV (LUB)
C
C	{C}
	IF(ICASE .EQ. 5) THEN
	   CALL FMSCM (LUC)
	ELSE
	   CALL FMSCV (LUC)
	END IF
	IF(ASK('Do you want another solution?')) GO TO 1
	CALL FMSPOP(MYNAME)
	CALL FMSEND
 2001	FORMAT (/'You selected case 1, MOD',I2/
     1	 '   C{',I7,',',I7,'} = A[',I7,',',I7,'] * B{',I7,',',I7,'}')
 2003	FORMAT (/'You selected case 3, MOD',I2/
     1	 '   C{',I7,',',I7,'} = A{',I7,',',I7,'} * B{',I7,',',I7,'}')
 2004	FORMAT (/'You selected case 4, MOD',I2/
     1	 '   C{',I7,',',I7,'} = A{',I7,',',I7,'}T * B{',I7,',',I7,'}')
 2005	FORMAT (/'You selected case 5, MOD',I2/
     1	 '   C[',I7,',',I7,'] = A{',I7,',',I7,'}T * B{',I7,',',I7,'}')
	END
C=======================================================================
	LOGICAL FUNCTION ASK(QUESTION)
C=======================================================================
	CHARACTER* (*) QUESTION
	CHARACTER*1 IYN
	WRITE(6,2000) QUESTION
	READ (5,1000) IYN
	IF( (IYN .EQ. 'Y') .OR. (IYN .EQ. 'y') ) THEN
	   ASK = .TRUE.
	ELSE
	   ASK = .FALSE.
	END IF
	RETURN
 1000	FORMAT (A)
 2000	FORMAT (1X,A,' (y,n)>')
	END
C=======================================================================
	INTEGER FUNCTION ASK_I(STRING)
C=======================================================================
	CHARACTER* (*) STRING
	WRITE(6,2000) STRING
	READ (5,*) ASK_I
	RETURN
 2000	FORMAT (1X,A,'>')
	END
	
C	================================================================
	SUBROUTINE RSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
     1                     JCOL1, JCOL2, IJSTEP)
C	================================================================
	INTEGER    IROW1, IROW2, JCOL1, JCOL2, IJSTEP
	INTEGER    LOWEQ(*), LOCEQ(*)
	REAL*8     A(0:*), D(*)
	COMMON /MYDATA/ERROR, IOPT
	REAL*8 ERROR
	REAL*8 EI
C
	IF(IOPT .EQ. 0) THEN
C          Populate the diagonal with test data:
	   IF(IROW2 .EQ. JCOL2) THEN
C	      This is a diagonal block:
	      DO I = IROW1,IROW2
	         D(I) = RANAIJ(I,I)
	      END DO
	   END IF
C
C	   Populate [AL] with random numbers:
C	   The term A(I,J) is addressed as A(LOCEQ(I)+IJSTEP*J)
	   DO I = IROW1,IROW2
	      J1 = MAX0(LOWEQ(I),JCOL1)
	      J2 = MIN0(I-1,JCOL2)
	      L  = LOCEQ(I) + IJSTEP*J1
	      DO J = J1,J2
	         A(L) = RANAIJ(I,J)
	         L = L + IJSTEP
	      END DO
	   END DO
	ELSE
C	   Determine the maximum element in [AL]
	   IF(IROW2 .EQ. JCOL2) THEN
C	      This is a diagonal block:
	      DO I = IROW1,IROW2
	         EI = ABS(D(I))
	         IF(EI .GT. ERROR) ERROR = EI
	      END DO
	   END IF
	   DO I = IROW1,IROW2
	      J1 = MAX0(LOWEQ(I),JCOL1)
	      J2 = MIN0(I-1,JCOL2)
	      L  = LOCEQ(I) + IJSTEP*J1
	      DO J = J1,J2
	         EI = ABS(A(L))
	         IF(EI .GT. ERROR) ERROR = EI
	         L = L + IJSTEP
	      END DO
	   END DO
	END IF
	RETURN
	END
C	================================================================
	SUBROUTINE RNUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
     1                     JCOL1, JCOL2, IJSTEP)
C	================================================================
	INTEGER    IROW1, IROW2, JCOL1, JCOL2, IJSTEP
	INTEGER    LOWEQ(*), LOCEQ(*)
	REAL*8     A(0:*), D(*)
	COMMON /MYDATA/ERROR, IOPT
	REAL*8     ERROR
	REAL*8     EI
C
	IF(IOPT .EQ. 0) THEN
C	   Populate [AU] with test data:
C	   The term A(I,J) is addressed as A(LOCEQ(J)+IJSTEP*I)
	   DO J = JCOL1,JCOL2
	      I1 = MAX0(LOWEQ(J),IROW1)
	      I2 = MIN0(J-1,IROW2)
	      L  = LOCEQ(J) + IJSTEP*I1
	      DO I = I1,I2
	         A(L) = RANAIJ(I,J)
	         L = L + IJSTEP
	      END DO
	   END DO
	ELSE
C	   Find the largest element in [AU]:
	   DO J = JCOL1,JCOL2
	      I1 = MAX0(LOWEQ(J),IROW1)
	      I2 = MIN0(J-1,IROW2)
	      L  = LOCEQ(J) + IJSTEP*I1
	      DO I = I1,I2
	         EI = ABS(A(L))
	         IF(EI .GT. ERROR) ERROR = EI
	         L = L + IJSTEP
	      END DO
	   END DO
	END IF
	RETURN
	END
C	================================================================
	SUBROUTINE CHUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
     1                     JCOL1, JCOL2, IJSTEP)
C	================================================================
	INTEGER    IROW1, IROW2, JCOL1, JCOL2, IJSTEP
	INTEGER    LOWEQ(*), LOCEQ(*)
	REAL*8     D(*)
	COMPLEX*16 A(0:*)
	COMMON /MYDATA/ERROR, IOPT
	REAL*8     ERROR
	REAL*8     EI
	REAL*8     AIJ
C
	IF(IOPT .EQ. 0) THEN
C          Populate the diagonal with test data:
	   IF(IROW2 .EQ. JCOL2) THEN
C	      This is a diagonal block:
	      DO I = IROW1,IROW2
	         D(I) = RANAIJ(I,I)
	      END DO
	   END IF
C
C	   Populate profile of [AL] with test data:
C	   The term A(I,J) is addressed as A(LOCEQ(I)+IJSTEP*J)
	   DO I = IROW1,IROW2
	      J1 = MAX0(LOWEQ(I),JCOL1)
	      J2 = MIN0(I-1,JCOL2)
	      L  = LOCEQ(I) + IJSTEP*J1
	      DO J = J1,J2
	         AIJ  = RANAIJ(I,J)
	         A(L) = DCMPLX(AIJ,AIJ)
	         L = L + IJSTEP
	      END DO
	   END DO
	ELSE
C          Determine the maximum element in [AL]:
	   IF(IROW2 .EQ. JCOL2) THEN
C	      This is a diagonal block:
	      DO I = IROW1,IROW2
	         EI = ABS(D(I))
	         IF(EI .GT. ERROR) ERROR = EI
	      END DO
	   END IF
C
	   DO I = IROW1,IROW2
	      J1 = MAX0(LOWEQ(I),JCOL1)
	      J2 = MIN0(I-1,JCOL2)
	      L  = LOCEQ(I) + IJSTEP*J1
	      DO J = J1,J2
	         EI = ABS(A(L))
	         IF(EI .GT. ERROR) ERROR = EI
	         L = L + IJSTEP
	      END DO
	   END DO
	END IF
	RETURN
	END
C	================================================================
	SUBROUTINE CSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
     1                     JCOL1, JCOL2, IJSTEP)
C	================================================================
	INTEGER    IROW1, IROW2, JCOL1, JCOL2, IJSTEP
	INTEGER    LOWEQ(*), LOCEQ(*)
	COMPLEX*16 A(0:*), D(*)
	COMMON /MYDATA/ERROR, IOPT
	REAL*8     ERROR
	REAL*8     EI
	REAL*8     AIJ
C
	IF(IOPT .EQ. 0) THEN
C          Populate the diagonal with test data:
	   IF(IROW2 .EQ. JCOL2) THEN
C	      This is a diagonal block:
	      DO I = IROW1,IROW2
	         AIJ  = RANAIJ(I,I)
	         D(I) = DCMPLX(AIJ,AIJ)
	      END DO
	   END IF
C
C	   Populate profile of [AL] with test data:
C	   The term A(I,J) is addressed as A(LOCEQ(I)+IJSTEP*J)
	   DO I = IROW1,IROW2
	      J1 = MAX0(LOWEQ(I),JCOL1)
	      J2 = MIN0(I-1,JCOL2)
	      L  = LOCEQ(I) + IJSTEP*J1
	      DO J = J1,J2
	         AIJ  = RANAIJ(I,J)
	         A(L) = DCMPLX(AIJ,AIJ)
	         L = L + IJSTEP
	      END DO
	   END DO
	ELSE
C          Find the maximum element in [AL]:
	   IF(IROW2 .EQ. JCOL2) THEN
C	      This is a diagonal block:
	      DO I = IROW1,IROW2
	         EI = ABS(D(I))
	         IF(EI .GT. ERROR) ERROR = EI
	      END DO
	   END IF
C
	   DO I = IROW1,IROW2
	      J1 = MAX0(LOWEQ(I),JCOL1)
	      J2 = MIN0(I-1,JCOL2)
	      L  = LOCEQ(I) + IJSTEP*J1
	      DO J = J1,J2
	         EI = ABS(A(L))
	         IF(EI .GT. ERROR) ERROR = EI
	         L = L + IJSTEP
	      END DO
	   END DO
	END IF
	RETURN
	END
C	================================================================
	SUBROUTINE CNUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
     1                     JCOL1, JCOL2, IJSTEP)
C	================================================================
	INTEGER    IROW1, IROW2, JCOL1, JCOL2, IJSTEP
	INTEGER    LOWEQ(*), LOCEQ(*)
	COMPLEX*16 A(0:*), D(*)
	COMMON /MYDATA/ERROR, IOPT
	REAL*8     EI
	REAL*8     ERROR
	REAL*8     AIJ
C
	IF(IOPT .EQ. 0) THEN
C	   Populate profile of [AL] with test data:
C	   The term A(I,J) is addressed as A(LOCEQ(J)+IJSTEP*I)
	   DO J = JCOL1,JCOL2
	      I1 = MAX0(LOWEQ(J),IROW1)
	      I2 = MIN0(J-1,IROW2)
	      L  = LOCEQ(J) + IJSTEP*I1
	      DO I = I1,I2
	         AIJ  = RANAIJ(I,J)
	         A(L) = DCMPLX(AIJ,AIJ)
	         L = L + IJSTEP
	      END DO
	   END DO
	ELSE
C	   Determine the maximum element in [AU]:
C	   The term A(I,J) is addressed as A(LOCEQ(J)+IJSTEP*I)
	   DO J = JCOL1,JCOL2
	      I1 = MAX0(LOWEQ(J),IROW1)
	      I2 = MIN0(J-1,IROW2)
	      L  = LOCEQ(J) + IJSTEP*I1
	      DO I = I1,I2
	         EI = ABS(A(L))
	         IF(EI .GT. ERROR) ERROR = EI
	         L = L + IJSTEP
	      END DO
	   END DO
	END IF
	RETURN
	END
C	================================================================
	SUBROUTINE RNUSLB (A, LOCI, LOCJ, LUFLAG, JEQN1, JEQN2, NUMEQ)
C	================================================================
	INTEGER    JEQN1, JEQN2, NUMEQ
	INTEGER    LOCI(NUMEQ), LOCJ(JEQN1:JEQN2,2), LUFLAG(NUMEQ)
	REAL*8     A(0:*)
	COMMON /MYDATA/ERROR, IOPT
	REAL*8     ERROR
	REAL*8     EI
C
	IF(IOPT .EQ. 0) THEN
C          Fill [A] with random numbers:
	   DO J = JEQN1,JEQN2
	      DO I = 1,NUMEQ
	         L = LOCI(I) + LOCJ(J,LUFLAG(I))
	         A(L) = RANAIJ(I,J)
	      END DO
	   END DO
	ELSE
C          Find the maximum element in [A]:
	   DO J = JEQN1,JEQN2
	      DO I = 1,NUMEQ
	         L = LOCI(I) + LOCJ(J,LUFLAG(I))
	         EI = ABS(A(L))
	         IF(EI .GT. ERROR) ERROR = EI
	      END DO
	   END DO
	END IF
	RETURN
	END
C	================================================================
	SUBROUTINE CNUSLB (A, LOCI, LOCJ, LUFLAG, JEQN1, JEQN2, NUMEQ)
C	================================================================
	INTEGER    JEQN1, JEQN2, NUMEQ
	INTEGER    LOCI(NUMEQ), LOCJ(JEQN1:JEQN2,2), LUFLAG(NUMEQ)
	COMPLEX*16 A(0:*), DIA, OFFDIA
	PARAMETER  (DIA   =( 1.0D0, 1.0D0))
	PARAMETER  (OFFDIA=(-1.0D0,-1.0D0))
	COMMON /MYDATA/ERROR, IOPT
	REAL*8     ERROR
	REAL*8     EI
	REAL*8     AIJ
C
	IF(IOPT .EQ. 0) THEN
C          Fill [A] with random numbers:
	   DO J = JEQN1,JEQN2
	      DO I = 1,NUMEQ
	         L = LOCI(I) + LOCJ(J,LUFLAG(I))
	         AIJ  = RANAIJ(I,J)
	         A(L) = DCMPLX(AIJ,AIJ)
	      END DO
	   END DO
	ELSE
C          Find the maximum element in [A]:
	   DO J = JEQN1,JEQN2
	      DO I = 1,NUMEQ
	         L = LOCI(I) + LOCJ(J,LUFLAG(I))
	         EI = ABS(A(L))
	         IF(EI .GT. ERROR) ERROR = EI
	      END DO
	   END DO
	END IF
	RETURN
	END
C=======================================================================
	SUBROUTINE VFILL (LUX)
C=======================================================================
C	Fill vector file with random data
	INTEGER LUX(25)
	CHARACTER*5 MYNAME
	PARAMETER (MYNAME='VFILL')
	INTEGER MAXMD
	INTEGER MDUSED
	INTEGER MDAVIL
	REAL*8  RAVIL
	REAL*8  RREQD
C
C	FMS memory management requires the following arrays:
	POINTER (RMD_PTR, RMD)
	REAL*8     RMD(0:1)
	CALL FMSPSH(MYNAME)
	CALL FMSIGT ('MEMPTR', RMD_PTR)
C
	NUMEQ  = LUX(3)
	LENVEC = LUX(4)
	IDTYPE = LUX(5)
	NUMVEC = LUX(6)
C
C       Allocate storage for a block of vectors:
C
C	How much memory is available?
	CALL FMSIGT ('MEMPTR', MAXMD)
	CALL FMSIGT ('MEMPTR', MDUSED)
	MDAVIL = MAXMD - MDUSED
	RAVIL  = DFLOAT(MDAVIL)
	RREQD  = DFLOAT(LENVEC)*DFLOAT(NUMVEC)
	IF(RAVIL .GT. RREQD) THEN
C	   The memory available is larger than the vector file.
	   MDAVIL = LENVEC*NUMVEC
	END IF
C
C	Allocate MDAVIL words of memory:
	CALL FMSRMG (RMD, LV1, MDAVIL)
C
C	Fill blocks of vectors with random numbers and write to LUX:
	NVINC  = MDAVIL/LENVEC
	LDISK  = 1
	DO NV1 = 1,NUMVEC,NVINC
	   NV2 = MIN0(NV1+NVINC-1,NUMVEC)
	   NVEC = NV2 - NV1 + 1
	   LV  = LV1
	   DO NV = NV1,NV2
	      CALL XRAN (RMD(LV), IDTYPE, NUMEQ, NV)
	      LV = LV + LENVEC
	   END DO
	   NWRITE = NVEC*LENVEC
	   CALL FMSWRT (LUX(1), LDISK, RMD(LV1), NWRITE)
	   LDISK = LDISK + NWRITE
	END DO
C
C	Return the memory
	CALL FMSRMR (RMD, LV1, MDAVIL)
	CALL FMSPOP(MYNAME)
	RETURN
	END
C=======================================================================
	SUBROUTINE VCHECK (LUX, ERROR)
C=======================================================================
C	Find the largest absolute value in file LUX
	INTEGER LUX(25)
	REAL*8  ERROR
	REAL*8  EI
	CHARACTER*6 MYNAME
	PARAMETER (MYNAME='VCHECK')
	INTEGER MAXMD
	INTEGER MDUSED
	INTEGER MDAVIL
	REAL*8  RAVIL
	REAL*8  RREQD
C
C	FMS memory management requires the following arrays:
	POINTER (RMD_PTR, RMD)
	REAL*8     RMD(0:1)
	CALL FMSPSH(MYNAME)
	CALL FMSIGT ('MEMPTR', RMD_PTR)
C
	NUMEQ  = LUX(3)
	LENVEC = LUX(4)
	IDTYPE = LUX(5)
	NUMVEC = LUX(6)
C
C       Allocate storage for a block of vectors:
C
C	How much memory is available?
	CALL FMSIGT ('MEMPTR', MAXMD)
	CALL FMSIGT ('MEMPTR', MDUSED)
	MDAVIL = MAXMD - MDUSED
	RAVIL  = DFLOAT(MDAVIL)
	RREQD  = DFLOAT(LENVEC)*DFLOAT(NUMVEC)
	IF(RAVIL .GT. RREQD) THEN
C	   The memory available is larger than the vector file.
	   MDAVIL = LENVEC*NUMVEC
	END IF
C
C	Allocate MDAVIL words of memory:
	CALL FMSRMG (RMD, LV1, MDAVIL)
C
C	Read blocks of vectors:
	NVINC  = MDAVIL/LENVEC
	LDISK  = 1
	DO NV1 = 1,NUMVEC,NVINC
	   NV2 = MIN0(NV1+NVINC-1,NUMVEC)
	   NVEC = NV2 - NV1 + 1
	   NREAD= NVEC*LENVEC
	   CALL FMSRED (LUX(1), LDISK, RMD(LV1), NREAD)
	   LDISK = LDISK + NREAD
	   LV  = LV1
	   DO NV = NV1,NV2
	      LVI = LV
	      IMAX = IDTYPE*NUMEQ
	      DO I = 1,IMAX
	         EI = ABS(RMD(LVI))
	         IF(EI .GT. ERROR) ERROR = EI
	         LVI = LVI + 1
	      END DO
	      LV = LV + LENVEC
	   END DO
	END DO
C
C	Return the memory
	CALL FMSRMR (RMD, LV1, MDAVIL)
	CALL FMSPOP(MYNAME)
	RETURN
	END
C=======================================================================
	FUNCTION RANAIJ (I,J)
C=======================================================================
	REAL*8    RANAIJ
	INTEGER   I,J
	INTEGER   IY
	INTEGER   ISEEDR, ISEEDC
	PARAMETER (ISEEDR=11111111)
	PARAMETER (ISEEDC=12345678)
        IY     = I*ISEEDR + J*ISEEDC
        IY     = MOD(IY,65536)
	IY     = IY * 11111 + 1
	RANAIJ = DFLOAT((IY/256)+1) / 2844374.0D0
	RANAIJ = RANAIJ - 0.5D0
	RETURN
	END
C=======================================================================
	SUBROUTINE XRAN (X, IDTYPE, NROWS, NCOL)
C=======================================================================
	REAL*8    X(IDTYPE, NROWS)
	INTEGER   IDTYPE, NROWS, NCOL
	REAL*8    AIJ
	INTEGER   ISEEDR, ISEEDC
	PARAMETER (ISEEDR=11111111)
	PARAMETER (ISEEDC=12345678)
	DO IROW = 1,NROWS
           IY     = IROW*ISEEDR + NCOL*ISEEDC
           IY     = MOD(IY,65536)
	   IY     = IY * 11111 + 1
	   AIJ    = DFLOAT((IY/256)+1) / 2844374.0D0
	   AIJ    = AIJ - 0.5D0
	   X(1,IROW) = AIJ
	   IF(IDTYPE .EQ. 2) X(2,IROW) = AIJ
	END DO
	RETURN
	END