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>|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) D CALL FMSPRM (LU_A11, 1) 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) D CALL FMSPRV (LU_A12, 1, N2, 10) 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) D CALL FMSPRV (LU_A13, 1, N3, 10) 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) D CALL FMSPRV (LU_A21, 1, N2, 10) 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) D CALL FMSPRV (LU_A31, 1, N3, 10) END IF C C Generate {B1}: CALL BGEN (RMD(L_TEMP), IDTYPE, 1, N1, NRHS, LU_B1) D CALL FMSPRV (LU_B1, 1, NRHS, 10) 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) D CALL FMSPRM(LU_A22, 1) C C Generate {B2}: CALL BGEN (RMD(L_TEMP), IDTYPE, N1+1, NUMEQ, NRHS, LU_B2) D CALL FMSPRV (LU_B2, 1, NRHS, 10) 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_A12, 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_A13, 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, TRANSA, 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 D WRITE(6,5001) LENV, JC1, LENV 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 D PRINT *,'Restoring first column' 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 D WRITE(6,2001) IROW, SUM 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 D WRITE(6,2000) IROW, JCOL, RANAIJ 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 D2000 FORMAT (' VRAN: A(',I5,',',I5,')....=',F10.5) D2001 FORMAT (' VRAN: SUM OF ROW',I5, '...=',F10.5) 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