Example 21. General Out-of-core Matrix Multiply
This example demonstrates the use of the
FMS general
matrix multiply routine
FMSMM.
-
FMS Module
RS = Real Symmetric
RN = Real Nonsymmetric
CH = Complex Hermitian
CS = Complex Symmetric
CN = Complex Nonsymmetric
-
Size and type of matrices
-
Any FMS Parameter
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