Go To Top Go Up Go Back Go Forward

Example 11. General Out-of-core Matrix

This example solves a system of equations where the following may be specified as input:

This example is designed to support pre-sales benchmarking of machines.

C       E X A M P L E   11
C
C	Subroutines called in parallel:
	EXTERNAL WRITEV
C
C	Input Data:
	LOGICAL ASK
	INTEGER ASK_I
C
C	Number of vectors to reduce during factoring:
	PARAMETER (NUMRED = 0)
C
C	Skip operations during solving (no):
	PARAMETER (ISKIP  = 0)
C
C	FMS matrix and vector file attributes:
	INTEGER     LUF(25)
	INTEGER     LUX(25)
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	Profile vector for a full matrix:
	INTEGER LOWEQ(1)
C
C	FMS Parameter values:
	INTEGER MAXCPU
C
C	Local variables:
	INTEGER     IDTYPE, I, LX, LENX, LDISK, NV
	REAL*8      EI, ERROR
C
C	Common block to communicate with fill routines:
	COMMON /MYDATA/MOD, N, NRHS
	COMMON /SHARE/IWORK
	DATA LOWEQ/-1/
C
C (1)   Initialize FMS:
	CALL FMSINI
	CALL FMSGST ('MAXCPU', MAXCPU)
	CALL FMSIST ('MDATAU', 2)
	CALL FMSIGT ('MEMPTR', IMD_PTR)
	CALL FMSIGT ('MEMPTR', RMD_PTR)
	CALL FMSIGT ('MEMPTR', CMD_PTR)
	CALL FMSIST ('LOGTIM',    3)
	CALL FMSIST ('IPRF'  , 1026)
    1	CONTINUE
	WRITE (6,*) 'Enter the FMS module number (1 to 5)'
	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'
	READ  (5,*) MOD
	IF( (MOD.LT.1) .OR. (MOD.GT.5) ) GO TO 1
	N    = ASK_I('Enter the number of equations')
	NRHS = ASK_I('Enter the number of solution vectors')
	WRITE (6,*) 'You may now alter any FMS parameter.'
	WRITE (6,*) 'When you are finished, type the letters RETURN'
	CALL FMSSET
	IF(MOD .LE. 2) THEN
	   IDTYPE = 1
	ELSE
	   IDTYPE = 2
	END IF
	NUMEQ  = N
C
C (2)   Open FMS files:
	IF(MOD.EQ.1) CALL RSDI (LOWEQ, N, 'LUF', LUF)
	IF(MOD.EQ.2) CALL RNDI (LOWEQ, N, 'LUF', LUF)
	IF(MOD.EQ.3) CALL CHDI (LOWEQ, N, 'LUF', LUF)
	IF(MOD.EQ.4) CALL CSDI (LOWEQ, N, 'LUF', LUF)
	IF(MOD.EQ.5) CALL CNDI (LOWEQ, N, 'LUF', LUF)
	CALL FMSOV (N, IDTYPE, NRHS, 'LUX', LUX)
C
C (3)   Write data to FMS files:
C
C	Fill the RHS vectors in parallel:
	WRITE (6,2000)
	IWORK = 0
	DO ICPU = 1,(MAXCPU-1)
	   CALL FMSPAR (1, WRITEV, LUX)
	END DO
	IF(MAXCPU .GT. 1) CALL FMSRUN
	CALL WRITEV(LUX)
	IF(MAXCPU .GT. 1) CALL FMSYNC
C
C (4)   Perform matrix algebra:
	IF(MOD.EQ.1) CALL RSDF (LUF, LUF, LUX, LUX,  NUMRED)
	IF(MOD.EQ.2) CALL RNDF (LUF, LUF, LUX, LUX,  NUMRED)
	IF(MOD.EQ.3) CALL CHDF (LUF, LUF, LUX, LUX,  NUMRED)
	IF(MOD.EQ.4) CALL CSDF (LUF, LUF, LUX, LUX,  NUMRED)
	IF(MOD.EQ.5) CALL CNDF (LUF, LUF, LUX, LUX,  NUMRED)
C
	CALL FMSIGT ('LOGTIM', LOGTIM)
	IF(LOGTIM .GE. 3) CALL FMSCST ('SHOW', 'CALLS')
	IF(MOD.EQ.1) CALL RSDS (LUF, LUX, LUX, NRHS, ISKIP)
	IF(MOD.EQ.2) CALL RNDS (LUF, LUX, LUX, NRHS, ISKIP)
	IF(MOD.EQ.3) CALL CHDS (LUF, LUX, LUX, NRHS, ISKIP)
	IF(MOD.EQ.4) CALL CSDS (LUF, LUX, LUX, NRHS, ISKIP)
	IF(MOD.EQ.5) CALL CNDS (LUF, LUX, LUX, NRHS, ISKIP)
C
C (5)   Read data from FMS files:
C       Check the answer:
	ERROR = 0.0D0
	LDISK  = 1
	DO 60 NV = 1,NRHS
	   IF(IDTYPE .EQ. 1) THEN
	      CALL FMSRED (LUX(1), LDISK, RMD(LX), N)
	      DO 50 I = 1,N
	         EI = ABS(RMD(LX-1+I) - 1.0D0)
	         IF(EI .GT. ERROR) ERROR = EI
   50         CONTINUE
	   ELSE
	      CALL FMSRED (LUX(1), LDISK, CMD(LX), 2*N)
	      DO 51 I = 1,N
	         EI = ABS(CMD(LX-1+I) - 1.0D0)
	         IF(EI .GT. ERROR) ERROR = EI
   51         CONTINUE
	   END IF
	   LDISK  = LDISK  + LUX(4)
   60   CONTINUE
	WRITE(6,*) 'MAXIMUM ERROR =', ERROR
C
C (6)   Close FMS files:
	CALL FMSCM (LUF)
	CALL FMSCV (LUX)
	IF(ASK('Do you want another solution?')) GO TO 1
	CALL FMSEND
 2000	FORMAT (/
     1	' Writing the RHS vectors in parallel'/
     2	' ===================================')
	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(*), DIA, OFFDIA
	PARAMETER  (DIA   = 1.0D0)
	PARAMETER  (OFFDIA=-1.0D0)
	COMMON /MYDATA/MOD, N, NRHS
C
C       Populate the diagonal with test data:
	IF(IROW2 .EQ. JCOL2) THEN
C	   This is a diagonal block:
	   DO 10 I = IROW1,IROW2
	      D(I) = DIA
   10	   CONTINUE
	   IF(IROW1 .EQ. 1) D(1) = DFLOAT(N)
	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 20 I = IROW1,IROW2
	   J = LOWEQ(I)
	   IF( (J .GE. JCOL1) .AND.
     1	       (J .LE. JCOL2) .AND.
     2	       (J .LT.     I) ) A(LOCEQ(I) + IJSTEP*J) = OFFDIA
   20	CONTINUE
	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(*), OFFDIA
	PARAMETER  (OFFDIA=-1.0D0)
C
C	Populate profile of [AL] with test data:
C	The term A(I,J) is addressed as A(LOCEQ(J)+IJSTEP*I)
	DO 20 J = JCOL1,JCOL2
	   I = LOWEQ(J)
	   IF( (I .GE. IROW1) .AND.
     1	       (I .LE. IROW2) .AND.
     2	       (I .LT.     J) ) A(LOCEQ(J) + IJSTEP*I) = OFFDIA
   20	CONTINUE
	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(*), DIA
	COMPLEX*16 A(0:*), OFFDIA
	PARAMETER  (DIA   =  1.0D0        )
	PARAMETER  (OFFDIA=(-1.0D0,0.0D0))
	COMMON /MYDATA/MOD, N, NRHS
C
C       Populate the diagonal with test data:
	IF(IROW2 .EQ. JCOL2) THEN
C	   This is a diagonal block:
	   DO 10 I = IROW1,IROW2
	      D(I) = DIA
   10	   CONTINUE
	   IF(IROW1 .EQ. 1) D(1) = DFLOAT(N)
	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 20 I = IROW1,IROW2
	   J = LOWEQ(I)
	   IF( (J .GE. JCOL1) .AND.
     1	       (J .LE. JCOL2) .AND.
     2	       (J .LT.     I) ) A(LOCEQ(I) + IJSTEP*J) = OFFDIA
   20	CONTINUE
	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(*), DIA, OFFDIA
	PARAMETER  (DIA   =( 1.0D0, 1.0D0))
	PARAMETER  (OFFDIA=(-1.0D0,-1.0D0))
	COMMON /MYDATA/MOD, N, NRHS
C
C       Populate the diagonal with test data:
	IF(IROW2 .EQ. JCOL2) THEN
C	   This is a diagonal block:
	   DO 10 I = IROW1,IROW2
	      D(I) = DIA
   10	   CONTINUE
	   IF(IROW1 .EQ. 1) D(1) = DCMPLX(N,N)
	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 20 I = IROW1,IROW2
	   J = LOWEQ(I)
	   IF( (J .GE. JCOL1) .AND.
     1	       (J .LE. JCOL2) .AND.
     2	       (J .LT.     I) ) A(LOCEQ(I) + IJSTEP*J) = OFFDIA
   20	CONTINUE
	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(*), OFFDIA
	PARAMETER  (OFFDIA=(-1.0D0,-1.0D0))
C
C	Populate profile of [AL] with test data:
C	The term A(I,J) is addressed as A(LOCEQ(J)+IJSTEP*I)
	DO 20 J = JCOL1,JCOL2
	   I = LOWEQ(J)
	   IF( (I .GE. IROW1) .AND.
     1	       (I .LE. IROW2) .AND.
     2	       (I .LT.     J) ) A(LOCEQ(J) + IJSTEP*I) = OFFDIA
   20	CONTINUE
	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:*), DIA, OFFDIA
	PARAMETER  (DIA   = 1.0D0)
	PARAMETER  (OFFDIA=-1.0D0)
C
C       Fill first row and diagonal:
	LU1    = LUFLAG(1)
	LI1    = LOCI(1)
	DO 10 J = JEQN1,JEQN2
	   LIJ = LI1 + LOCJ(J,LU1)
	   A(LIJ) = OFFDIA
	   LD  = LOCI(J) + LOCJ(J,LUFLAG(J))
	   A(LD ) = DIA
   10	CONTINUE
C
C	Fill column 1 if this is the first slab:
	IF(JEQN1 .EQ. 1) THEN
	   DO 20 I = 1,NUMEQ
	      LIJ = LOCI(I) + LOCJ(1,LUFLAG(I))
	      A(LIJ) = OFFDIA
   20	   CONTINUE
	   LD  = LOCI(1) + LOCJ(1,LUFLAG(1))
	   A(LD) = DFLOAT(NUMEQ)
	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))
C
C       Fill first row and diagonal:
	LU1    = LUFLAG(1)
	LI1    = LOCI(1)
	DO 10 J = JEQN1,JEQN2
	   LIJ = LI1 + LOCJ(J,LU1)
	   A(LIJ) = OFFDIA
	   LD  = LOCI(J) + LOCJ(J,LUFLAG(J))
	   A(LD ) = DIA
   10	CONTINUE
C
C	Fill column 1 if this is the first slab:
	IF(JEQN1 .EQ. 1) THEN
	   DO 20 I = 1,NUMEQ
	      LIJ = LOCI(I) + LOCJ(1,LUFLAG(I))
	      A(LIJ) = OFFDIA
   20	   CONTINUE
	   LD  = LOCI(1) + LOCJ(1,LUFLAG(1))
	   A(LD) = DCMPLX(NUMEQ,NUMEQ)
	END IF
	RETURN
	END
C=======================================================================
	SUBROUTINE WRITEV (LUX)
C=======================================================================
C
C	DESCRIPTION:
C	   This subroutine runs in parallel. For each vector
C	   that it finds, it fills in the terms and calls FMSWRT to
C	   write the vector to the vector file LUX.
C
C	FORMAL PARAMETERS:
C	   (R ) LUX(25) = FMS vector file attributes
C-----------------------------------------------------------------------
C
C	Formal Parameters:
	INTEGER    LUX(25)
C
C	Local Variables:
	INTEGER    MYNODE
	INTEGER    LUPR
	INTEGER    NUMEQ
	INTEGER    NUMVEC
	INTEGER    NROWS, NCOLS
	INTEGER    IEQN1, IEQN2, JEQN1, JEQN2
	INTEGER    MY_WORK
	INTEGER    MY_TOTAL
	INTEGER    LX
	REAL*8     RMD(0:1)
	COMPLEX*16 CMD(0:1)
	POINTER    (RMD_PTR,RMD)
	POINTER    (CMD_PTR,CMD)
	COMMON /MYDATA/MOD, N, NRHS
	COMMON/SHARE/IWORK
	INTEGER    PARENT
	PARAMETER (PARENT = 0)
C
C	Get your CPU number (0 to MAXCPU-1)
	CALL FMSIGT ('MYNODE', MYNODE)
	CALL FMSIGT ('LUPR',   LUPR)
	NUMEQ  = LUX( 3)
	LENVEC = LUX( 4)
	IDTYPE = LUX( 5)
	NUMVEC = LUX( 6)
	NROWS  = NUMEQ
	NCOLS  = 1
	IEQN1  = 1
	IEQN2  = NUMEQ
	LDX    = LENVEC/IDTYPE
C
C	Allocate a work array for storing the data you generate:
C	For this example, we will fill [B] vector by vector.
	IF(IDTYPE .EQ. 1) THEN
	   CALL FMSIGT ('MEMPTR', RMD_PTR)
	   CALL FMSRMG (RMD, LX, LDX)
	   DO I=1,NUMEQ
	      RMD(LX+I-1) = 0.0D0
	   END DO
	   RMD(LX) = 1.0D0
	ELSE
	   CALL FMSIGT ('MEMPTR', CMD_PTR)
	   CALL FMSCMG (CMD, LX, LDX)
	   DO I=1,NUMEQ
	      CMD(LX+I-1) = (0.0D0,0.0D0)
	   END DO
	   IF(MOD .EQ. 3) THEN
	      CMD(LX) = (1.0D0,0.0D0)
	   ELSE
	      CMD(LX) = (1.0D0,1.0D0)
	   END IF
	END IF
C
C	Find the next vector to fill:
	MY_TOTAL = 0
   10	CONTINUE
	   CALL FMSONE
C	      IWORK   = IWORK + 1
C	      MY_WORK = IWORK
	      MY_WORK = INTINC(IWORK)
	   CALL FMSALL
	   IF(MY_WORK .GT. NUMVEC) GO TO 100
	   MY_TOTAL = MY_TOTAL + 1
	   LDISK  = 1 + LENVEC*(MY_WORK-1)
C
	   IF(IDTYPE .EQ. 1) THEN
	      CALL FMSWRT (LUX, LDISK, RMD(LX), LENVEC)
	   ELSE
	      CALL FMSPUT (LUX, LDISK, CMD(LX), LENVEC)
	   END IF
	GO TO 10
  100	CONTINUE
C	Your part of the work is done:
C
C	Return work array:
	IF(IDTYPE .EQ. 1) THEN
	   CALL FMSRMR (RMD, LX, LDX)
	ELSE
	   CALL FMSCMR (CMD, LX, LDX)
	END IF
C
C	Report your total work:
	CALL FMSONE
	   WRITE(LUPR,2000) MYNODE, MY_TOTAL
	CALL FMSALL
	RETURN
 2000	FORMAT (' Process',I3,' computed',I5,' vectors.')
 2001	FORMAT (' Process',I3,' is computing vector ',I5)
	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=======================================================================
        INTEGER FUNCTION INTINC (I)
C-----------------------------------------------------------------------
C	This function increments a volatile shared variable.  It is
C	placed in a subroutine to prevent some compilers from storing
C	the value in a register and not updating it.
	INTEGER I
	I = I + 1
	INTINC = I
	RETURN
	END

Go To Top Go Up Go Back Go Forward
Copyright © Multipath Corporation