+- + -+ +- -+ +- -+ | A11 | A12 | | X1 | | B1 | +-----+-----+ +----+ = +----+ | A21 | A22 | | X2 | | B2 | +- + -+ +- -+ +- -+The system of equations is partitioned as shown above. On the first solution, all terms are computed. On following solutions, only those terms which must be recomputed due to changes in the matrix elements are reevaluated. The main work is performed by subroutine RESOLV, which is designed to be cut and pasted into your application.
The following input is used to test subroutine RESOLV:
- Number of equations in [A11]
- Number of equations in [A22]
- Number of solution vectors
- Any FMS Parameter
- Do you want another solution?
- Do you want a new A11?
- Do you want a new A12?
- Do you want a new A21?
- Do you want a new A22?
- Do you want a new X1?
- Do you want a new X2?
If you answer yes to any of these questions, the problem is solved again, recomputing only those terms that are necessary.
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 FMSCV (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 ('N', -1, LU_F12, LU_X2, LU_S1, LU_X1)
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