C T E S T R E S O L V
C
C This program tests subroutine RESOLV.
C The test matrix is of the form
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
C Program name:
CHARACTER*10 MYNAME
PARAMETER (MYNAME='EXAMPLE_14')
C
C FMS matrix and vector file attributes:
INTEGER LU_A11(25)
INTEGER LU_A12(25)
INTEGER LU_A21(25)
INTEGER LU_A22(25), LU_F22(25)
INTEGER LU_B1 (25), LU_X1(25)
INTEGER LU_B2 (25), LU_S2(25), LU_X2(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 Data Type:
INTEGER IDTYPE
PARAMETER (IDTYPE=2)
C
C Input Data:
LOGICAL ASK
INTEGER ASK_I
C
C LOCAL VARIABLES:
C Number of equations in [A11]:
INTEGER N1
C Number of equations in [A22]:
INTEGER N2
C Maximum of N1 and N2:
INTEGER NMAX
C Number of solution vectors:
INTEGER NRHS
C Form a new [A11] matrix:
LOGICAL NEWA11
C Form a new [A12] matrix:
LOGICAL NEWA12
C Form a new [A21] matrix:
LOGICAL NEWA21
C Form a new [A22] matrix:
LOGICAL NEWA22
C Form a new {B1} matrix:
LOGICAL NEWB1
C Form a new {B2} matrix:
LOGICAL NEWB2
C Check the values of {X2}:
LOGICAL ICK_X2
C Error on current term:
REAL*8 EI
C Maximum overall error:
REAL*8 E2MAX, ERROR
C CMD Pointer to scratch vector:
INTEGER L_TEMP
C
C Profile vector for a full matrix:
INTEGER LOWEQ(1)
DATA LOWEQ/-1/
C
C CONSTANTS:
COMPLEX*16 C_ONE
DATA C_ONE/(1.0D0,0.0D0)/
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH (MYNAME)
C
C Use full column partial pivoting (SLAB format):
CALL FMSIST ('MFMAT', 3)
C
C Read in problem size parameters:
N1 = ASK_I('Enter the number of equations in A11')
N2 = ASK_I('Enter the number of equations in A22')
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
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
C
C Get a temporary vector for generating columns and RHS's:
NMAX = MAX0(N1,N2)
CALL FMSCMG (CMD, L_TEMP, NMAX)
C
C (2) Open FMS files:
CALL CNDI (LOWEQ, N1, 'LU_A11', LU_A11)
CALL FMSOV (N1, IDTYPE, N2 , 'LU_A12', LU_A12)
CALL FMSOV (N1, IDTYPE, N2 , 'LU_A21', LU_A21)
CALL CNDI (LOWEQ, N2, 'LU_A22', LU_A22)
CALL FMSOM (LU_A22, 'LU_F22', LU_F22)
CALL FMSOV (N1, IDTYPE, NRHS, 'LU_B1' , LU_B1 )
CALL FMSOV (N1, IDTYPE, NRHS, 'LU_X1' , LU_X1 )
CALL FMSOV (N2, IDTYPE, NRHS, 'LU_B2' , LU_B2 )
CALL FMSOV (N2, IDTYPE, NRHS, 'LU_S2' , LU_S2 )
CALL FMSOV (N2, IDTYPE, NRHS, 'LU_X2' , LU_X2 )
C
C (3) Write data to FMS files:
NEWA11 = .TRUE.
NEWA12 = .TRUE.
NEWA21 = .TRUE.
NEWA22 = .TRUE.
NEWB1 = .TRUE.
NEWB2 = .TRUE.
ICK_X2 = .TRUE.
C
C Branch back to here for refinement solutions:
30 CONTINUE
C
C Generate [A11]:
IF(NEWA11) CALL A11GEN (CMD(L_TEMP), N1, N2, LU_A11)
C
C Generate [A12]:
IF(NEWA12) CALL A12GEN (CMD(L_TEMP), N1, N2, LU_A12)
C
C Generate [A21]:
IF(NEWA21) CALL A21GEN (CMD(L_TEMP), N1, N2, LU_A21)
C
C Generate [A22]:
IF(NEWA22) CALL A22GEN (CMD(L_TEMP), N2, LU_A22)
C
C Generate {B1}:
IF(NEWB1 ) CALL B1GEN (CMD(L_TEMP), N1, NRHS, LU_B1 )
C
C Generate {B2}:
IF(NEWB2 ) CALL B2GEN (CMD(L_TEMP), N2, NRHS, LU_B2 )
C
C (4) Perform matrix algebra:
CALL RESOLV
1 (LU_A11, LU_A11,
2 LU_A12, LU_A12,
3 LU_A21,
4 LU_A22, LU_F22,
5 LU_B1, LU_B1, LU_X1,
6 LU_B2, LU_S2, LU_X2,
7 NEWA11, NEWA12,
8 NEWA21, NEWA22,
9 NEWB1, NEWB2,
1 ICK_X2, E2MAX)
C
C (5) Read data from FMS files:
C Check the answer:
ERROR = E2MAX
C
C Initialize FMSGET:
CALL FMSGET (LU_X1, 0, 0, 0, 0, CMD, 0)
DO 60 IRHS = 1,NRHS
C
C Check {X1}:
CALL FMSGET (LU_X1, N1, 1, 1, IRHS, CMD(L_TEMP), N1)
DO 51 I = 0,(N1-1)
EI = ABS(CMD(L_TEMP+I) - C_ONE)
IF(EI .GT. ERROR) ERROR = EI
51 CONTINUE
60 CONTINUE
C
C End FMSGET:
CALL FMSGET (LU_X1, 0, 0, N1+1, 1, CMD, 0)
WRITE(6,*) 'MAXIMUM ERROR =', ERROR
C
C Modify the matrix and solve it again:
IF(.NOT. ASK('Do you want another solution?')) GO TO 100
NEWA11 = ASK('Do you want a new A11?')
NEWA12 = ASK('Do you want a new A12?')
NEWA21 = ASK('Do you want a new A21?')
NEWA22 = ASK('Do you want a new A22?')
NEWB1 = ASK('Do you want a new B1?' )
NEWB2 = ASK('Do you want a new B2?' )
C
C We have chosen to overlay factored data on input matrices.
C Therefore if new factors must be recomputed, the arrays must
C be reset to the unfactored values.
C Based on the options selected above, these conditions become:
IF( NEWA11 ) THEN
NEWA12 = .TRUE.
NEWA22 = .TRUE.
NEWB1 = .TRUE.
NEWB2 = .TRUE.
END IF
IF( NEWA12 ) NEWA22 = .TRUE.
IF( NEWB1 ) NEWB2 = .TRUE.
IF( NEWA11 .OR.
1 NEWA12 .OR.
2 NEWA21 .OR.
3 NEWA22 .OR.
4 NEWB1 .OR.
5 NEWB2) GO TO 30
C
C (6) Close FMS files:
100 CONTINUE
CALL FMSCM (LU_A11)
CALL FMSCV (LU_A12)
CALL FMSCV (LU_A21)
CALL FMSCM (LU_A22)
CALL FMSCM (LU_F22)
CALL FMSCV (LU_B1 )
CALL FMSCV (LU_X1 )
CALL FMSCV (LU_B2 )
CALL FMSCV (LU_S2 )
CALL FMSCV (LU_X2 )
CALL FMSCMR (CMD, L_TEMP, NMAX)
CALL FMSPOP (MYNAME)
CALL FMSEND
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 A11GEN (ACOL, N1, N2, LU_A11)
C=======================================================================
COMPLEX*16 ACOL(N1)
INTEGER N1, N2
INTEGER LU_A11(25)
CHARACTER*6 MYNAME
PARAMETER (MYNAME='A11GEN')
REAL*8 R_ZERO, A11_Re
DATA R_ZERO/0.0D0/
COMPLEX*16 C_MONE, C_ZERO, C_ONE
DATA C_MONE/(-1.0D0,0.0D0)/
DATA C_ZERO/( 0.0D0,0.0D0)/
DATA C_ONE /( 1.0D0,0.0D0)/
C
CALL FMSPSH (MYNAME)
WRITE(6,2000)
C
C Initialize FMSCOL
CALL FMSCOL (0, ACOL, LU_A11)
C
C Write first column
DO 10 I = 2,N1
ACOL(I) = C_MONE
10 CONTINUE
A11_Re = DFLOAT(N1+N2)
ACOL(1) = DCMPLX(A11_Re,R_ZERO)
CALL FMSCOL (1, ACOL, LU_A11)
C
C Write remaining columns:
DO 20 I = 2,N1
ACOL(I) = C_ZERO
20 CONTINUE
ACOL(1) = C_MONE
DO 30 I = 2,N1
ACOL(I) = C_ONE
CALL FMSCOL (I, ACOL, LU_A11)
ACOL(I) = C_ZERO
30 CONTINUE
C
C End FMSCOL
CALL FMSCOL (N1+1, ACOL, LU_A11)
CALL FMSPOP (MYNAME)
RETURN
2000 FORMAT (' Populating [A11]:')
END
C=======================================================================
SUBROUTINE A12GEN (ACOL, N1, N2, LU_A12)
C=======================================================================
COMPLEX*16 ACOL(N1)
INTEGER N1, N2
INTEGER LU_A12(25)
CHARACTER*6 MYNAME
PARAMETER (MYNAME='A12GEN')
COMPLEX*16 C_MONE, C_ZERO
DATA C_MONE/(-1.0D0,0.0D0)/
DATA C_ZERO/( 0.0D0,0.0D0)/
C
CALL FMSPSH (MYNAME)
WRITE(6,2000)
C
DO 10 I=2,N1
ACOL(I) = C_ZERO
10 CONTINUE
ACOL(1) = C_MONE
C
C Initialize FMSPUT:
CALL FMSPUT (LU_A12, 0, 0, 0, 0, ACOL, 0)
C
C Write out columns of [A12]:
DO 20 I = 1,N2
CALL FMSPUT (LU_A12, N1, 1, 1, I, ACOL, N1)
20 CONTINUE
C
C End FMSPUT:
CALL FMSPUT (LU_A12, 0, 0, N1+1, 0, ACOL, 0)
C
C Set status:
LU_A12(7) = 0
CALL FMSPOP (MYNAME)
RETURN
2000 FORMAT (' Populating [A12]:')
END
C=======================================================================
SUBROUTINE A21GEN (AROW, N1, N2, LU_A21)
C=======================================================================
COMPLEX*16 AROW(N1)
INTEGER N1, N2
INTEGER LU_A21(25)
CHARACTER*6 MYNAME
PARAMETER (MYNAME='A21GEN')
COMPLEX*16 C_MONE, C_ZERO
DATA C_MONE/(-1.0D0,0.0D0)/
DATA C_ZERO/( 0.0D0,0.0D0)/
C
CALL FMSPSH (MYNAME)
WRITE(6,2000)
C
DO 10 I=2,N1
AROW(I) = C_ZERO
10 CONTINUE
AROW(1) = C_MONE
C
C Initialize FMSPUT:
CALL FMSPUT (LU_A21, 0, 0, 0, 0, AROW, 0)
C
C Write out rows of [A21], which is stored transposed:
DO 20 I = 1,N2
CALL FMSPUT (LU_A21, N1, 1, 1, I, AROW, N1)
20 CONTINUE
C
C End FMSPUT:
CALL FMSPUT (LU_A21, 0, 0, N1+1, 0, AROW, 0)
C
C Set status:
LU_A21(7) = 0
CALL FMSPOP (MYNAME)
RETURN
2000 FORMAT (' Populating [A21]:')
END
C=======================================================================
SUBROUTINE A22GEN (ACOL, N2, LU_A22)
C=======================================================================
COMPLEX*16 ACOL(N2)
INTEGER N2
INTEGER LU_A22(25)
CHARACTER*6 MYNAME
PARAMETER (MYNAME='A22GEN')
COMPLEX*16 C_ZERO, C_ONE
DATA C_ZERO/( 0.0D0,0.0D0)/
DATA C_ONE /( 1.0D0,0.0D0)/
C
CALL FMSPSH (MYNAME)
WRITE(6,2000)
C
C Initialize FMSCOL
CALL FMSCOL (0, ACOL, LU_A22)
C
DO 10 I = 1,N2
ACOL(I) = C_ZERO
10 CONTINUE
C
DO 20 I = 1,N2
ACOL(I) = C_ONE
CALL FMSCOL (I, ACOL, LU_A22)
ACOL(I) = C_ZERO
20 CONTINUE
C
C End FMSCOL
CALL FMSCOL (N2+1, ACOL, LU_A22)
C
C Set status:
LU_A22(7) = 0
CALL FMSPOP (MYNAME)
RETURN
2000 FORMAT (' Populating [A22]:')
END
C=======================================================================
SUBROUTINE B1GEN (RHS, N1, NRHS, LU_B1)
C=======================================================================
COMPLEX*16 RHS(N1)
INTEGER N1, NRHS
INTEGER LU_B1(25)
CHARACTER*5 MYNAME
PARAMETER (MYNAME='B1GEN')
COMPLEX*16 C_ZERO, C_ONE
DATA C_ZERO/(0.0D0,0.0D0)/
DATA C_ONE /(1.0D0,0.0D0)/
C
CALL FMSPSH (MYNAME)
WRITE(6,2000)
C
DO 10 I=2,N1
RHS(I) = C_ZERO
10 CONTINUE
RHS(1) = C_ONE
C
C Initialize FMSPUT:
CALL FMSPUT (LU_B1, 0, 0, 0, 0, RHS, 0)
C
C Write out columns of {B1}:
DO 20 I = 1,NRHS
CALL FMSPUT (LU_B1, N1, 1, 1, I, RHS, N1)
20 CONTINUE
C
C End FMSPUT:
CALL FMSPUT (LU_B1, 0, 0, N1+1, 0, RHS, 0)
C
C Set status:
LU_B1(7) = 0
CALL FMSPOP (MYNAME)
RETURN
2000 FORMAT (' Populating {B1}:')
END
C=======================================================================
SUBROUTINE B2GEN (RHS, N2, NRHS, LU_B2)
C=======================================================================
COMPLEX*16 RHS(N2)
INTEGER N2, NRHS
INTEGER LU_B2(25)
CHARACTER*5 MYNAME
PARAMETER (MYNAME='B2GEN')
COMPLEX*16 C_ZERO, C_ONE
DATA C_ZERO/(0.0D0,0.0D0)/
DATA C_ONE /(1.0D0,0.0D0)/
C
CALL FMSPSH (MYNAME)
WRITE(6,2000)
C
DO 10 I=1,N2
RHS(I) = C_ZERO
10 CONTINUE
C
C Initialize FMSPUT:
CALL FMSPUT (LU_B2, 0, 0, 0, 0, RHS, 0)
C
C Write out columns of {B2}:
DO 20 I = 1,NRHS
CALL FMSPUT (LU_B2, N2, 1, 1, I, RHS, N2)
20 CONTINUE
C
C End FMSPUT:
CALL FMSPUT (LU_B2, 0, 0, N2+1, 0, RHS, 0)
C
C Set status:
LU_B2(7) = 0
CALL FMSPOP (MYNAME)
RETURN
2000 FORMAT (' Populating {B2}:')
END
C=======================================================================
SUBROUTINE RESOLV
1 (LU_A11, LU_F11,
2 LU_A12, LU_F12,
2 LU_A21,
4 LU_A22, LU_F22,
3 LU_B1, LU_S1, LU_X1,
4 LU_B2, LU_S2, LU_X2,
5 NEWA11, NEWA12,
6 NEWA21, NEWA22,
7 NEWB1, NEWB2,
6 ICK_X2, E2MAX)
C=======================================================================
C
C DESCRIPTION:
C This subroutine solves the linear system
C
C +- + -+ +- -+ +- -+
C | A11 | A12 | | X1 | | B1 |
C +-----+-----+ +----+ = +----+
C | A21 | A22 | | X2 | | B2 |
C +- + -+ +- -+ +- -+
C
C when all or only part of the matrix terms change.
C
C
C On input, the following matrix terms are specified:
C
C +- + -+ +- -+ +- -+
C | A11 | A12 | | | | B1 |
C +-----+-----+ +----+ = +----+
C | A21 | A22 | | | | B2 |
C +- + -+ +- -+ +- -+
C
C
C The following intermediate work arrays are used:
C
C +- + -+ +- -+ +- -+
C | | | | S1 | | |
C +-----+-----+ +----+ = +----+
C | | | | S2 | | |
C +- + -+ +- -+ +- -+
C
C
C On output, the following terms are returned:
C
C +- + -+ +- -+ +- -+
C | F11 | F12 | | X1 | | |
C +-----+-----+ +----+ = +----+
C | | F22 | | X2 | | |
C +- + -+ +- -+ +- -+
C
C You may overlay any of the corresponding input, intermediate
C or output terms to save storage by specifying the same actual
C parameter on the subroutine call.
C
C Arrays containing input data may share storage with
C intermediate and output arrays if these conditions are met:
C 1) If any value changes in the input array, all values in that
C input array are redefined.
C 2) In defining the new values for an array, the previous
C values are not required. (They are no longer available
C because they have been overwritten)
C If these conditions are true, then the following arrays can
C share the same storage:
C [A11] & [F11],
C [A12] & [F12],
C [A22] & [F22],
C {B1} & {S1},
C {B2} & {S2}.
C
C Ths storage for {S1} and {S2} can be shared with {X1} and {X2}.
C No additional computation is required if the following
C conditions are met:
C {S1} & {X1} : if [A11] or {B1} change,
C {S2} & {X2} : if [A11], [A12], {B1} or {B2} change.
C If these conditions are not met and the storage is shared, then
C additional computation will be performed.
C
C You may also overlay {B1},{S1},{X1} or {B2},{S2},{X2} if
C all the above conditions are met.
C
C FORMAL PARAMETERS:
C (R ) LU_A11(25) = INTEGER ARRAY
C FMS matrix file attributes for [A11]
C
C (R ) LU_F11(25) = INTEGER ARRAY
C FMS matrix file attributes for [F11]
C May be the same as LU_A11.
C
C (R ) LU_A12(25) = INTEGER ARRAY
C FMS vector file attributes for [A12]
C
C (R ) LU_F12(25) = INTEGER ARRAY
C FMS vector file attributes for [F12]
C May be the same as LU_A12.
C
C (R ) LU_A21(25) = INTEGER ARRAY
C FMS vector file attributes for [A21]
C Data on this file is not changed.
C
C (R ) LU_A22(25) = INTEGER ARRAY
C FMS mtarix file attributes for [A22]
C
C (R ) LU_F22(25) = INTEGER ARRAY
C FMS matrix file attributes for [F22]
C
C (R ) LU_B1(25) = INTEGER ARRAY
C FMS vector file attributes for {B1}
C
C (R ) LU_S1(25) = INTEGER ARRAY
C FMS vector file attributes for {S1}.
C This is an intermediate work array.
C May be the same value as LU_B1.
C
C (R ) LU_X1(25) = INTEGER ARRAY
C FMS vector file attributes for {X1}.
C May be the same as LU_S1.
C
C (R ) LU_B2(25) = INTEGER ARRAY
C FMS vector file attributes for {B2}
C
C (R ) LU_S2(25) = INTEGER ARRAY
C FMS vector file attributes for {S2}.
C This is an intermediate work array.
C May be the same value as LU_B2.
C
C (R ) LU_X2(25) = INTEGER ARRAY
C FMS vector file attributes for {X2}.
C May be the same as LU_B2.
C
C (R ) NEWA11 = LOGICAL
C Indicates that [A11] has changed.
C
C (R ) NEWA12 = LOGICAL
C Indicated that [A12] has changed.
C
C (R ) NEWA21 = LOGICAL
C Indicated that [A21] has changed.
C
C (R ) NEWA22 = LOGICAL
C Indicated that [A22] has changed.
C
C (R ) NEWB1 = Logical
C Indicates that {B1} has changed.
C
C (R ) NEWB2 = Logical
C Indicates that {B2} has changed.
C
C (R ) ICK_X2 = Logical
C Check the values in {X2}.
C This requires access to [A22] and {B2}.
C Requires the following:
C LU_S2 or LU_X2 do not overlay LU_B2
C
C (R ) E2MAX = REAL*8
C Maximum error in {X2}.
C
C FMSCOM PARAMETERS:
C (RW) IACCOM = Multiply accumulate flag
C (R ) LUPR = FORTRAN unit for printing
C
C CALLED FROM:
C EXAMPLE_14
C
C SUBPROGRAMS USED:
C CNDF
C CNDS
C FMSMM
C
C ERROR CONDITIONS:
C None.
C
C HISTORY:
C
C NOTES:
C
C Expanding the partitioned system gives:
C
C [A11]{X1} + [A12]{X2} = {B1} (1)
C [A21]{X1} + [A22]{X2} = {B2} (2)
C
C Premultiplying (1) by inv[A11] gives
C
C {X1} = {S1} - [F12]{X2} (3)
C
C where
C {S1} = inv[A11]{B1} (4)
C [F12] = inv[A11][A12] (5)
C
C Substituting (3) into (2) gives:
C
C [F22]{X2} = {S2} (6)
C
C where
C [F22] = [A22] - [A21][F12] (7)
C {S2} = {B2} - [A21]{S1} (8)
C
C
C The error in {X2} may be computed from:
C
C E2MAX = | {B2} - [A21]{X1} - [A22]{X2} |
C
C The algorithm consists of the following 8 steps:
C
C Solve large system:
C 1. Factor [A11] into [F11]=[L11][U11]
C 2. Solve [F11][F12] = [A12] (eq. 5)
C 3. Solve [F11]{S1} = {B1} (eq. 4)
C
C Project to small system:
C 4. Multiply [F22'] = [A22] - [A21][F12] (eq. 7)
C 5. Factor [F22'] into [F22]=[L22][U22]
C 6. Multiply {S2} = {B2} - [A21]{S1} (eq. 8)
C
C Solve small system:
C 7. Solve [F22]{X2} = {S2} (eq. 6)
C
C Project back to large system:
C 8. Multiply {X1} = {S1} - [F12]*{X2} (eq. 3)
C
C Check {X2}: (optional)
C 9. Multiply {ERR} = {B2} - [A21]{X1} - [A22]{X2}
C Then compute the maximum absolute value in {ERR}.
C
C Summary of arrays used:
C
C STEP A11 F11 A12 F12 A21 A22 F22 B1 S1 X1 B2 S2 X2
C ==== ======= ======= === ======= =========== ===========
C 1 R-->W
C 2 R R-->W
C 3 R R-->W
C 4 R R R-->W
C 5 RW
C 6 R R R-->W
C 7 R R-->W
C 8 R R-->W R
C
C Steps required:
C
C STEP1 = (NEWA11 )
C STEP2 = (NEWA12 | STEP1 )
C STEP3 = (NEWB1 | STEP1 )
C STEP4 = (NEWA22 | NEWA21 | STEP2 )
C STEP5 = (STEP4 )
C STEP6 = (NEWB2 | NEWA21 | STEP3 )
C STEP7 = (STEP5 | STEP6 )
C STEP8 = (STEP3 | STEP2 | STEP7 )
C
C or, substituting,
C
C STEP1 = (NEWA11 )
C STEP2 = (NEWA11 | NEWA12 )
C STEP3 = (NEWA11 | NEWB1 )
C STEP4 = (NEWA11 | NEWA12 | NEWA21 | NEWA22 )
C STEP5 = (NEWA11 | NEWA12 | NEWA21 | NEWA22 )
C STEP6 = (NEWA11 | NEWA21 | NEWB1 | NEWB2 )
C STEP7 = (NEWA11 | NEWA12 | NEWA21 | NEWA22 | NEWB1 | NEWB2 )
C STEP8 = (NEWA11 | NEWA12 | NEWA21 | NEWA22 | NEWB1 | NEWB2 )
C
C Also if the intermediate vector storage is overlaid on the
C storage for the solution vectors, it is assumed that it no
C longer contains the intermediate arrays. Therefore,
C STEP3 = (STEP3 | (LU_S1 = LU_X1) )
C STEP6 = (STEP6 | (LU_S2 = LU_X2) )
C-----------------------------------------------------------------------
C FORMAL PARAMETERS
C-----------------------------------------------------------------------
INTEGER LU_A11(25)
INTEGER LU_F11(25)
C
INTEGER LU_A12(25)
INTEGER LU_F12(25)
C
INTEGER LU_A21(25)
C
INTEGER LU_A22(25)
C
INTEGER LU_F22(25)
C
INTEGER LU_B1 (25)
INTEGER LU_S1 (25)
INTEGER LU_X1 (25)
C
INTEGER LU_B2 (25)
INTEGER LU_S2 (25)
INTEGER LU_X2 (25)
C
LOGICAL NEWA11
LOGICAL NEWA12
LOGICAL NEWA21
LOGICAL NEWA22
LOGICAL NEWB1
LOGICAL NEWB2
LOGICAL ICK_X2
REAL*8 E2MAX
C-----------------------------------------------------------------------
C LOCAL VARIABLES
C-----------------------------------------------------------------------
C
CHARACTER*6 MYNAME
PARAMETER (MYNAME='RESOLV')
C
C FMS Memory management:
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
COMMON /MYFMS/ L_A22, LDA22
C
C Data type = complex:
INTEGER IDTYPE
PARAMETER (IDTYPE = 2)
C
C Skip operations during solving (no):
INTEGER ISKIP
PARAMETER (ISKIP=0)
C
C Constants:
COMPLEX*16 C_ONE(1)
DATA C_ONE /(1.0D0,0.0D0)/
REAL*8 R_ZERO
DATA R_ZERO / 0.0D0/
REAL*8 R_ONE
DATA R_ONE / 1.0D0/
REAL*8 R_THREE
DATA R_THREE/ 3.0D0/
REAL*8 R_EIGHT
DATA R_EIGHT/ 8.0D0/
C
C FMS Parameters:
INTEGER LUPR
C
C Matrix dimensions:
INTEGER N1
INTEGER N2
INTEGER NRHS
REAL*8 RN1
REAL*8 RN2
REAL*8 RNRHS
C
C Temporary file:
INTEGER LU_ERR(25)
INTEGER LUA0 (25)
DATA LUA0(1)/0/
INTEGER LUX0 (25)
INTEGER LUS0 (25)
C
C CPU, Wall and I/O time:
REAL*8 TCPU(10)
REAL*8 TWALL(10)
REAL*8 TIO(10)
REAL*8 T1, T2, T3
C
C Steps to perform:
LOGICAL STEP1
LOGICAL STEP2
LOGICAL STEP3
LOGICAL STEP4
LOGICAL STEP5
LOGICAL STEP6
LOGICAL STEP7
LOGICAL STEP8
LOGICAL STEP9
C
C Floating point operations:
REAL*8 OPS(9), TOPS
C-----------------------------------------------------------------------
C FORTRAN PROCEDURE
C-----------------------------------------------------------------------
CALL FMSPSH (MYNAME)
C
C Get the FORTRAN unit for printng:
CALL FMSIGT ('LUPR', LUPR)
C
C Get the memory pointer for arrays IMD, RMD and CMD:
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
C
C Determine the size of the matrices:
N1 = LU_A11(8)
N2 = LU_A22(8)
NRHS = LU_B1 (6)
NUMEQ = N1 + N2
RN1 = DFLOAT(N1)
RN2 = DFLOAT(N2)
RNRHS = DFLOAT(NRHS)
WRITE(LUPR,2000) NUMEQ, N1, N2, NRHS
C
C Determine which steps to do:
STEP1 = (NEWA11 )
STEP2 = (NEWA12 .OR. STEP1 )
STEP3 = (NEWB1 .OR. STEP1 )
STEP4 = (NEWA22 .OR. NEWA21 .OR. STEP2 )
STEP5 = (STEP4 )
STEP6 = (NEWB2 .OR. NEWA21 .OR. STEP3 )
STEP7 = (STEP5 .OR. STEP6 )
STEP8 = (STEP2 .OR. STEP3 .OR. STEP7 )
STEP9 = (ICK_X2 )
C
C Check overlaid storage:
STEP3 = (STEP3 .OR. (LU_S1(1) .EQ. LU_X1(1)) )
STEP6 = (STEP6 .OR. (LU_S2(2) .EQ. LU_X2(2)) )
IF( (LU_B2(1) .EQ. LU_S2(1) ) .OR.
1 (LU_B2(1) .EQ. LU_X2(1) ) ) STEP9 = .FALSE.
C
C Initialize operation counts for each step:
OPS(1) = R_ZERO
OPS(2) = R_ZERO
OPS(3) = R_ZERO
OPS(4) = R_ZERO
OPS(5) = R_ZERO
OPS(6) = R_ZERO
OPS(7) = R_ZERO
OPS(8) = R_ZERO
OPS(9) = R_ZERO
C
C (1) Factor [A11] into [L11][U11]:
CALL FMSTIM (TCPU(1), TWALL(1), TIO(1))
IF(STEP1) THEN
WRITE(LUPR,2001)
OPS(1) = R_EIGHT*RN1*RN1*RN1/R_THREE
CALL CNDAF (LU_A11, C_ONE, 1, LUS0, 0, LUA0, LU_F11,
1 LUX0, LUX0, 0)
END IF
C
C (2) Solve [F11][F12] = [A12]:
CALL FMSTIM (TCPU(2), TWALL(2), TIO(2))
IF(STEP2) THEN
WRITE(LUPR,2002)
OPS(2) = R_EIGHT*RN1*RN1*RN2
CALL CNDS (LU_F11, LU_A12, LU_F12, N2, 0)
END IF
C
C (3) Solve [F11]{S1} = {B1}:
CALL FMSTIM (TCPU(3), TWALL(3), TIO(3))
IF(STEP3) THEN
WRITE(LUPR,2003)
OPS(3) = R_EIGHT*RN1*RN1*RNRHS
CALL CNDS (LU_F11, LU_B1, LU_S1, NRHS, 0)
END IF
C
C (4) Multiply [F22'] = [A22] - [A21][F12]
CALL FMSTIM (TCPU(4), TWALL(4), TIO(4))
IF(STEP4) THEN
WRITE(LUPR,2004)
OPS(4) = R_EIGHT*RN1*RN2*RN2
CALL FMSMM (LU_F22, LU_A22, -1, 'T', LU_A21, LU_F12)
END IF
C
C (5) Factor [F22'] into [F22]:
CALL FMSTIM (TCPU(5), TWALL(5), TIO(5))
IF(STEP5) THEN
WRITE(LUPR,2005)
OPS(5) = R_EIGHT*RN2*RN2*RN2/R_THREE
CALL CNDAF (LU_F22, C_ONE, 1, LUS0, 0, LUA0, LU_F22,
1 LUX0, LUX0, 0)
END IF
C
C (6) Multiply {S2} = {B2} - [A21]{S1}
CALL FMSTIM (TCPU(6), TWALL(6), TIO(6))
IF(STEP6) THEN
WRITE(LUPR,2006)
OPS(6) = R_EIGHT*RN2*RN1*RNRHS
CALL FMSMM (LU_S2, LU_B2, -1, 'T', LU_A21, LU_S1)
END IF
C
C (7) Solve [F22]{X2} = {S2}:
CALL FMSTIM (TCPU(7), TWALL(7), TIO(7))
IF(STEP7) THEN
WRITE(LUPR,2007)
OPS(7) = R_EIGHT*RN2*RN2*RNRHS
CALL CNDS (LU_F22, LU_S2, LU_X2, NRHS, 0)
END IF
C
C (8) Multiply {X1} = {S1} - [F12]{X2}
CALL FMSTIM (TCPU(8), TWALL(8), TIO(8))
IF(STEP8) THEN
WRITE(LUPR,2008)
OPS(8) = R_EIGHT*RN1*RN2*RNRHS
CALL FMSMM (LU_X1, LU_S1, -1, 'N', LU_F12, LU_X2)
END IF
C
C (9) Check {X2}: {ERR} = {B2} - [A21]{X1} - [A22]{X2}
CALL FMSTIM (TCPU(9), TWALL(9), TIO(9))
IF(STEP9) THEN
WRITE(LUPR,2009)
OPS(9) = R_EIGHT*(RN1+RN2)*RN2*RNRHS
C
C Open a vector file to store the error array ERR(N2,NRHS):
CALL FMSOV (N2, IDTYPE, NRHS, 'LU_ERR' , LU_ERR)
C
C Multiply {ERR} = {B2} - [A21]{X1}
CALL FMSMM (LU_ERR, LU_B2, -1, 'T', LU_A21, LU_X1)
C
C Multiply {ERR} = {ERR} - [A22]{X2}:
CALL FMSMM (LU_ERR, LU_ERR, -1, 'N', LU_A22, LU_X2)
C
C Read LU_ERR and check the error:
C
C Allocate storage for one vector:
CALL FMSCMG (CMD, L_ERR, N2)
C
C Initialize FMSGET:
CALL FMSGET (LU_ERR, 0, 0, 0, 0, CMD, 0)
C
C Compute the maximum error:
E2MAX = R_ZERO
DO 95 N = 1,NRHS
C Read in next vector:
CALL FMSGET (LU_ERR, N2, 1, 1, N, CMD(L_ERR), N2)
L = L_ERR
DO 90 I = 1,N2
T1 = ABS(CMD(L))
IF(T1 .GT. E2MAX) E2MAX = T1
L = L + 1
90 CONTINUE
95 CONTINUE
C
C End FMSGET:
CALL FMSGET (LU_ERR, 0, 0, N2+1, 1, CMD, 0)
C
C Close the vector file LU_ERR:
CALL FMSCV(LU_ERR)
C
C Delete the storage for {ERR}:
CALL FMSCMR (CMD, L_ERR, N2)
END IF
C
C Print timings:
WRITE(LUPR,2010)
CALL FMSTIM (TCPU(10), TWALL(10), TIO(10))
TOPS = R_ZERO
DO 100 ISTEP = 1,9
IF(OPS(ISTEP) .EQ. R_ZERO) THEN
C
C Operation skipped:
WRITE(LUPR,2012) ISTEP
ELSE
TOPS = TOPS + OPS(ISTEP)
C
C CPU time and Mflops:
T1 = TCPU(ISTEP+1) - TCPU(ISTEP)
IF(T1 .GT. R_ZERO) THEN
MF1 = ((1.0D-6)*OPS(ISTEP)/T1)
ELSE
MF1 = 0
END IF
C
C Wall time, Mflops and ratio:
T2 = TWALL(ISTEP+1) - TWALL(ISTEP)
T3 = TIO(ISTEP+1) - TIO(ISTEP)
IF((T2 .GT. R_ZERO) .AND.
1 (T2 .GT. T3 ) ) THEN
MF2 = ((1.0D-6)*OPS(ISTEP)/T2)
RATIO = T1/(T2-T3)
ELSE
MF2 = 0
RATIO = R_ONE
END IF
WRITE(LUPR,2011) ISTEP, T1, T2, T3, MF1, MF2, RATIO
END IF
100 CONTINUE
T1 = TCPU(10) - TCPU(1)
T2 = TWALL(10) - TWALL(1)
T3 = TIO(10) - TIO(1)
IF(T1 .GT. R_ZERO) THEN
MF1 = ((1.0D-6)*TOPS/T1)
ELSE
MF1 = 0
END IF
IF(T2 .GT. R_ZERO) THEN
MF2 = ((1.0D-6)*TOPS/T2)
RATIO = T1/(T2-T3)
ELSE
MF2 = 0
RATIO = R_ONE
END IF
WRITE(LUPR,2013) T1, T2, T3, MF1, MF2, RATIO
CALL FMSPOP (MYNAME)
RETURN
2000 FORMAT (/
1 ' Start of Subroutine RESOLV',/
2 ' Rank of full matrix.....=',I6/
3 ' Rank of submatrix #1....=',I6/
4 ' Rank of submatrix #2....=',I6/
5 ' # of right-hand-sides...=',I6/)
2001 FORMAT (/
1 ' STEP 1. Factor [A11] into [F11]'/
2 ' =================================')
2002 FORMAT (/
1 ' STEP 2. Solve [F11][F12] = [A12]'/
2 ' ===================================')
2003 FORMAT (/
1 ' STEP 3. Solve [F11]{S1} = {B1}'/
2 ' ==================================')
2004 FORMAT (/
1 ' STEP 4. Multiply [F22] = [A22] - [A21][F12]'/
2 ' ===========================================')
2005 FORMAT (/
1 ' STEP 5. Factor [F22]'/
2 ' ======================')
2006 FORMAT (/
1 ' STEP 6. Multiply {S2} = {B2} - [A21]{S1}'/
2 ' ==========================================')
2007 FORMAT (/
1 ' STEP 7. Solve [F22]{X2} = {S2}'/
2 ' =================================')
2008 FORMAT (/
1 ' STEP 8. Multiply {X1} = {S1} - [A12]*{X2}'/
2 ' =========================================')
2009 FORMAT (/
1 ' STEP 9. Checking values in {X2}'/
2 ' =========================================')
2010 FORMAT (/
1 ' Summary of Subroutine RESOLV:'//
2 ' <----------TIMES (Sec.)----------->'
2 ' <--MFLOPS-->'/
2 ' STEP CPU WALL I/O CPU WALL',
2 ' RATIO'/
3 ' ===== ========= ========= ========= ==== ====',
3 ' =====')
2011 FORMAT (I6,3F13.3,2I8,F9.3)
2012 FORMAT (I6,6X,'Skipped')
2013 FORMAT (' TOTAL',3F13.3,2I8,F9.3)
END