This example illustrates how to fill a matrix in parallel
using subroutine FMSPUT.
C E X A M P L E 13
C
C Program name:
CHARACTER*10 MYNAME
PARAMETER (MYNAME='EXAMPLE_13')
C
C Subroutines called in parallel:
EXTERNAL WRITEM
EXTERNAL WRITEV
C
C FMS Module number (1=RS, 2=RN, 3=CH, 4=CS, 5=CN)
INTEGER MOD
C
C Number of equations:
INTEGER N
C
C Number of R.H.S. vectors:
INTEGER NRHS
C
C Size of block to fill matrix:
INTEGER NROWS, NCOLS
C
C Problem data:
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
C Work queue index (shared):
COMMON/SHARE/IWORK
INTEGER IWORK
C
C FMS matrix and vector file attributes:
INTEGER LUA(15), LUF(25), LUA0(25)
DATA LUA0(1)/0/
INTEGER LUX(25)
INTEGER LUS0(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 Input Data:
LOGICAL ASK
INTEGER ASK_I
C
C FMS Parameters:
INTEGER MDATAU
C
C Profile vector for a full matrix:
INTEGER LOWEQ(1)
DATA LOWEQ/-1/
C
C Local variables:
INTEGER MAXCPU
INTEGER IDTYPE
INTEGER ISTYPE
INTEGER I, LX, LDISK, NV
REAL*8 EI
LOGICAL ITEST
C
C Number of vectors to reduce during factoring:
PARAMETER (NUMRED = 0)
C
C Skip operations during solving (no):
PARAMETER (ISKIP = 0)
C
C Constants:
REAL*8 RALPHA(1)
DATA RALPHA/1.0D0/
COMPLEX*16 CALPHA(1)
DATA CALPHA/(1.0D0,0.0D0)/
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH(MYNAME)
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')
ISCALE = 1
DO I=1,15
IF(10**ISCALE .GT. N) GO TO 11
ISCALE = ISCALE + 1
END DO
11 CONTINUE
NRHS = ASK_I('Enter the number of solution vectors')
NROWS= ASK_I('Enter the number of rows per block (0=NEQBIO)')
NCOLS= ASK_I('Enter the number of columns per block (0=NEQBIO)')
IF(ASK('Do you want to check the matrix data?')) THEN
MDATAU = 2
ELSE
MDATAU = 0
END IF
CALL FMSIST ('MDATAU', MDATAU)
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)
CALL FMSIGT ('MAXCPU', MAXCPU)
IF(MOD .LE. 2) THEN
IDTYPE = 1
ELSE
IDTYPE = 2
END IF
ISTYPE = 1
IF(MOD .EQ. 2) ISTYPE = 2
IF(MOD .EQ. 3) ISTYPE = 3
IF(MOD .EQ. 5) ISTYPE = 2
C
C (2) Open FMS files:
CALL FMSIST ('IOKIDS', MAXCPU-1)
IF(MOD.EQ.1) CALL RSDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.2) CALL RNDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.3) CALL CHDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.4) CALL CSDI (LOWEQ, N, 'LUA', LUA)
IF(MOD.EQ.5) CALL CNDI (LOWEQ, N, 'LUA', LUA)
CALL FMSOV (N, IDTYPE, NRHS, 'LUX', LUX)
C
C (3) Write data to FMS files:
C
C Fill the RHS vectors in parallel:
C Initialize FMSPUT:
CALL FMSPUT (LUX, 0, 0, 0, 0, RMD, 1)
C Generate and write the vector data in parallel:
IWORK = 0
DO 30 ICPU=2,MAXCPU
CALL FMSPAR (1, WRITEV, LUX)
30 CONTINUE
IF(MAXCPU .GT. 1) CALL FMSRUN
CALL WRITEV (LUX)
IF(MAXCPU .GT. 1) CALL FMSYNC
C
C End FMSPUT:
CALL FMSPUT (LUX, 0, 0, N+1, 0, RMD, 1)
C
C If testing the matrix fill, first fill with one processor:
IF(MDATAU .NE. 0) THEN
IF(ISTYPE .EQ. 2) THEN
C This matrix is nonsymmetric:
CALL FMSPUT (LUA, 0, 0, 0, 0, RMD, 1)
ELSE
C This matrix is symmetric or Hermitian:
C Initialize FMSPUT to define the matrix terms in the upper
C triangle.
CALL FMSPUT (LUA, 0, 0, 0, 0, RMD, 1)
END IF
ISIGN = -1
IWORK = 0
CALL WRITEM (NROWS, NCOLS, LUA)
C End FMSPUT:
CALL FMSPUT (LUA, 0, 0, N+1, 0, RMD, 1)
END IF
C
C Fill the matrix in parallel:
C
C Initialize FMSPUT:
IF(ISTYPE .EQ. 2) THEN
C This matrix is nonsymmetric:
CALL FMSPUT (LUA, 0, 0, 0, 0, RMD, 1)
ELSE
C This matrix is symmetric or Hermitian:
C Initialize FMSPUT to define the matrix terms in the upper
C triangle.
CALL FMSPUT (LUA, 0, 0, 0, 0, RMD, 1)
END IF
ISIGN = 1
IWORK = 0
DO 35 ICPU=2,MAXCPU
CALL FMSPAR (3, WRITEM, NROWS, NCOLS, LUA)
35 CONTINUE
IF(MAXCPU .GT. 1) CALL FMSRUN
CALL WRITEM (NROWS, NCOLS, LUA)
IF(MAXCPU .GT. 1) CALL FMSYNC
C
C End FMSPUT:
CALL FMSPUT (LUA, 0, 0, N+1, 0, RMD, 1)
C
C (4) Perform matrix algebra:
C
C This example offers two testing options:
C MDATAU=0; Factor the matrix with valid test data;
C MDATAU=2; Check the matrix data and quit.
ERROR = 0.0D0
IF(MDATAU .EQ. 0) THEN
DO I=1,25
LUF(I) = LUA(I)
END DO
ELSE
LUF(1) = 0
END IF
IF(MOD.EQ.1) CALL RSDAF (LUA, RALPHA, 1, LUS0, 0, LUA0, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.2) CALL RNDAF (LUA, RALPHA, 1, LUS0, 0, LUA0, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.3) CALL CHDAF (LUA, RALPHA, 1, LUS0, 0, LUA0, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.4) CALL CSDAF (LUA, CALPHA, 1, LUS0, 0, LUA0, LUF, LUX, LUX, NUMRED)
IF(MOD.EQ.5) CALL CNDAF (LUA, CALPHA, 1, LUS0, 0, LUA0, LUF, LUX, LUX, NUMRED)
C
IF(MDATAU .EQ. 2) GO TO 70
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:
C
C Allocate a vector for reading {X}:
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMG (RMD, LX, N)
ELSE
CALL FMSCMG (CMD, LX, N)
END IF
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
70 CONTINUE
WRITE(6,*) 'MAXIMUM ERROR =', ERROR
C
C (6) Close FMS files:
CALL FMSCM (LUF)
CALL FMSCV (LUX)
IF(MDATAU .EQ. 0) THEN
IF(IDTYPE .EQ. 1) THEN
CALL FMSRMR (RMD, LX, N)
ELSE
CALL FMSCMR (CMD, LX, N)
END IF
END IF
IF(ASK('Do you want another solution?')) GO TO 1
CALL FMSPOP(MYNAME)
CALL FMSEND
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 FMSPUT 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:
CHARACTER*6 MYNAME
PARAMETER (MYNAME='WRITEV')
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/SHARE/IWORK
INTEGER PARENT
PARAMETER (PARENT = 0)
C
CALL FMSPSH(MYNAME)
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, NUMEQ)
ELSE
CALL FMSIGT ('MEMPTR', CMD_PTR)
CALL FMSCMG (CMD, LX, NUMEQ)
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)
C Write message under protection of lock:
IF(MY_WORK .LE. NUMVEC) WRITE(LUPR,2001) MYNODE, MY_WORK
CALL FMSALL
IF(MY_WORK .GT. NUMVEC) GO TO 100
MY_TOTAL = MY_TOTAL + 1
JEQN1 = MY_WORK
JEQN2 = JEQN1
C
IF(IDTYPE .EQ. 1) THEN
CALL RVFILL (RMD(LX), IEQN1, IEQN2, JEQN1, JEQN2, NUMEQ)
CALL FMSPUT (LUX, NROWS, NCOLS, IEQN1, JEQN1, RMD(LX),
1 LDX)
ELSE
CALL CVFILL (CMD(LX), IEQN1, IEQN2, JEQN1, JEQN2, NUMEQ)
CALL FMSPUT (LUX, NROWS, NCOLS, IEQN1, JEQN1, CMD(LX),
1 LDX)
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, N)
ELSE
CALL FMSCMR (CMD, LX, N)
END IF
C
C Report your total work:
CALL FMSONE
WRITE(LUPR,2000) MYNODE, MY_TOTAL
CALL FMSALL
CALL FMSPOP(MYNAME)
RETURN
2000 FORMAT (' Process',I3,' computed',I5,' vectors.')
2001 FORMAT (' Process',I3,' is computing vector ',I5)
END
C=======================================================================
SUBROUTINE WRITEM (NROWS, NCOLS, LUA)
C=======================================================================
C
C DESCRIPTION:
C This subroutine runs in parallel. For each block
C A(NROWS, NCOLS) that it finds, it fills in the terms and
C calls FMSPUT to write the block to the matrix file LUA.
C
C FORMAL PARAMETERS:
C (R ) NROWS = Number of rows per block
C (R ) NCOLS = Number of columns per block
C (R ) LUA(25) = FMS matrix file attributes
C-----------------------------------------------------------------------
C
C Formal Parameters:
INTEGER NROWS
INTEGER NCOLS
INTEGER LUA(25)
C
C Local Variables:
CHARACTER*6 MYNAME
PARAMETER (MYNAME='WRITEM')
INTEGER NEQROW, NEQCOL
INTEGER NBKROW, NBKCOL
INTEGER IEQN1, IEQN2, JEQN1, JEQN2
INTEGER IBLOCK, JBLOCK
INTEGER NUMBLK
INTEGER MYNODE
INTEGER MY_WORK
INTEGER MY_TOTAL
INTEGER LUPR
INTEGER NUMEQ
INTEGER NEQBIO
INTEGER NUMSEG
INTEGER IDTYPE
INTEGER LA
REAL*8 RMD(0:1)
COMPLEX*16 CMD(0:1)
POINTER (RMD_PTR,RMD)
POINTER (CMD_PTR,CMD)
COMMON/SHARE/IWORK
INTEGER PARENT
PARAMETER (PARENT = 0)
C
CALL FMSPSH(MYNAME)
C Get your CPU number (0 to MAXCPU-1)
CALL FMSIGT ('MYNODE', MYNODE)
CALL FMSIGT ('LUPR', LUPR)
NUMEQ = LUA( 8)
NUMSEG = LUA( 6)
IDTYPE = LUA(11)
ISTYPE = LUA(12)
MFMAT = LUA(13)
NEQBIO = LUA(21)
NEQROW = NROWS
NEQCOL = NCOLS
IF(NEQROW .EQ. 0) THEN
IF(MFMAT .EQ. 1) THEN
NEQROW = NUMEQ/5
ELSE
NEQROW = NEQBIO
END IF
END IF
IF(NEQCOL .EQ. 0) THEN
IF(MFMAT .EQ. 1) THEN
NEQCOL = NUMEQ/5
ELSE
NEQCOL = NEQBIO
END IF
END IF
NBKROW = (NUMEQ+NEQROW-1)/NEQROW
NBKCOL = (NUMEQ+NEQCOL-1)/NEQCOL
NUMBLK = NBKROW*NBKCOL
LDA = NEQROW
C
C Allocate a work array for storing the data you generate:
CALL FMSIGT ('MDATAU', MDATAU)
IF(IDTYPE .EQ. 1) THEN
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSRMG (RMD, LA, NEQROW*NEQCOL)
ELSE
CALL FMSIGT ('MEMPTR', CMD_PTR)
CALL FMSCMG (CMD, LA, NEQROW*NEQCOL)
END IF
C
C Fill the matrix:
MY_TOTAL = 0
10 CONTINUE
CALL FMSONE
C IWORK = IWORK + 1
C MY_WORK = IWORK
MY_WORK = INTINC(IWORK)
C Write message under protection of lock:
JBLOCK = (MY_WORK + NBKROW - 1)/NBKROW
IBLOCK = MY_WORK - NBKROW*(JBLOCK-1)
IF(MY_WORK .LE. NUMBLK)
1 WRITE(LUPR,2001) MYNODE, IBLOCK, JBLOCK
CALL FMSALL
IF(MY_WORK .GT. NUMBLK) GO TO 100
MY_TOTAL = MY_TOTAL + 1
C
IEQN1 = 1 + NEQROW*(IBLOCK-1)
JEQN1 = 1 + NEQCOL*(JBLOCK-1)
IEQN2 = IEQN1 + NEQROW - 1
IEQN2 = MIN0(IEQN2,NUMEQ)
JEQN2 = JEQN1 + NEQCOL - 1
JEQN2 = MIN0(JEQN2,NUMEQ)
NR = IEQN2 - IEQN1 + 1
NC = JEQN2 - JEQN1 + 1
IF(ISTYPE .NE. 2) THEN
C This problem is symmetric.
C Skip blocks entirely in [L]:
IF(IEQN1 .GT. JEQN2) GO TO 10
END IF
IF(IDTYPE .EQ. 1) THEN
CALL RMFILL (RMD(LA), NR, NC, IEQN1, JEQN1, NUMEQ)
IF(MDATAU .GT. 0)
1 CALL RCHECK (RMD(LA), NR, NC, IEQN1, JEQN1)
CALL FMSPUT (LUA, NR, NC, IEQN1, JEQN1, RMD(LA), NR)
ELSE
CALL CMFILL (CMD(LA), NR, NC, IEQN1, JEQN1, NUMEQ)
IF(MDATAU .GT. 0)
1 CALL CCHECK (CMD(LA), NR, NC, IEQN1, JEQN1)
CALL FMSPUT (LUA, NR, NC, IEQN1, JEQN1, CMD(LA), NR)
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, NEQROW*NEQCOL)
ELSE
CALL FMSCMR (CMD, LX, NEQROW*NEQCOL)
END IF
C
C Report your total work:
CALL FMSONE
WRITE(LUPR,2000) MYNODE, MY_TOTAL
CALL FMSALL
CALL FMSPOP(MYNAME)
RETURN
2000 FORMAT (' Process',I3,' computed',I5,' blocks.')
2001 FORMAT (' Process',I3,' is computing block (',I3,',',I3,')')
END
C=======================================================================
SUBROUTINE RVFILL (X, IROW1, IROW2, JCOL1, JCOL2, NUMEQ)
C=======================================================================
C Fill a test vector with REAL*8 data:
REAL*8 X(IROW1:IROW2,JCOL1:JCOL2)
CHARACTER*6 MYNAME
PARAMETER (MYNAME='RVFILL')
INTEGER IROW1
INTEGER IROW2
INTEGER JCOL1
INTEGER JCOL2
REAL*8 ZERO, X1
DATA ZERO/0.0D0/
DATA X1 /1.0D0/
C
CALL FMSPSH(MYNAME)
C Zero out [X]
DO IROW=IROW1,IROW2
DO JCOL=JCOL1,JCOL2
X(IROW,JCOL) = ZERO
END DO
END DO
C
C Set first row to X1:
IF(IROW1 .EQ. 1) THEN
DO JCOL=JCOL1,JCOL2
X(1,JCOL) = X1
END DO
END IF
CALL FMSPOP(MYNAME)
RETURN
END
C=======================================================================
SUBROUTINE CVFILL (X, IROW1, IROW2, JCOL1, JCOL2, NUMEQ)
C=======================================================================
C Fill a test vector with COMPLEX*16 data:
CHARACTER*6 MYNAME
PARAMETER (MYNAME='CVFILL')
COMPLEX*16 X(IROW1:IROW2,JCOL1:JCOL2)
INTEGER IROW1
INTEGER IROW2
INTEGER JCOL1
INTEGER JCOL2
COMPLEX*16 ZERO, X1
DATA ZERO/(0.0D0, 0.0D0)/
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
CALL FMSPSH(MYNAME)
IF(MOD .EQ. 3) THEN
X1 = (1.0D0,0.0D0)
ELSE
X1 = (1.0D0,1.0D0)
END IF
C
C Zero out [X]
DO IROW=IROW1,IROW2
DO JCOL=JCOL1,JCOL2
X(IROW,JCOL) = ZERO
END DO
END DO
C
C Set first row to X1:
IF(IROW1 .EQ. 1) THEN
DO JCOL=JCOL1,JCOL2
X(1,JCOL) = X1
END DO
END IF
CALL FMSPOP(MYNAME)
RETURN
END
C=======================================================================
SUBROUTINE RMFILL (A, NR, NC, IROW1, JCOL1, NUMEQ)
C=======================================================================
C Fill a matrix block with REAL*8 data:
REAL*8 A(NR,NC)
INTEGER NR
INTEGER NC
INTEGER IROW1
INTEGER JCOL1
INTEGER I, J, IROW, JCOL
REAL*8 OFFDIA, ZERO, DIA, AI, AJ, SCALE
CHARACTER*6 MYNAME
PARAMETER (MYNAME='RMFILL')
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
DATA OFFDIA/-1.0D0/
DATA ZERO / 0.0D0/
DATA DIA / 1.0D0/
C
CALL FMSPSH(MYNAME)
C Get MDATAU:
CALL FMSIGT ('MDATAU', MDATAU)
C
IF(MDATAU .EQ. 0) THEN
C Fill the block with test data:
C Zero out [A]
DO I=1,NR
DO J=1,NC
A(I,J) = ZERO
END DO
END DO
C
C Set first column to OFFDIA:
IF(JCOL1 .EQ. 1) THEN
DO I=1,MIN0(NR,NUMEQ-IROW1+1)
A(I,1) = OFFDIA
END DO
END IF
C
C Set first row to OFFDIA:
IF(IROW1 .EQ. 1) THEN
DO J=1,MIN0(NC,NUMEQ-JCOL1+1)
A(1,J) = OFFDIA
END DO
END IF
C
C Set diagonal to DIA:
IDIA1 = MAX0(IROW1,JCOL1)
IDIA2 = MIN0(IROW1+NR-1,JCOL1+NC-1,NUMEQ)
IF(IDIA1 .LE. IDIA2) THEN
DO IDIA = IDIA1,IDIA2
A(IDIA-IROW1+1,IDIA-JCOL1+1) = DIA
END DO
END IF
C
C Set A(1,1) to NUMEQ
IF( (IROW1 .EQ. 1) .AND.
1 (JCOL1 .EQ. 1) ) THEN
A(1,1) = DFLOAT(NUMEQ)
END IF
ELSE
C
C Fill the block with (row,column) numbers:
SCALE = 10**ISCALE
SCALE = ISIGN * SCALE
IF(MOD .EQ. 1) THEN
C Matrix is real symmetric:
DO J=1,NC
JCOL = J + JCOL1 - 1
AJ = DFLOAT(JCOL)
DO I=1,NR
IROW = I + IROW1 - 1
AI = DFLOAT(IROW)
IF(JCOL .LE. IROW) THEN
A(I,J) = SCALE*AI + AJ
ELSE
A(I,J) = SCALE*AJ + AI
END IF
END DO
END DO
ELSE
C Matrix is real nonsymmetric:
DO J=1,NC
JCOL = J + JCOL1 - 1
AJ = DFLOAT(JCOL)
DO I=1,NR
IROW = I + IROW1 - 1
AI = SCALE*DFLOAT(IROW)
A(I,J) = AI + AJ
END DO
END DO
END IF
END IF
CALL FMSPOP(MYNAME)
RETURN
END
C=======================================================================
SUBROUTINE CMFILL (A, NR, NC, IROW1, JCOL1, NUMEQ)
C=======================================================================
C Fill a matrix block with COMPLEX*16 data:
COMPLEX*16 A(NR,NC)
INTEGER NR
INTEGER NC
INTEGER IROW1
INTEGER JCOL1
INTEGER NUMEQ
INTEGER I, J, IROW, JCOL
REAL*8 RN, SCALE, A_Re, A_Im
COMPLEX*16 ZERO, DIA, OFFDIA, A11
DATA ZERO/( 0.0D0, 0.0D0)/
CHARACTER*6 MYNAME
PARAMETER (MYNAME='CMFILL')
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
CALL FMSPSH(MYNAME)
C Get MDATAU:
CALL FMSIGT ('MDATAU', MDATAU)
C
IF(MDATAU .EQ. 0) THEN
C Fill the block with test data:
DIA = ( 1.0D0, 1.0D0)
OFFDIA = (-1.0D0,-1.0D0)
RN = DFLOAT(NUMEQ)
IF(MOD .EQ. 3) THEN
C Hermitian test data:
DIA = ( 1.0D0, 0.0D0)
OFFDIA = (-1.0D0, 0.0D0)
A11 = DCMPLX(RN,0.0D0)
ELSE
C Complex symmetric or nonsymmetric test data:
DIA = ( 1.0D0, 1.0D0)
OFFDIA = (-1.0D0,-1.0D0)
A11 = DCMPLX(RN,RN)
END IF
C
C Zero out [A]
DO I=1,NR
DO J=1,NC
A(I,J) = ZERO
END DO
END DO
C
C Set first column to (-1,-1):
IF(JCOL1 .EQ. 1) THEN
DO I=1,MIN0(NR,NUMEQ-IROW1+1)
A(I,1) = OFFDIA
END DO
END IF
C
C Set first row to (-1,-1):
IF(IROW1 .EQ. 1) THEN
DO J=1,MIN0(NC,NUMEQ-JCOL1+1)
A(1,J) = OFFDIA
END DO
END IF
C
C Set diagonal to (1,1):
IDIA1 = MAX0(IROW1,JCOL1)
IDIA2 = MIN0(IROW1+NR-1,JCOL1+NC-1,NUMEQ)
IF(IDIA1 .LE. IDIA2) THEN
DO IDIA = IDIA1,IDIA2
A(IDIA-IROW1+1,IDIA-JCOL1+1) = DIA
END DO
END IF
C
C Set A(1,1) to (NUMEQ,NUMEQ)
IF( (IROW1 .EQ. 1) .AND.
1 (JCOL1 .EQ. 1) ) THEN
A(1,1) = A11
END IF
ELSE
C
C Fill the block with (row,column) numbers:
IF(MOD .EQ. 3) THEN
C Hermitian test data; Real diagonals.
SCALE = 10**ISCALE
DO J=1,NC
JCOL = J + JCOL1 - 1
A_Im = DFLOAT(ISIGN*JCOL)
DO I=1,NR
IROW = I + IROW1 - 1
A_Re = DFLOAT(ISIGN*IROW)
IF(JCOL .LT. IROW) THEN
C This term is in the lower triangle:
A(I,J) = DCMPLX(A_Re,A_Im)
ELSE IF(JCOL .EQ. IROW) THEN
C This term is on the diagonal:
RN = SCALE * A_Re + A_Im
A(I,J) = DCMPLX(RN,0.0D0)
ELSE
C This term is in the upper triangle:
A(I,J) = DCMPLX(A_Im,-A_Re)
END IF
END DO
END DO
ELSE IF(MOD .EQ. 4) THEN
C Complex Symmetric:
DO J=1,NC
JCOL = J + JCOL1 - 1
A_Im = DFLOAT(ISIGN*JCOL)
DO I=1,NR
IROW = I + IROW1 - 1
A_Re = DFLOAT(ISIGN*IROW)
IF(JCOL .LE. IROW) THEN
A(I,J) = DCMPLX(A_Re,A_Im)
ELSE
A(I,J) = DCMPLX(A_Im,A_Re)
END IF
END DO
END DO
ELSE
C Complex Nonsymmetric:
DO J=1,NC
JCOL = J + JCOL1 - 1
A_Im = DFLOAT(ISIGN*JCOL)
DO I=1,NR
IROW = I + IROW1 - 1
A_Re = DFLOAT(ISIGN*IROW)
A(I,J) = DCMPLX(A_Re,A_Im)
END DO
END DO
END IF
END IF
CALL FMSPOP(MYNAME)
RETURN
END
C=======================================================================
SUBROUTINE RCHECK (A, NR, NC, IROW1, JCOL1)
C=======================================================================
C Check the matrix block:
REAL*8 A(NR,NC)
INTEGER NR
INTEGER NC
INTEGER IROW1
INTEGER JCOL1
INTEGER I, J, IROW, JCOL
REAL*8 SCALE, ATEST
CHARACTER*6 MYNAME
PARAMETER (MYNAME='RCHECK')
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
CALL FMSPSH(MYNAME)
SCALE = 10**ISCALE
SCALE = ISIGN * SCALE
DO I=1,NR
IROW = I + IROW1 - 1
AI = DFLOAT(IROW)
DO J=1,NC
JCOL = J + JCOL1 - 1
AJ = DFLOAT(JCOL)
IF(JCOL .LE. IROW) THEN
ATEST = SCALE*AI + AJ
ELSE
IF(MOD .EQ. 1) THEN
ATEST = SCALE*AJ + AI
ELSE
ATEST = SCALE*AI + AJ
END IF
END IF
IF(A(I,J) .NE. ATEST) CALL RERROR (IROW, JCOL, A(I,J), ATEST)
END DO
END DO
CALL FMSPOP(MYNAME)
RETURN
END
C=======================================================================
SUBROUTINE CCHECK (A, NR, NC, IROW1, JCOL1)
C=======================================================================
C Fill a matrix block with COMPLEX*16 data:
COMPLEX*16 A(NR,NC)
INTEGER NR
INTEGER NC
INTEGER IROW1
INTEGER JCOL1
INTEGER I, J, IROW, JCOL
REAL*8 RN, SCALE, A_Re, A_Im
COMPLEX*16 ATEST
CHARACTER*6 MYNAME
PARAMETER (MYNAME='CCHECK')
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
CALL FMSPSH(MYNAME)
SCALE = 10**ISCALE
DO I=1,NR
IROW = I + IROW1 - 1
A_Re = DFLOAT(ISIGN*IROW)
DO J=1,NC
JCOL = J + JCOL1 - 1
A_Im = DFLOAT(ISIGN*JCOL)
IF(JCOL .LT. IROW) THEN
ATEST = DCMPLX(A_Re,A_Im)
ELSE IF(JCOL .EQ. IROW) THEN
IF(MOD .EQ. 3) THEN
RN = SCALE * A_Re + A_Im
ATEST = DCMPLX(RN,0.0D0)
ELSE
ATEST = DCMPLX(A_Re,A_Im)
END IF
ELSE
IF(MOD .EQ. 3) THEN
ATEST = DCMPLX(A_Im,-A_Re)
ELSE IF(MOD .EQ. 4) THEN
ATEST = DCMPLX(A_Im,A_Re)
ELSE
ATEST = DCMPLX(A_Re,A_Im)
END IF
END IF
IF(A(I,J) .NE. ATEST) CALL CERROR (IROW,JCOL,A(I,J),ATEST)
END DO
END DO
CALL FMSPOP(MYNAME)
RETURN
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
C=======================================================================
SUBROUTINE RERROR (IROW, JCOL, AIJ, AVALUE)
C-----------------------------------------------------------------------
C This subroutine prints a real matrix term that is incorrect.
INTEGER IROW, JCOL
REAL*8 AIJ, AVALUE
CHARACTER*55 FMT
DATA FMT( 1:36)/'(''CPU___:A('',I_,'','',I_,'')='',F__.0,'';'/
DATA FMT(37:55)/' Should be '',F__.0)'/
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
INTEGER MYNODE
INTEGER NERRORS
DATA NERRORS/0/
CALL FMSIGT ('MYNODE', MYNODE)
ISCALE2 = MAX0(2*ISCALE,10)
WRITE(FMT( 6: 8), 2003) MYNODE
WRITE(FMT(15:15), 2001) ISCALE
WRITE(FMT(22:22), 2001) ISCALE
WRITE(FMT(30:31), 2002) ISCALE2
WRITE(FMT(51:52), 2002) ISCALE2
IF(ERROR .EQ. 0.0D0) ERROR = AVALUE
WRITE(6, FMT) IROW, JCOL, AIJ, AVALUE
NERRORS = NERRORS + 1
IF(NERRORS .GT. 100) STOP 'Fatal error checking Matrix Data'
2001 FORMAT (I1)
2002 FORMAT (I2)
2003 FORMAT (I3.3)
RETURN
END
C=======================================================================
SUBROUTINE CERROR (IROW, JCOL, AIJ, AVALUE)
C-----------------------------------------------------------------------
C This subroutine prints a real matrix term that is incorrect.
INTEGER IROW, JCOL
COMPLEX*16 AIJ, AVALUE
CHARACTER*82 FMT
DATA FMT( 1:34)/'(''CPU___:A('',I_,'','',I_,'')=('',F__.0'/
DATA FMT(35:48)/','','',F__.0,'');'/
DATA FMT(49:82)/'_Should_be_('',F__.0,'','',F__.0,'')'')'/
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
INTEGER MYNODE
INTEGER NERRORS
DATA NERRORS/0/
CALL FMSIGT ('MYNODE', MYNODE)
ISCALE1 = ISCALE + 1
WRITE(FMT( 6: 8), 2003) MYNODE
WRITE(FMT(15:15), 2001) ISCALE
WRITE(FMT(22:22), 2001) ISCALE
WRITE(FMT(31:32), 2002) ISCALE1
WRITE(FMT(41:42), 2002) ISCALE1
WRITE(FMT(64:65), 2002) ISCALE1
WRITE(FMT(74:75), 2002) ISCALE1
IF(ERROR .EQ. 0.0D0) ERROR = AVALUE
WRITE(6, FMT) IROW, JCOL, AIJ, AVALUE
NERRORS = NERRORS + 1
IF(NERRORS .GT. 100) STOP 'Fatal error checking Matrix Data'
2001 FORMAT (I1)
2002 FORMAT (I2)
2003 FORMAT (I3.3)
RETURN
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(*), AIJ, ATEST, SCALE
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
C Test the matrix data:
SCALE = 10**ISCALE
DO IROW = IROW1,IROW2
JSTART = MAX0(LOWEQ(IROW), JCOL1)
JEND = MIN0(IROW, JCOL2)
LOCI = LOCEQ(IROW)
DO JCOL = JSTART, JEND
IF(JCOL .EQ. IROW) THEN
C This term is on the diagonal:
AIJ = D(IROW)
ELSE
C This term is on or within the profile:
AIJ = A(LOCI + IJSTEP*JCOL)
END IF
ATEST = SCALE * DFLOAT(IROW) + DFLOAT(JCOL)
IF(AIJ .NE. ATEST) CALL RERROR (IROW, JCOL, AIJ, ATEST)
END DO
END DO
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(*), AIJ, ATEST, SCALE
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
C Test the matrix data:
SCALE = 10**ISCALE
DO JCOL = JCOL1,JCOL2
ISTART = MAX0(LOWEQ(JCOL), IROW1)
IEND = MIN0(JCOL-1, IROW2)
LOCJ = LOCEQ(JCOL)
DO IROW = ISTART, IEND
ATEST = SCALE * DFLOAT(IROW) + DFLOAT(JCOL)
AIJ = A(LOCJ + IJSTEP*IROW)
IF(AIJ .NE. ATEST) CALL RERROR (IROW, JCOL, AIJ, ATEST)
END DO
END DO
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(*), DII, DTEST, SCALE
COMPLEX*16 A(0:*), AIJ, ATEST
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
C Test the matrix data:
SCALE = 10**ISCALE
DO IROW = IROW1,IROW2
JSTART = MAX0(LOWEQ(IROW), JCOL1)
JEND = MIN0(IROW, JCOL2)
LOCI = LOCEQ(IROW)
DO JCOL = JSTART, JEND
IF(JCOL .EQ. IROW) THEN
C This term is on the diagonal:
DII = D(IROW)
DTEST = SCALE * DFLOAT(IROW) + DFLOAT(JCOL)
IF(DII .NE. DTEST) CALL RERROR (IROW, JCOL, DII, DTEST)
ELSE
C This term is on or within the profile:
AIJ = A(LOCI + IJSTEP*JCOL)
IF(JCOL .LT. IROW) THEN
C This term is in the lower triangle:
ATEST = DCMPLX(IROW,JCOL)
ELSE
C This term is in the upper triangle:
c Use complex conjugate:
ATEST = DCMPLX(JCOL,-IROW)
END IF
IF(AIJ .NE. ATEST) CALL CERROR (IROW, JCOL, AIJ, ATEST)
END IF
END DO
END DO
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(*), AIJ, ATEST
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
C Test the matrix data:
DO IROW = IROW1,IROW2
JSTART = MAX0(LOWEQ(IROW), JCOL1)
JEND = MIN0(IROW, JCOL2)
LOCI = LOCEQ(IROW)
DO JCOL = JSTART, JEND
IF(JCOL .EQ. IROW) THEN
C This term is on the diagonal:
AIJ = D(IROW)
ELSE
C This term is on or within the profile:
AIJ = A(LOCI + IJSTEP*JCOL)
END IF
ATEST = DCMPLX(IROW,JCOL)
IF(AIJ .NE. ATEST) CALL CERROR (IROW, JCOL, AIJ, ATEST)
END DO
END DO
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(*), AIJ, ATEST
C
C Test the matrix data:
DO JCOL = JCOL1,JCOL2
ISTART = MAX0(LOWEQ(JCOL), IROW1)
IEND = MIN0(JCOL-1, IROW2)
LOCJ = LOCEQ(JCOL)
DO IROW = ISTART, IEND
AIJ = A(LOCJ + IJSTEP*IROW)
ATEST = DCMPLX(IROW,JCOL)
IF(AIJ .NE. ATEST) CALL CERROR (IROW, JCOL, AIJ, ATEST)
END DO
END DO
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:*), AIJ, ATEST, SCALE
COMMON /MYDATA/ERROR, MOD, N, NRHS, NROWS, NCOLS, ISCALE, ISIGN
REAL*8 ERROR
C
C Test the matrix data:
SCALE = 10**ISCALE
DO JCOL = JEQN1,JEQN2
DO IROW = 1,NUMEQ
AIJ = A(LOCI(IROW) + LOCJ(1,LUFLAG(IROW)))
ATEST = SCALE * DFLOAT(IROW) + DFLOAT(JCOL)
IF(AIJ .NE. ATEST) CALL RERROR (IROW, JCOL, AIJ, ATEST)
END DO
END DO
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:*), AIJ, ATEST
C
C Test the matrix data:
DO JCOL = JEQN1,JEQN2
DO IROW = 1,NUMEQ
AIJ = A(LOCI(IROW) + LOCJ(1,LUFLAG(IROW)))
ATEST = DCMPLX(IROW,JCOL)
IF(AIJ .NE. ATEST) CALL CERROR (IROW, JCOL, AIJ, ATEST)
END DO
END DO
RETURN
END