C The FMS Out-of-Core matrix multiply routine FMSMM forms the
C product of two matrixes [A][B] and stores the result in matrix
C [C]. Under some matrix partitioning conditions, the size of
C the product [AB] = [A][B] will be different than the matrix [C].
C Normally, FMS will store the product [AB] starting at the upper
C left corner of [C], C(1,1). However, it may be necessary to
C shift the starting position of the product [AB] in [C] to some
C other position. The FMS Parameters MMROW and MMCOL specify
C the number of rows and columns to shift the product [AB] in
C [C], storing [AB] in [C] starting at C(1+MMROW,1+MMCOL). The
C default values of MMROW and MMCOL are 0.
C
C Another shifting requirement occurs when the matrix [A] has
C fewer columns than the matrix [B] has rows. In this case, FMS
C provides the shifting parameter MMKA, which has the effect of
C adding MMKA columns of 0's to matrix [A]. The matrix product
C then becomes
C
C [AB] = [0|A][B]
C
C The following example illustrates the use of FMSMM and these
C FMS Parameters.
C
C The system of linear equations is partitioned as follows:
C
C -- +- + + -+ +- -+ +- -+
C N1 | A11 | A12 | A13 | | X1 | | B1 |
C -- +-----+-----+-----+ +----+ +----+
C N2 | A21 | | | | = | |
C -- +-----+ A22 + + X2 + + B2 +
C N3 | A31 | | | | | |
C -- +- +- -+ +- -+ +- -+
C
C |< N1>|< N2>|< N3>| <NRHS> <NRHS>
C |<---N22--->|
C
C Because the matrices [A31] and [A13] do not lie on a boundary
C of [A22], it is necessary to use the matrix multiply shifting
C parameters during the partitioned solution of this system.
C The following table lists the products computed by FMSMM and
C the required values of the shifting parameters.
C
C Matrix Product Stored In MMROW MMCOL MMKA
C ============== ========= ===== ===== ====
C [A21][A12] [A22] 0 0 0
C [A21][A13] [A22] 0 N2 0
C [A21][X1] [X2] 0 0 0
C [A31][A12] [A22] N2 0 0
C [A31][A13] [A22] N2 N2 0
C [A31][X1] [X2] N2 0 0
C [0|A13][X2] [X1] 0 0 N2
C [A12][X2] [X1] 0 0 0
C
C Normally FMS will overlay matrix factors and solution vectors
C on the input matrix data. However, for this example, we want
C to save the original matrices to check the result. Therefore
C separate files will be used for the original matrix, factors
C and the solution vectors.
C
C The following FMS files will be used:
C
C File File
C Name Type CONTENTS
C ====== ====== ========
C LU_A11 Matrix [A11]
C LU_F11 Matrix [F11], factors of [A11]
C
C LU_A12 Vector [A12]
C LU_F12 Vector [F12], solution of [F11][F12]=[A12]
C
C LU_A13 Vector [A13]
C LU_F13 Vector [F13], solution of [F11][F13]=[A13]
C
C LU_A21 Vector [A21]T
C
C LU_A31 Vector [A31]T
C
C LU_A22 Matrix [A22]
C LU_F22 Matrix [F22], factors of [A22]-[A21|A31][F12|F13]
C
C LU_B1 Vector {B1}
C LU_X1 Vector {X1}
C
C LU_B2 Vector {B2}
C LU_X2 Vector {X2}
C
C For symmetric problems, the files LU_A12 and LU_A13 are not
C required, because they are the same as LU_A21 and LU_A31.
C
C INPUT:
C ======
C The following input is used:
C FMS Module (1=RS,2=RN,3=CH,4=CS,5=CN)
C Number of equations in [A11], N1
C Number of equations in [A21], N2
C Number of equations in [A31], N2
C Number of solution vectors, NRHS
C Matrix test data (1=known, 2=random), MTEST
C Any FMS Parameter
C
C
C For known matrix test data (MTEST=1), the following pattern
C is used:
C
C [A] {X} = {B}
C +- -+ + + + +
C | N -1 -1 -1 -1| | 1 | | 1 |
C |-1 1 0 0 0| | 1 | | 0 |
C |-1 0 1 0 0| | 1 | = | 0 |
C |-1 0 0 1 0| | 1 | | 0 |
C |-1 0 0 0 1| | 1 | | 0 |
C +- -+ + + + +
C
C where N=N1+N2+N3
C
C Program name:
CHARACTER*10 MYNAME
PARAMETER (MYNAME='EXAMPLE_15')
C
C FMS matrix and vector file attributes:
INTEGER LU_A11(25), LU_F11(25)
INTEGER LU_A12(25), LU_F12(25)
INTEGER LU_A13(25), LU_F13(25)
INTEGER LU_A21(25)
INTEGER LU_A31(25)
INTEGER LU_A22(25), LU_F22(25)
INTEGER LU_B1 (25), LU_X1 (25)
INTEGER LU_B2 (25), LU_X2 (25)
INTEGER LUA0 (25)
DATA LUA0(1)/0/
INTEGER LUX0 (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 LOCAL VARIABLES:
C FMS module number (1 to 5)
INTEGER MOD
C Number of equations in [A11]:
INTEGER N1
C Number of equations in [A21]:
INTEGER N2
C Number of equations in [A31]:
INTEGER N3
C Number of equations in [A22]: (=N2+N3)
C INTEGER N22
C Number of solution vectors:
INTEGER NRHS
C FMS Data Type (1 or 2)
INTEGER IDTYPE
C FMS Matrix Symmetry (1, 2 or 3):
INTEGER ISTYPE
C Test matrix data type (1=fixed, 2=random):
INTEGER MTEST
C Error values:
REAL*8 EI
REAL*8 ERROR
C Size of scratch work vector:
INTEGER N_TEMP
C RMD Pointer to scratch vector:
INTEGER L_TEMP
C Transpose of matrix A:
CHARACTER*1 TRANSA
C
C CONSTANTS:
REAL*8 R_ZERO, R_ONE(1), R_BIG
DATA R_ZERO/0.0D0/
DATA R_ONE /1.0D0/
DATA R_BIG /1.0D30/
COMPLEX*16 C_ONE(1)
DATA C_ONE /(1.0D0,0.0D0)/
C
C Profile vector for a full matrix:
INTEGER LOWEQ(1)
DATA LOWEQ/-1/
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH(MYNAME)
CALL FMSIGT ('LUPR' , LUPR)
C
C Use block format matrix:
CALL FMSIST ('MFMAT', 2)
C
C Read in problem size parameters:
10 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 10
END IF
N1 = ASK_I('Enter the number of equations in A11, N1')
N2 = ASK_I('Enter the number of columns in A12, N2')
N3 = ASK_I('Enter the number of columns in A13, N3')
NRHS = ASK_I('Enter the number of solution vectors, NRHS')
MTEST= ASK_I('Enter the test data (1=fixed,2=random)')
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)
N22 = N2 + N3
NUMEQ = N1 + N22
IF(MOD .LE. 2) THEN
C Real*8 data:
IDTYPE = 1
ELSE
C Complex*16 data:
IDTYPE = 2
END IF
LENDIA = IDTYPE*NUMEQ
TRANSA = 'T'
IF ((MOD .EQ. 1) .OR. (MOD .EQ. 4)) THEN
C Symmetric matrix:
ISTYPE = 1
ELSE IF((MOD .EQ. 2) .OR. (MOD .EQ. 5)) THEN
C Nonsymmetric matrix:
ISTYPE = 2
ELSE
C Hermitian matrix:
ISTYPE = 3
END IF
C
C Get a temporary vector for generating columns and RHS's:
N_TEMP = IDTYPE*NUMEQ
CALL FMSRMG (RMD, L_TEMP, N_TEMP)
C
C (2) Open FMS files:
C
C Matrix files:
C
C [A11]:
IF(N1 .GT. 0) THEN
IF (MOD.EQ.1) THEN
CALL RSDI (LOWEQ, N1 , 'LU_A11', LU_A11)
ELSE IF(MOD.EQ.2) THEN
CALL RNDI (LOWEQ, N1 , 'LU_A11', LU_A11)
ELSE IF(MOD.EQ.3) THEN
CALL CHDI (LOWEQ, N1 , 'LU_A11', LU_A11)
ELSE IF(MOD.EQ.4) THEN
CALL CSDI (LOWEQ, N1 , 'LU_A11', LU_A11)
ELSE
CALL CNDI (LOWEQ, N1 , 'LU_A11', LU_A11)
END IF
CALL FMSOM (LU_A11, 'LU_F11', LU_F11)
END IF
C
C [A22]:
IF(N22 .GT. 0) THEN
IF (MOD.EQ.1) THEN
CALL RSDI (LOWEQ, N22, 'LU_A22', LU_A22)
ELSE IF(MOD.EQ.2) THEN
CALL RNDI (LOWEQ, N22, 'LU_A22', LU_A22)
ELSE IF(MOD.EQ.3) THEN
CALL CHDI (LOWEQ, N22, 'LU_A22', LU_A22)
ELSE IF(MOD.EQ.4) THEN
CALL CSDI (LOWEQ, N22, 'LU_A22', LU_A22)
ELSE
CALL CNDI (LOWEQ, N22, 'LU_A22', LU_A22)
END IF
CALL FMSOM (LU_A22, 'LU_F22', LU_F22)
END IF
C
C Vector and off-diagonal matrix files:
IF(N1 .GT. 0) THEN
IF(ISTYPE .NE. 1) THEN
IF(N2 .GT. 0) THEN
CALL FMSOV (N1, IDTYPE, N2, 'LU_A12', LU_A12)
END IF
IF(N3 .GT. 0) THEN
CALL FMSOV (N1, IDTYPE, N3, 'LU_A13', LU_A13)
END IF
END IF
IF(N2 .GT. 0) THEN
CALL FMSOV (N1, IDTYPE, N2, 'LU_A21', LU_A21)
CALL FMSOV (N1, IDTYPE, N2, 'LU_F12', LU_F12)
END IF
IF(N3 .GT. 0) THEN
CALL FMSOV (N1, IDTYPE, N3, 'LU_A31', LU_A31)
CALL FMSOV (N1, IDTYPE, N3, 'LU_F13', LU_F13)
END IF
CALL FMSOV (N1 , IDTYPE, NRHS, 'LU_B1' , LU_B1 )
CALL FMSOV (N1 , IDTYPE, NRHS, 'LU_X1' , LU_X1 )
END IF
IF(N22 .GT. 0) THEN
CALL FMSOV (N22, IDTYPE, NRHS, 'LU_B2' , LU_B2)
CALL FMSOV (N22, IDTYPE, NRHS, 'LU_X2' , LU_X2)
END IF
C
C (3) Write data to FMS files:
C
C Generate [A11]:
IF(N1 .GT. 0) THEN
CALL AGEN (RMD(L_TEMP), IDTYPE, 1, N1, 1, N1, NUMEQ, MTEST,
1 LU_A11, ISTYPE)
C
IF(ISTYPE .NE. 1) THEN
C
C Generate [A12]:
IF(N2 .GT. 0) THEN
CALL AGEN (RMD(L_TEMP), IDTYPE, 1, N1, N1+1, N1+N2,
1 NUMEQ, MTEST, LU_A12, ISTYPE)
END IF
C
C Generate [A13]:
IF(N3 .GT. 0) THEN
CALL AGEN (RMD(L_TEMP), IDTYPE, 1, N1, N1+N2+1, NUMEQ,
1 NUMEQ, MTEST, LU_A13, ISTYPE)
END IF
END IF
C
C Generate [A21]:
IF(N2 .GT. 0) THEN
CALL AGEN (RMD(L_TEMP), IDTYPE, N1+1, N1+N2, 1, N1,
1 NUMEQ, MTEST, LU_A21, ISTYPE)
END IF
C
C Generate [A31]:
IF(N3 .GT. 0) THEN
CALL AGEN (RMD(L_TEMP), IDTYPE, N1+N2+1, NUMEQ, 1, N1,
1 NUMEQ, MTEST, LU_A31, ISTYPE)
END IF
C
C Generate {B1}:
CALL BGEN (RMD(L_TEMP), IDTYPE, 1, N1, NRHS, LU_B1)
END IF
C
C Generate [A22]:
IF(N22 .GT. 0) THEN
CALL AGEN (RMD(L_TEMP), IDTYPE, N1+1, NUMEQ, N1+1, NUMEQ,
1 NUMEQ, MTEST, LU_A22, ISTYPE)
C
C Generate {B2}:
CALL BGEN (RMD(L_TEMP), IDTYPE, N1+1, NUMEQ, NRHS, LU_B2)
END IF
C
C (4) Perform matrix algebra:
C
C (4.1) Factor [A11] into [F11]:
IF(N1 .GT. 0) THEN
WRITE(LUPR,2001)
IF (MOD .EQ. 1) THEN
CALL RSDAF
1 (LU_A11, R_ONE, 1, LUS0, 0, LUA0, LU_F11, LU_B1, LU_X1, 0)
ELSE IF(MOD .EQ. 2) THEN
CALL RNDAF
1 (LU_A11, R_ONE, 1, LUS0, 0, LUA0, LU_F11, LU_B1, LU_X1, 0)
ELSE IF(MOD .EQ. 3) THEN
CALL CHDAF
1 (LU_A11, R_ONE, 1, LUS0, 0, LUA0, LU_F11, LU_B1, LU_X1, 0)
ELSE IF(MOD .EQ. 4) THEN
CALL CSDAF
1 (LU_A11, C_ONE, 1, LUS0, 0, LUA0, LU_F11, LU_B1, LU_X1, 0)
ELSE
CALL CNDAF
1 (LU_A11, C_ONE, 1, LUS0, 0, LUA0, LU_F11, LU_B1, LU_X1, 0)
END IF
END IF
C
C (4.2) Solve [F11]{F12} = {A12}
IF( (N1 .GT. 0) .AND.
1 (N2 .GT. 0) ) THEN
WRITE(LUPR,2002)
IF (MOD .EQ. 1) THEN
CALL RSDS (LU_F11, LU_A21, LU_F12, N2, 0)
ELSE IF(MOD .EQ. 2) THEN
CALL RNDS (LU_F11, LU_A21, LU_F12, N2, 0)
ELSE IF(MOD .EQ. 3) THEN
CALL CHDS (LU_F11, LU_A12, LU_F12, N2, 0)
ELSE IF(MOD .EQ. 4) THEN
CALL CSDS (LU_F11, LU_A21, LU_F12, N2, 0)
ELSE
CALL CNDS (LU_F11, LU_A12, LU_F12, N2, 0)
END IF
END IF
C
C (4.3) Solve [F11]{F13} = {A13}
IF( (N1 .GT. 0) .AND.
1 (N3 .GT. 0) ) THEN
WRITE(LUPR,2003)
IF (MOD .EQ. 1) THEN
CALL RSDS (LU_F11, LU_A31, LU_F13, N3, 0)
ELSE IF(MOD .EQ. 2) THEN
CALL RNDS (LU_F11, LU_A31, LU_F13, N3, 0)
ELSE IF(MOD .EQ. 3) THEN
CALL CHDS (LU_F11, LU_A13, LU_F13, N3, 0)
ELSE IF(MOD .EQ. 4) THEN
CALL CSDS (LU_F11, LU_A31, LU_F13, N3, 0)
ELSE
CALL CNDS (LU_F11, LU_A13, LU_F13, N3, 0)
END IF
END IF
C
C (4.4) Solve [F11]{X1} = {B1}
IF(N1 .GT. 0) THEN
WRITE(LUPR,2004)
IF (MOD .EQ. 1) THEN
CALL RSDS (LU_F11, LU_B1, LU_X1, NRHS, 0)
ELSE IF(MOD .EQ. 2) THEN
CALL RNDS (LU_F11, LU_B1, LU_X1, NRHS, 0)
ELSE IF(MOD .EQ. 3) THEN
CALL CHDS (LU_F11, LU_B1, LU_X1, NRHS, 0)
ELSE IF(MOD .EQ. 4) THEN
CALL CSDS (LU_F11, LU_B1, LU_X1, NRHS, 0)
ELSE
CALL CNDS (LU_F11, LU_B1, LU_X1, NRHS, 0)
END IF
END IF
C
C (4.5) Multiply [F22] = [A22] - [A21]T[F12]
IF(N22 .GT. 0) THEN
IF( (N1 .GT. 0) .AND.
1 (N2 .GT. 0) ) THEN
WRITE(LUPR,2005)
CALL FMSMM (LU_F22, LU_A22, -1, TRANSA, LU_A21, LU_F12)
ELSE
CALL FMSCPY (LU_A22, ' ', LU_F22)
END IF
END IF
C
C (4.6) Multiply [F22] = [F22] - [A21]T[F13]
IF( (N1 .GT. 0) .AND.
1 (N3 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2006)
IF(ISTYPE .NE. 1) THEN
CALL FMSIST ('MMCOL', N2)
CALL FMSMM (LU_F22, LU_F22, -1, TRANSA, LU_A21, LU_F13)
CALL FMSIST ('MMCOL', 0)
END IF
END IF
C
C (4.7) Multiply [F22] = [F22] - [A31]T[F12]
IF( (N1 .GT. 0) .AND.
1 (N2 .GT. 0) .AND.
2 (N3 .GT. 0) .AND.
3 (N22 .GT. 0) ) THEN
WRITE(LUPR,2007)
CALL FMSIST ('MMROW', N2)
CALL FMSMM (LU_F22, LU_F22, -1, TRANSA, LU_A31, LU_F12)
END IF
C
C (4.8) Multiply [F22] = [F22] - [A31]T[F13]
IF( (N1 .GT. 0) .AND.
1 (N3 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2008)
CALL FMSIST ('MMCOL', N2)
CALL FMSMM (LU_F22, LU_F22, -1, 'T', LU_A31, LU_F13)
CALL FMSIST ('MMCOL', 0)
CALL FMSIST ('MMROW', 0)
END IF
C
C (4.9) Multiply {X2} = {B2} - [A21]T{X1}
IF(N22 .GT. 0) THEN
IF( (N1 .GT. 0) .AND.
1 (N2 .GT. 0) ) THEN
WRITE(LUPR,2009)
CALL FMSMM (LU_X2 , LU_B2 , -1, TRANSA, LU_A21, LU_X1 )
ELSE
CALL FMSCPY (LU_B2, ' ', LU_X2)
END IF
END IF
C
C (4.10)Multiply {X2} = {X2} - [A31]{X1}
IF( (N1 .GT. 0) .AND.
1 (N3 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2010)
CALL FMSIST ('MMROW', N2)
CALL FMSMM (LU_X2, LU_X2, -1, TRANSA, LU_A31, LU_X1 )
CALL FMSIST ('MMROW', 0)
END IF
C
C (4.11)Factor [A22] into [F22]
IF(N22 .GT. 0) THEN
WRITE(LUPR,2011)
IF (MOD .EQ. 1) THEN
CALL RSDAF
1 (LU_F22, R_ONE, 1, LUS0, 0, LUA0, LU_F22, LU_X2, LU_X2, 0)
ELSE IF(MOD .EQ. 2) THEN
CALL RNDAF
1 (LU_F22, R_ONE, 1, LUS0, 0, LUA0, LU_F22, LU_X2, LU_X2, 0)
ELSE IF(MOD .EQ. 3) THEN
CALL CHDAF
1 (LU_F22, R_ONE, 1, LUS0, 0, LUA0, LU_F22, LU_X2, LU_X2, 0)
ELSE IF(MOD .EQ. 4) THEN
CALL CSDAF
1 (LU_F22, C_ONE, 1, LUS0, 0, LUA0, LU_F22, LU_X2, LU_X2, 0)
ELSE IF(MOD .EQ. 5) THEN
CALL CNDAF
1 (LU_F22, C_ONE, 1, LUS0, 0, LUA0, LU_F22, LU_X2, LU_X2, 0)
END IF
END IF
C
C (4.l2)Solve [F22]{X2} = {X2'}
IF(N22 .GT. 0) THEN
WRITE(LUPR,2012)
IF (MOD .EQ. 1) THEN
CALL RSDS (LU_F22, LU_X2, LU_X2, NRHS, 0)
ELSE IF(MOD .EQ. 2) THEN
CALL RNDS (LU_F22, LU_X2, LU_X2, NRHS, 0)
ELSE IF(MOD .EQ. 3) THEN
CALL CHDS (LU_F22, LU_X2, LU_X2, NRHS, 0)
ELSE IF(MOD .EQ. 4) THEN
CALL CSDS (LU_F22, LU_X2, LU_X2, NRHS, 0)
ELSE
CALL CNDS (LU_F22, LU_X2, LU_X2, NRHS, 0)
END IF
END IF
C
C (4.13)Multiply {X1} = {X1} - [F13]{X2}
IF( (N1 .GT. 0) .AND.
1 (N3 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2013)
CALL FMSIST ('MMKA', N2)
CALL FMSMM (LU_X1, LU_X1, -1, 'N', LU_F13, LU_X2)
CALL FMSIST ('MMKA', 0)
END IF
C
C (4.14)Multiply {X1} = {X1} - [F12]{X2}
IF( (N1 .GT. 0) .AND.
1 (N2 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2014)
CALL FMSMM (LU_X1, LU_X1, -1, 'N', LU_F12, LU_X2)
END IF
C
C Compute the error in {X1}:
C {E1} = {B1} - [A11]{X1} - [A12|A13]{X2}:
C
C (4.15)Multiply {B1} = {B1} - [A11]{X1}
IF(N1 .GT. 0) THEN
WRITE(LUPR,2015)
CALL FMSMM (LU_B1, LU_B1, -1, 'N', LU_A11, LU_X1 )
END IF
C
C (4.16)Multiply {B1} = {B1} - [A12]{X2}
IF( (N1 .GT. 0) .AND.
1 (N2 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2016)
IF(ISTYPE .NE. 1) THEN
CALL FMSMM (LU_B1, LU_B1, -1, 'N', LU_A12, LU_X2 )
ELSE
CALL FMSMM (LU_B1, LU_B1, -1, 'N', LU_A21, LU_X2 )
END IF
END IF
C
C (4.17)Multiply {B1} = {B1} - [A13]{X2}
IF( (N1 .GT. 0) .AND.
1 (N3 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2017)
CALL FMSIST ('MMKA', N2)
IF(ISTYPE .NE. 1) THEN
CALL FMSMM (LU_B1, LU_B1, -1, 'N', LU_A13, LU_X2 )
ELSE
CALL FMSMM (LU_B1, LU_B1, -1, 'N', LU_A31, LU_X2 )
END IF
CALL FMSIST ('MMKA', 0)
END IF
C
C Compute the error in {X2}:
C {E2} = {B2} - [A22]{X2} - [A21]{X1} - [A31]{X1}:
C
C (4.18)Multiply {B2} = {B2} - [A22]{X2}
IF(N22 .GT. 0) THEN
WRITE(LUPR,2018)
CALL FMSMM (LU_B2, LU_B2, -1, 'N', LU_A22, LU_X2)
END IF
C
C (4.19)Multiply {B2} = {B2} - [A21]T{X1}
IF( (N1 .GT. 0) .AND.
1 (N2 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2019)
CALL FMSMM (LU_B2, LU_B2, -1, TRANSA, LU_A21, LU_X1)
END IF
C
C (4.20)Multiply {B2} = {B2} - [A31]T{X1}
IF( (N1 .GT. 0) .AND.
1 (N3 .GT. 0) .AND.
2 (N22 .GT. 0) ) THEN
WRITE(LUPR,2020)
CALL FMSIST ('MMROW', N2)
CALL FMSMM (LU_B2, LU_B2, -1, TRANSA, LU_A31, LU_X1)
CALL FMSIST ('MMROW', 0)
END IF
C
C (5) Read data from FMS files:
C
IF(N1 .GT. 0) THEN
C Check {X1}:
ERROR = R_ZERO
C
C Initialize FMSGET:
CALL FMSGET (LU_B1, 0, 0, 0, 0, RMD, 0)
DO IRHS = 1,NRHS
CALL FMSGET (LU_B1, N1, 1, 1, IRHS, RMD(L_TEMP), N1)
L = L_TEMP
DO I = 1,N1
IF(IDTYPE .EQ. 1) THEN
EI = ABS(RMD(L))
L = L + 1
ELSE
EI = ABS(DCMPLX(RMD(L),RMD(L+1)))
L = L + 2
END IF
IF(EI .GT. ERROR) ERROR = EI
END DO
END DO
C
C End FMSGET:
CALL FMSGET (LU_B1, 0, 0, N1+1, 0, RMD, 0)
WRITE(6,*) 'MAXIMUM ERROR IN {X1} =', ERROR
END IF
C
IF(N22 .GT. 0) THEN
C Check {X2}:
ERROR = R_ZERO
C
C Initialize FMSGET:
CALL FMSGET (LU_B2, 0, 0, 0, 0, RMD, 0)
DO IRHS = 1,NRHS
CALL FMSGET (LU_B2, N22, 1, 1, IRHS, RMD(L_TEMP), N22)
L = L_TEMP
DO I = 1,N22
IF(IDTYPE .EQ. 1) THEN
EI = ABS(RMD(L))
L = L + 1
ELSE
EI = ABS(DCMPLX(RMD(L),RMD(L+1)))
L = L + 2
END IF
IF(EI .GT. ERROR) ERROR = EI
END DO
END DO
C
C End FMSGET:
CALL FMSGET (LU_B2, 0, 0, N22+1, 0, RMD, 0)
WRITE(6,*) 'MAXIMUM ERROR IN {X2} =', ERROR
END IF
C
C (6) Close FMS files:
IF(N1 .GT. 0) THEN
CALL FMSCM (LU_A11)
CALL FMSCM (LU_F11)
IF(ISTYPE .NE. 1) THEN
IF(N2 .GT. 0) CALL FMSCV (LU_A12)
IF(N3 .GT. 0) CALL FMSCV (LU_A13)
END IF
IF(N2 .GT. 0) THEN
CALL FMSCV (LU_A21)
CALL FMSCV (LU_F12)
END IF
IF(N3 .GT. 0) THEN
CALL FMSCV (LU_A31)
CALL FMSCV (LU_F13)
END IF
CALL FMSCV (LU_B1 )
CALL FMSCV (LU_X1 )
END IF
C
IF(N22 .GT. 0) THEN
CALL FMSCM (LU_A22)
CALL FMSCM (LU_F22)
CALL FMSCV (LU_B2 )
CALL FMSCV (LU_X2 )
END IF
CALL FMSRMR (RMD, L_TEMP, N_TEMP)
IF(ASK('Do you want another solution?')) GO TO 10
CALL FMSPOP(MYNAME)
CALL FMSEND
2001 FORMAT (/
1 ' ( 1) Factor [A11] into [F11]'/
2 ' ==============================')
2002 FORMAT (/
1 ' ( 2) Solve [F11]{F12} = {A12}'/
2 ' ================================')
2003 FORMAT (/
1 ' ( 3) Solve [F11]{F13} = {A13}'/
2 ' ================================')
2004 FORMAT (/
1 ' ( 4) Solve [F11]{X1a} = {B1}'/
2 ' ===============================')
2005 FORMAT (/
1 ' ( 5) Multiply [F22a] = [A22] - [A21]T[F12]'/
2 ' ===========================================')
2006 FORMAT (/
1 ' ( 6) Multiply [F22b] = [F22a] - [A21]T[F13]'/
2 ' ===========================================')
2007 FORMAT (/
1 ' ( 7) Multiply [F22c] = [F22b] - [A31]T[F12]'/
2 ' ===========================================')
2008 FORMAT (/
1 ' ( 8) Multiply [F22d] = [F22d] - [A31]T[F13]'/
2 ' ===========================================')
2009 FORMAT (/
1 ' ( 9) Multiply {X2a} = {B2} - [A21]T{X1a}'/
2 ' ===========================================')
2010 FORMAT (/
1 ' (10) Multiply {X2b} = {X2a} - [A31]T{X1a}'/
2 ' ===========================================')
2011 FORMAT (/
1 ' (11) Factor [F22d] into [F22]'/
2 ' ===============================')
2012 FORMAT (/
1 ' (12) Solve [F22]{X2} = {X2b}'/
2 ' ===============================')
2013 FORMAT (/
1 ' (13) Multiply {X1b} = {X1a} - [F13]{X2}'/
2 ' =========================================')
2014 FORMAT (/
1 ' (14) Multiply {X1} = {X1b} - [F12]{X2}'/
2 ' =========================================')
2015 FORMAT (/
1 ' (15) Multiply {E1a} = {B1} - [A11]{X1}'/
2 ' =========================================')
2016 FORMAT (/
1 ' (16) Multiply {E1b} = {E1a} - [A12]{X2}'/
2 ' =========================================')
2017 FORMAT (/
1 ' (17) Multiply {E1} = {E1b} - [A13]{X2}'/
2 ' =========================================')
2018 FORMAT (/
1 ' (18) Multiply {E2a} = {B2} - [A22]{X2}'/
2 ' =========================================')
2019 FORMAT (/
1 ' (19) Multiply {E2b} = {E2a} - [A21]{X1}'/
2 ' =========================================')
2020 FORMAT (/
1 ' (20) Multiply {E2} = {E2b} - [A31]{X1}'/
2 ' =========================================')
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 AGEN (A, IDTYPE, IROW1, IROW2, JCOL1, JCOL2,
1 NUMEQ, MTEST, LUA, ISTYPE)
C=======================================================================
REAL*8 A(IDTYPE,NUMEQ)
INTEGER IROW1, IROW2, JCOL1, JCOL2, NUMEQ
INTEGER LUA(25)
REAL*8 A_1, A_DIAG
REAL*8 R_ZERO, R_ONE
DATA R_ZERO/0.0D0/
DATA R_ONE /1.0D0/
C
WRITE(6,2000) IROW1, IROW2, JCOL1, JCOL2
NROWS = IROW2 - IROW1 + 1
NCOLS = JCOL2 - JCOL1 + 1
IEDGE = 0
IF(IROW1 .GT. JCOL1) THEN
C This is a off-diagonal block in the lower triangle:
C Write it out transposed:
LDU = 1
I1 = JCOL1
I2 = JCOL2
J1 = IROW1
J2 = IROW2
IF(JCOL1 .EQ. 1) IEDGE = 1
ELSE IF(IROW1 .EQ. JCOL1) THEN
C This block is on the diagonal:
C Write it out by columns:
LDU = 2
I1 = IROW1
I2 = IROW2
J1 = JCOL1
J2 = JCOL2
IF(IROW1 .EQ. 1) IEDGE = 1
ELSE
C This block is in the upper triangle:
C Write it out by columns:
LDU = 3
I1 = IROW1
I2 = IROW2
J1 = JCOL1
J2 = JCOL2
IF(IROW1 .EQ. 1) IEDGE = 1
END IF
LENV = I2 - I1 + 1
C
IF(MTEST .EQ. 1) THEN
C Fixed data:
C Initialize A to default values (generic column):
DO I = 1,IDTYPE
DO J=1,LENV
A(I,J) = R_ZERO
END DO
END DO
A_DIAG = R_ONE
IF(IEDGE .EQ. 1) THEN
A_1 = -R_ONE
ELSE
A_1 = R_ZERO
END IF
A(1,1) = A_1
END IF
C
C Initialize FMSPUT:
CALL FMSPUT (LUA, 0, 0, 0, 0, A, 0)
C
C Loop over the columns:
DO J = J1,J2
JC1 = J - J1 + 1
IF(MTEST .EQ. 1) THEN
C Use fixed data:
IF(LDU .EQ. 2) THEN
C Diagonal block:
IF(J .EQ. 1) THEN
C Special case. This is the first global column:
DO I=1,NROWS
A(1,I) = -R_ONE
END DO
A_DIAG = DFLOAT(NUMEQ)
END IF
A(1,JC1) = A_DIAG
END IF
ELSE
C Use random data:
IF(LDU .EQ. 1) THEN
C Lower triangle off-diagonal block; Generate transposed
CALL VRAN (A, IDTYPE, J, J, JCOL1, JCOL2,
1 NUMEQ, ISTYPE)
ELSE
C Diagonal or upper triangle block:
CALL VRAN (A, IDTYPE, IROW1, IROW2, J, J,
1 NUMEQ, ISTYPE)
END IF
END IF
D5001 FORMAT (' CALL FMSPUT(LUA,',I5,', 1, 1,',I5,', A,',I5,')')
CALL FMSPUT (LUA, LENV, 1, 1, JC1, A, LENV)
IF(MTEST .EQ. 1) THEN
C Restore A to default:
IF(LDU .EQ. 2) THEN
C Diagonal block:
IF(J .EQ. 1) THEN
C Special case. This is the first global column:
DO I=1,NROWS
A(1,I) = R_ZERO
END DO
A_DIAG = R_ONE
END IF
A(1,JC1) = R_ZERO
A(1,1) = A_1
END IF
END IF
END DO
C
C End FMSPUT:
NEND = LENV+1
CALL FMSPUT (LUA, 0, 0, NEND, 0, A, 0)
RETURN
2000 FORMAT (' Populating A(',I5,':',I5,',',I5,':',I5,')')
END
C=======================================================================
SUBROUTINE VRAN (A, IDTYPE, IROW1, IROW2, JCOL1, JCOL2, NUMEQ,
1 ISTYPE)
C=======================================================================
REAL*8 A(IDTYPE, IROW1:IROW2,JCOL1:JCOL2)
INTEGER IROW1, IROW2, JCOL1, JCOL2, NUMEQ, ISTYPE
REAL*8 RANAIJ, SUM
INTEGER ISEEDR, ISEEDC
PARAMETER (ISEEDR=11111111)
PARAMETER (ISEEDC=12345678)
DO JCOL = JCOL1,JCOL2
DO IROW = IROW1,IROW2
IF(JCOL .LT. IROW) THEN
C This term is in the lower triangle:
IY = IROW*ISEEDR + JCOL*ISEEDC
IY = MOD(IY,65536)
IY = IY * 11111 + 1
RANAIJ = DFLOAT((IY/256)+1) / 2844374.0D0
RANAIJ = RANAIJ - 0.5D0
ELSE IF(JCOL .EQ. IROW) THEN
C This term is on the diagonal:
C Set the value to (-) the sum of the off-diagonal terms:
SUM = 0.0D0
DO II = 1,NUMEQ
IF(II .NE. IROW) THEN
IF(II .LT. IROW) THEN
C This term is in the lower triangle:
IY = IROW*ISEEDR + II*ISEEDC
ELSE
C This term is in the upper triangle:
IF(ISTYPE .EQ. 2) THEN
C This matrix is nonsymmetric:
IY = IROW*ISEEDR + II*ISEEDC
ELSE
C This matrix is symmetric or Hermitian:
IY = II*ISEEDR + IROW*ISEEDC
END IF
END IF
IY = MOD(IY,65536)
IY = IY * 11111 + 1
RANAIJ = DFLOAT((IY/256)+1) / 2844374.0D0
RANAIJ = RANAIJ - 0.5D0
SUM = SUM + RANAIJ
END IF
END DO
RANAIJ = -SUM
ELSE
C This term is in the upper triangle:
IF(ISTYPE .EQ. 2) THEN
C This matrix is nonsymmetric:
IY = IROW*ISEEDR + JCOL*ISEEDC
ELSE
C This matrix is symmetric or Hermitian:
IY = JCOL*ISEEDR + IROW*ISEEDC
END IF
IY = MOD(IY,65536)
IY = IY * 11111 + 1
RANAIJ = DFLOAT((IY/256)+1) / 2844374.0D0
RANAIJ = RANAIJ - 0.5D0
END IF
A(1,IROW,JCOL) = RANAIJ
IF((IROW .EQ. 1) .AND.
1 (JCOL .EQ. 1)) A(1,1,1) = RANAIJ + 1.0D0
IF(IDTYPE .EQ. 2) THEN
IF(ISTYPE .EQ. 3) THEN
C This matrix is Hermitian:
IF(JCOL .LT. IROW) THEN
A(2,IROW,JCOL) = RANAIJ
ELSE IF(JCOL .EQ. IROW) THEN
A(1,IROW,JCOL) = 2*RANAIJ
A(2,IROW,JCOL) = 0.0D0
ELSE
A(2,IROW,JCOL) = -RANAIJ
END IF
ELSE
C This matrix is Complex Symmetric or Nonsymmetric:
A(2,IROW,JCOL) = RANAIJ
END IF
END IF
END DO
END DO
RETURN
END
C=======================================================================
SUBROUTINE BGEN (RHS, IDTYPE, IROW1, IROW2, NRHS, LUB)
C=======================================================================
REAL*8 RHS(IDTYPE, *)
INTEGER IDTYPE, IROW1, IROW2, NRHS
INTEGER LUB(25)
REAL*8 R_ZERO, R_ONE
DATA R_ZERO/0.0D0/
DATA R_ONE /1.0D0/
C
WRITE(6,2000) IROW1, IROW2, NRHS
C
NROWS = IROW2 - IROW1 + 1
DO I=1,IDTYPE
DO J=1,NROWS
RHS(I,J) = R_ZERO
END DO
END DO
IF(IROW1 .EQ. 1) RHS(1,1) = R_ONE
C
C Initialize FMSPUT:
CALL FMSPUT (LUB, 0, 0, 0, 0, RHS, 0)
C
C Write out columns of {B}:
DO IRHS = 1,NRHS
CALL FMSPUT (LUB, NROWS, 1, 1, IRHS, RHS, NROWS)
END DO
C
C End FMSPUT:
CALL FMSPUT (LUB, 0, 0, NROWS+1, 0, RHS, 0)
RETURN
2000 FORMAT (' Populating B(',I5,':',I5,',',I5,')')
END