+- + -+ +- -+ +- -+ | A11 | A12 | | X1 | | B1 | +-----+-----+ +----+ = +----+ | A21 | A22 | | X2 | | B2 | +- + -+ +- -+ +- -+The system of equations is partitioned as shown above. The diagonal blocks [A11] and [A22] are factored but the off-diagonal blocks [A12] and [A21] are not. A sequence of iterative solutions is performed using Jacobi or Gauss-Seidel iteration to obtain approximate solutions {X1} and {X2}.
This approach can save computer time over the exact solution under the following conditions:
- The matrix values in the diagonal blocks [A11] and [A22] are large compared to the values in the off-diagonal blocks [A12] and [A21].
- An approximate solution to {X1} and {X2} is acceptable.
- The number of vectors, times the number of iterations, is less than the dimension of [A22].
A variation of the Jacobi method is also presented where the solution values are the change in {DX1} and {DX2}. This approach is numerically equivalent to the standard Jacobi method that solves directly for the solution {X1} and {X2}.
Subroutine J_GS provides both Jacobi and Gauss-Seidel solutions. Subroutine Jacobi2 uses the incremental Jacobi approach.
The following input is used to test these subroutines:
- Number of equations in [A11]
- Number of equations in [A22]
- Number of solution vectors
- Maximum number of iterations to perform
- Maximum solution error to terminate iterations
- Solution method
- 1=Jacobi
- 2=Gauss-Seidel
- 3=Incremental Jacobi
- Vector norm used to test for convergence
- 0=Infinity norm (maximum absolute value)
- 1=Sum of the absolute values
- 2=Square root of the sum of the values squared
- If FMS output is required during the iterations.
- Any FMS Parameter
- Do you want another solution?
- Do you want a new A11?
- Do you want a new A22?
If you answer yes to any of these questions, the problem is solved again, recomputing only those terms that are necessary.
When all solutions of the given problem are complete, you are asked if you want to solve another problem, in which any of the input may change.
C T E S T J A C O B I A N D G A U S S - S E I D E L
C
C This program tests subroutines J_GS and Jacobi2, which
C illustrate how the FMS subroutines may be used for a block
C iterative solution.
C
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_17')
C
C FMS matrix and vector file attributes:
INTEGER LU_A11(25)
INTEGER LU_A12(25)
INTEGER LU_A21(25)
INTEGER LU_A22(25)
INTEGER LU_B1 (25), LU_X1(25)
INTEGER LU_B2 (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
REAL*8 ASK_R
C
C Solution methods:
INTEGER METHOD_JACOBI
PARAMETER (METHOD_JACOBI =1)
INTEGER METHOD_GAUSS_SEIDEL
PARAMETER (METHOD_GAUSS_SEIDEL =2)
INTEGER METHOD_JACOBI_INCREMENT
PARAMETER (METHOD_JACOBI_INCREMENT =3)
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 [A22] matrix:
LOGICAL NEWA22
C Type of solution method:
INTEGER METHOD
C Print FMS output during interation:
LOGICAL NOFMSP
C Maximum number of iterations:
INTEGER MAXITR
C Maximum iteration error:
REAL*8 ERRMAX
C Error on current term:
REAL*8 EI
C Maximum overall error:
REAL*8 ERROR
C CMD Pointer to scratch vector:
INTEGER L_TEMP
C Profile vector for a full matrix:
INTEGER LOWEQ(1)
DATA LOWEQ/-1/
C
C CONSTANTS:
REAL*8 R_ZERO
DATA R_ZERO/0.0D0/
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:
10 CONTINUE
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.....')
MAXITR = ASK_I('Enter the maximum number of iterations...')
ERRMAX = ASK_R('Enter the maximum error..................')
WRITE(6,*) 'The following solution methods are available:'
WRITE(6,*) ' 1 = Jacobi,'
WRITE(6,*) ' 2 = Gauss-Seidel,'
WRITE(6,*) ' 3 = Incremental Jacobi.'
METHOD = ASK_I('Enter solution method (1 | 2 | 3 ).......')
WRITE(6,*) 'The following vector norms are available:'
WRITE(6,*) ' 0 = Infinity norm (MAX(ABS(X(i)))'
WRITE(6,*) ' 1 = SUM( ABS( X(i) ) )'
WRITE(6,*) ' 2 =( SUM( ABS( X(i) )**2) )**(1/2)'
INORM = ASK_I('Enter the norm number (0 | 1 | 2)........')
NOFMSP = ASK ('Print FMS output during iterations?')
NOFMSP = (.NOT.NOFMSP)
WRITE (6,*) 'You may now alter any FMS parameter.'
WRITE (6,*) 'When you are finished, type the letters RETURN'
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 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_X2' , LU_X2 )
C
C (3) Write data to FMS files:
NEWA11 = .TRUE.
NEWA22 = .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]:
CALL A12GEN (CMD(L_TEMP), N1, N2, LU_A12)
C
C Generate [A21]:
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}:
CALL B1GEN (CMD(L_TEMP), N1, NRHS, LU_B1 )
C
C Generate {B2}:
CALL B2GEN (CMD(L_TEMP), N2, NRHS, LU_B2 )
C
C (4) Perform matrix algebra:
IF( (METHOD .EQ. METHOD_JACOBI) .OR.
1 (METHOD .EQ. METHOD_GAUSS_SEIDEL) ) THEN
CALL J_GS
1 (LU_A11, LU_A11,
2 LU_A12,
3 LU_A21,
4 LU_A22, LU_A22,
5 LU_B1, LU_X1,
6 LU_B2, LU_X2,
7 NEWA11, NEWA22,
8 INORM, MAXITR, ERRMAX, METHOD, NOFMSP)
C
ELSE IF(METHOD .EQ. METHOD_JACOBI_INCREMENT) THEN
CALL JACOBI2
1 (LU_A11, LU_A11,
2 LU_A12,
3 LU_A21,
4 LU_A22, LU_A22,
5 LU_B1, LU_X1,
6 LU_B2, LU_X2,
7 NEWA11, NEWA22,
8 INORM, MAXITR, ERRMAX, NOFMSP)
END IF
C
C (5) Read data from FMS files:
C Check the answer:
ERROR = R_ZERO
C
C Check {X1}:
CALL FMSGET (LU_X1, 0, 0, 0, 0, CMD, 0)
DO IRHS = 1,NRHS
CALL FMSGET (LU_X1, N1, 1, 1, IRHS, CMD(L_TEMP), N1)
DO I = 0,(N1-1)
EI = ABS(CMD(L_TEMP+I) - C_ONE)
IF(EI .GT. ERROR) ERROR = EI
END DO
END DO
CALL FMSGET (LU_X1, 0, 0, N1+1, 1, CMD, 0)
C
C Check {X2}:
CALL FMSGET (LU_X2, 0, 0, 0, 0, CMD, 0)
DO IRHS = 1,NRHS
CALL FMSGET (LU_X2, N2, 1, 1, IRHS, CMD(L_TEMP), N2)
DO I = 0,(N2-1)
EI = ABS(CMD(L_TEMP+I) - C_ONE)
IF(EI .GT. ERROR) ERROR = EI
END DO
END DO
CALL FMSGET (LU_X2, 0, 0, N2+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 with this matrix?'))
1 GO TO 100
NEWA11 = ASK('Do you want a new A11?')
NEWA22 = ASK('Do you want a new A22?')
C
IF( NEWA11 .OR.
1 NEWA22) 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 FMSCV (LU_B1 )
CALL FMSCV (LU_X1 )
CALL FMSCV (LU_B2 )
CALL FMSCV (LU_X2 )
CALL FMSCMR (CMD, L_TEMP, NMAX)
IF(ASK('Do you want to solve another problem?')) GO TO 10
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=======================================================================
REAL*8 FUNCTION ASK_R(STRING)
C=======================================================================
CHARACTER* (*) STRING
WRITE(6,2000) STRING
READ (5,*) ASK_R
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 A11_Re, R_ZERO
COMPLEX*16 C_ZERO, C_ONE, C_MONE
DATA R_ZERO/0.0D0/
DATA C_ZERO/( 0.0D0,0.0D0)/
DATA C_MONE/(-1.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_ZERO, C_MONE
DATA C_ZERO/( 0.0D0,0.0D0)/
DATA C_MONE/(-1.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_ZERO, C_MONE
DATA C_ZERO/( 0.0D0,0.0D0)/
DATA C_MONE/(-1.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 J_GS
1 (LU_A11, LU_F11,
2 LU_A12,
2 LU_A21,
4 LU_A22, LU_F22,
3 LU_B1, LU_X1,
4 LU_B2, LU_X2,
5 NEWA11, NEWA22,
6 INORM, MAXITR, ERRMAX, METHOD, NOFMSP)
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 using block Jacobi or Gauss_Seidel iteration.
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 | | | | X1_OLD | | |
C +-----+-----+ +--------+ = +----+
C | | | | X2_OLD | | |
C +- + -+ +- -+ +- -+
C
C
C On output, the following terms are returned:
C
C +- + -+ +- -+ +- -+
C | F11 | | | X1 | | |
C +-----+-----+ +----+ = +----+
C | | F22 | | X2 | | |
C +- + -+ +- -+ +- -+
C
C You may overlay any of the corresponding input or output terms
C to save storage by specifying the same actual parameter on the
C subroutine call.
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 Data on this file is not changed.
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 matrix 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_X1(25) = INTEGER ARRAY
C FMS vector file attributes for {X1}.
C May be the same as LU_B1.
C
C (R ) LU_B2(25) = INTEGER ARRAY
C FMS vector file attributes for {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 ) NEWA22 = LOGICAL
C Indicates that [A22] has changed.
C
C (R ) INORM = Integer
C Norm for computing convergence:
C = 0, Infinity norm (max. abs. value)
C = 1, First norm
C = 2, Second norm
C
C (R ) MAXITR = Integer
C The maximum number of iterations to
C perform.
C
C (R ) ERRMAX = REAL*8
C = Maximum error for convergence
C
C (R ) METHOD = Integer
C = 1, Jacobi
C = 2, Gauss-Seidel
C
C (R ) NOFMSP = Logical
C Print FMS output in iteration loop.
C
C FMSCOM PARAMETERS:
C (RW) IACCOM = Multiply accumulate flag
C (R ) LUPR = FORTRAN unit for printing
C
C CALLED FROM:
C EXAMPLE_17
C
C SUBPROGRAMS USED:
C CNDF
C CNDS
C FMSMM
C
C ERROR CONDITIONS:
C None.
C
C HISTORY:
C
C NOTES:
C
C-----------------------------------------------------------------------
C FORMAL PARAMETERS
C-----------------------------------------------------------------------
INTEGER LU_A11(25)
INTEGER LU_F11(25)
C
INTEGER LU_A12(25)
C
INTEGER LU_A21(25)
C
INTEGER LU_A22(25)
INTEGER LU_F22(25)
C
INTEGER LU_B1 (25)
INTEGER LU_X1 (25)
C
INTEGER LU_B2 (25)
INTEGER LU_X2 (25)
C
LOGICAL NEWA11
LOGICAL NEWA22
INTEGER INORM
INTEGER MAXITR
REAL*8 ERRMAX
INTEGER METHOD
LOGICAL NOFMSP
C-----------------------------------------------------------------------
C LOCAL VARIABLES
C-----------------------------------------------------------------------
C
CHARACTER*4 MYNAME
PARAMETER (MYNAME='J_GS')
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
C Data type = complex:
INTEGER IDTYPE
PARAMETER (IDTYPE = 2)
C
C Skip operations during solving (no):
INTEGER ISKIP
PARAMETER (ISKIP=0)
C
C Number of solution steps:
INTEGER NSTEPS
PARAMETER (NSTEPS=11)
C
C Method of solution:
INTEGER METHOD_JACOBI
PARAMETER (METHOD_JACOBI=1)
INTEGER METHOD_GAUSS_SEIDEL
PARAMETER (METHOD_GAUSS_SEIDEL=2)
C
C Constants:
COMPLEX*16 C_ONE(1)
DATA C_ONE /( 1.0D0,0.0D0)/
COMPLEX*16 C_MONE
DATA C_MONE /(-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_TWO
DATA R_TWO /2.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
INTEGER IPRF
INTEGER IPRS
INTEGER IPRVV
C
C Matrix dimensions:
INTEGER N1
INTEGER N2
INTEGER NRHS
REAL*8 RN1
REAL*8 RN2
REAL*8 RNRHS
C
C Convergence norms:
REAL*8 R1
REAL*8 R2
C
C Temporary files:
INTEGER LU_X1_OLD(25)
INTEGER LU_X2_OLD(25)
INTEGER LUA0 (25)
DATA LUA0(1)/0/
INTEGER LUS0 (25)
INTEGER LUX0 (25)
C
C CPU, Wall and I/O time:
REAL*8 TCPU(NSTEPS)
REAL*8 T1, TCPU_T1 ,TCPU_T2
REAL*8 TWALL(NSTEPS)
REAL*8 T2, TWALL_T1 ,TWALL_T2
REAL*8 TIO(NSTEPS)
REAL*8 T3, TIO_T1 ,TIO_T2
C
C Floating point operations:
REAL*8 OPS(NSTEPS), TOPS
INTEGER NCALLS(NSTEPS)
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, MAXITR, ERRMAX
C
C Initialize operation counts and timers for each step:
DO I=1,NSTEPS
OPS (I) = R_ZERO
TCPU (I) = R_ZERO
TWALL (I) = R_ZERO
TIO (I) = R_ZERO
NCALLS(I) = 0
END DO
C
C Open temporary files for previous solution:
CALL FMSOV (N1, IDTYPE, NRHS, 'LU_X1_OLD' , LU_X1_OLD )
C
CALL FMSTIM (TCPU_T1, TWALL_T1, TIO_T1)
CALL TIMER (0, R_ZERO, NCALLS, OPS, TCPU, TWALL, TIO)
C
C Initialize {X1_OLD} and {X2_OLD}:
C
C (1) Factor [A11] into [F11]=[L11][U11]:
ISTEP = 1
IF(NEWA11) THEN
WRITE(LUPR,2001)
CALL CNDAF (LU_A11, C_ONE, 1, LUS0, 0, LUA0, LU_F11,
1 LUX0, LUX0, 0)
CALL TIMER (1, R_EIGHT*RN1*RN1*RN1/R_THREE,
1 NCALLS, OPS, TCPU, TWALL, TIO)
END IF
C
C (2) Solve [F11]{X1_OLD} = {B1}:
WRITE(LUPR,2002)
CALL CNDS (LU_F11, LU_B1, LU_X1_OLD, NRHS, 0)
CALL TIMER (2, R_EIGHT*RN1*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (3) Factor [A22] into [F22]=[L22][U22]:
IF(NEWA22) THEN
WRITE(LUPR,2003)
CALL CNDAF (LU_A22, C_ONE, 1, LUS0, 0, LUA0, LU_F22,
1 LUX0, LUX0, 0)
CALL TIMER (3, R_EIGHT*RN2*RN2*RN2/R_THREE,
1 NCALLS, OPS, TCPU, TWALL, TIO)
END IF
C
IF(METHOD .EQ. METHOD_JACOBI) THEN
C
C (4) Not performed:
C
C (5) Solve [F22]{X2_OLD} = {B2}:
WRITE(LUPR,2105)
CALL CNDS (LU_F22, LU_B2, LU_X2_OLD, NRHS, 0)
CALL TIMER (5, R_EIGHT*RN2*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
ELSE IF(METHOD .EQ. METHOD_GAUSS_SEIDEL) THEN
C
C (4) Multiply {X2_OLD} = {B2} - [A21]{X1_OLD}:
WRITE(LUPR,2204)
CALL FMSMM (LU_X2_OLD, LU_B2, -1, 'T', LU_A21, LU_X1_OLD)
LU_X1_OLD(7) = 0
CALL TIMER (4, R_EIGHT*RN2*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (5) Solve [F22]{X2_OLD} = {X2_OLD}:
WRITE(LUPR,2205)
CALL CNDS (LU_F22, LU_X2_OLD, LU_X2_OLD, NRHS, 0)
CALL TIMER (5, R_EIGHT*RN2*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
END IF
C
C Turn off FMS printing in loop, if required:
IF(NOFMSP) THEN
CALL FMSIGT ('IPRF, IPRF)
CALL FMSIST ('IPRF, 0)
CALL FMSIGT ('IPRS, IPRS)
CALL FMSIST ('IPRS, 0)
CALL FMSIGT ('IPRVV, IPRVV)
CALL FMSIST ('IPRVV, 0)
END IF
C
C Loop over iterations:
DO 100 NITER = 1,MAXITR
C
C (6) Multiply {X1} = {B1} - [A12]{X2_OLD}:
IF(.NOT.NOFMSP) WRITE(LUPR,2006)
CALL FMSMM (LU_X1, LU_B1, -1, 'N', LU_A12, LU_X2_OLD)
LU_X1(7) = 0
CALL TIMER (6, R_EIGHT*RN1*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (7) Solve [F11]{X1} = {X1}:
IF(.NOT.NOFMSP) WRITE(LUPR,2007)
CALL CNDS (LU_F11, LU_X1, LU_X1, NRHS, 0)
CALL TIMER (7, R_EIGHT*RN1*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
IF(METHOD .EQ. METHOD_JACOBI) THEN
C
C (8) Multiply {X2} = {B2} - [A21]{X1_OLD}:
IF(.NOT.NOFMSP) WRITE(LUPR,2108)
CALL FMSMM (LU_X2, LU_B2, -1, 'T', LU_A21, LU_X1_OLD)
C
ELSE IF(METHOD .EQ. METHOD_GAUSS_SEIDEL) THEN
C
C (8) Multiply {X2} = {B2} - [A21]{X1}:
IF(.NOT.NOFMSP) WRITE(LUPR,2208)
CALL FMSMM (LU_X2, LU_B2, -1, 'T', LU_A21, LU_X1)
C
END IF
LU_X2(7) = 0
CALL TIMER (8, R_EIGHT*RN2*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (9) Solve [F22]{X2} = {X2}:
IF(.NOT.NOFMSP) WRITE(LUPR,2009)
CALL CNDS (LU_F22, LU_X2, LU_X2, NRHS, 0)
CALL TIMER (9, R_EIGHT*RN2*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (10) R1=|| ( {X1} - {X1_OLD} )/ {X1} ||; {X1}->{X1_OLD}:
IF(.NOT.NOFMSP) WRITE(LUPR,2010)
CALL FMSVAN (LU_X1_OLD, -1, LU_X1, INORM, R1, 3)
IF(INORM .GT. 1) THEN
R1 = R1/RN1
ELSE IF(INORM .EQ. 2) THEN
R1 = R1/SQRT(RN1)
END IF
CALL TIMER (10, R_TWO*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (11) R2=|| ( {X2} - {X2_OLD} )/ {X2} ||; {X2}->{X2_OLD}:
IF(.NOT.NOFMSP) WRITE(LUPR,2011)
CALL FMSVAN (LU_X2_OLD, -1, LU_X2, INORM, R2, 3)
IF(INORM .GT. 1) THEN
R2 = R2/RN2
ELSE IF(INORM .EQ. 2) THEN
R2 = R2/SQRT(RN2)
END IF
CALL TIMER (11, R_TWO*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
WRITE(LUPR,2050) NITER, R1, R2
C
IF( ( R1 .LE. ERRMAX) .AND.
1 ( R2 .LE. ERRMAX) ) GO TO 101
100 CONTINUE
101 CONTINUE
C
C Close temporary files:
CALL FMSCV (LU_X1_OLD)
CALL FMSCV (LU_X2_OLD)
C
C Print timings:
WRITE(LUPR,2020)
CALL FMSTIM (TCPU_T2, TWALL_T2, TIO_T2)
TOPS = R_ZERO
DO 200 ISTEP = 1,NSTEPS
NCALL = NCALLS(ISTEP)
IF(NCALL .EQ. 0) THEN
C
C Operation skipped:
WRITE(LUPR,2022) ISTEP
ELSE
TOPS = TOPS + OPS(ISTEP)
C
C CPU time and Mflops:
T1 = 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)
T3 = 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,2021) ISTEP, NCALL, T1, T2, T3, MF1, MF2, RATIO
END IF
200 CONTINUE
T1 = TCPU_T2 - TCPU_T1
T2 = TWALL_T2 - TWALL_T1
T3 = TIO_T2 - TIO_T1
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
IF(NOFMSP) THEN
C Restore print settings:
CALL FMSIST ('IPRF, IPRF)
CALL FMSIST ('IPRS, IPRS)
CALL FMSIST ('IPRVV, IPRVV)
END IF
WRITE(LUPR,2023) T1, T2, T3, MF1, MF2, RATIO
CALL FMSPOP(MYNAME)
RETURN
2000 FORMAT (/
1 ' Start of Subroutine J_GS',/
2 ' Rank of full matrix............=',I8/
3 ' Rank of submatrix #1...........=',I8/
4 ' Rank of submatrix #2...........=',I8/
5 ' # of right-hand-sides..........=',I8/
6 ' Maximum number of iterations...=',I8/
7 ' Maximum error norm.............=',E8.1/)
2001 FORMAT (/
1 ' STEP 1. Factor [A11] into [F11]'/
2 ' ==================================')
2002 FORMAT (/
1 ' STEP 2. Solve [F11]{X1_OLD} = {B1}'/
2 ' =======================================')
2003 FORMAT (/
1 ' STEP 3. Factor [A22] into [F22]'/
2 ' ==================================')
2204 FORMAT (/
1 ' STEP 4. Multiply {X2_OLD} = {B2} - [A21]{X1_OLD}'/
2 ' =================================================')
2105 FORMAT (/
1 ' STEP 5. Solve [F22]{X2_OLD} = {B2}'/
2 ' ======================================')
2205 FORMAT (/
1 ' STEP 5. Solve [F22]{X2_OLD} = {X2_OLD}'/
2 ' ===========================================')
2006 FORMAT (/
1 ' STEP 6. Multiply {X1} = {B1} - [A12]{X2_OLD}'/
2 ' =============================================')
2007 FORMAT (/
1 ' STEP 7. Solve [F11]{X1} = {X1}'/
2 ' ==================================')
2108 FORMAT (/
1 ' STEP 8. Multiply {X2} = {B2} - [A21]{X1_OLD}'/
2 ' =============================================')
2208 FORMAT (/
1 ' STEP 8. Multiply {X2} = {B2} - [A21]{X1}'/
2 ' =========================================')
2009 FORMAT (/
1 ' STEP 9. Solve [F22]{X2} = {X2}'/
2 ' ==================================')
2010 FORMAT (/
1 ' STEP 10. |( {X1} - {X1_OLD} )/{X1}; {X1}->{X1_OLD}'/,
2 ' ==================================================')
2011 FORMAT (/
1 ' STEP 11. |( {X2} - {X2_OLD} )/{X2}; {X2}->{X2_OLD}'/,
2 ' ===================================================')
2020 FORMAT (/
1 ' Summary of Subroutine J_GS:'//
2 ' <----------TIMES (Sec.)----------->',
2 ' <--MFLOPS-->'/
2 ' STEP RUN CPU WALL I/O',
2 ' CPU WALL',
2 ' RATIO'/
3 ' ===== === ========= ========= =========',
3 ' ==== ====',
3 ' =====')
2021 FORMAT (2I6,3F13.3,2I8,F9.3)
2022 FORMAT (I6,12X,'Skipped')
2023 FORMAT (' TOTAL',6X,3F13.3,2I8,F9.3)
2050 FORMAT (
1 ' ITERATION=',I3,
2 ', {DX1}NORM=',E12.3,
3 ', {DX2}NORM=',E12.3)
END
C=======================================================================
SUBROUTINE JACOBI2
1 (LU_A11, LU_F11,
2 LU_A12,
2 LU_A21,
4 LU_A22, LU_F22,
3 LU_B1, LU_X1,
4 LU_B2, LU_X2,
5 NEWA11, NEWA22,
6 INORM, MAXITR, ERRMAX, NOFMSP)
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 using an incremental form of block Jacobi iteration.
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 | | | | DX1(2) | | |
C +-----+-----+ +--------+ = +----+
C | | | | DX2(2) | | |
C +- + -+ +- -+ +- -+
C
C
C On output, the following terms are returned:
C
C +- + -+ +- -+ +- -+
C | F11 | | | X1 | | |
C +-----+-----+ +----+ = +----+
C | | F22 | | X2 | | |
C +- + -+ +- -+ +- -+
C
C You may overlay any of the corresponding input or output terms
C to save storage by specifying the same actual parameter on the
C subroutine call.
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_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 matrix 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_X1(25) = INTEGER ARRAY
C FMS vector file attributes for {X1}.
C May be the same as LU_B1.
C
C (R ) LU_B2(25) = INTEGER ARRAY
C FMS vector file attributes for {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 ) NEWA22 = LOGICAL
C Indicated that [A22] has changed.
C
C (R ) INORM = Integer
C Norm for computing convergence:
C = 0, Infinity norm (max. abs. value)
C = 1, First norm
C = 2, Second norm
C
C (R ) MAXITR = Integer
C The maximum number of iterations to
C perform.
C
C (R ) ERRMAX = REAL*8
C = Maximum error for convergence
C
C (R ) NOFMSP = Logical
C Print FMS output in iteration loop.
C
C FMSCOM PARAMETERS:
C (RW) IACCOM = Multiply accumulate flag
C (R ) LUPR = FORTRAN unit for printing
C
C CALLED FROM:
C EXAMPLE_17
C
C SUBPROGRAMS USED:
C CNDF
C CNDS
C FMSMM
C
C ERROR CONDITIONS:
C None.
C
C HISTORY:
C
C NOTES:
C
C In this iterative approach, let {X1} and {X2} represent the
C exact solution and {DX1} and {DX2} the error.
C
C Expanding the partitioned system gives:
C
C [A11]({X1}+{DX1}) + [A12]({X2}+{DX2}) = {B1} (1)
C [A21]({X1}+{DX1}) + [A22]({X2}+{DX2}) = {B2} (2)
C
C or
C
C [A11]{X1} + [A11]{DX1} + [A12]{X2} + [A12]{DX2} = {B1} (3)
C [A21]{X1} + [A21]{DX1} + [A22]{X2} + [A22]{DX2} = {B2} (4)
C
C or
C
C ([A11]{X1}+[A12]{X2}-{B1}) + [A11]{DX1} = -[A12]{DX2} (5)
C ([A21]{X1}+[A22]{X2}-{B2}) + [A22]{DX2} = -[A21]{DX1} (6)
C
C Noting that the first term in (5) and (6) is zero gives
C
C [A11]{DX1} = -[A12]{DX2} (7)
C [A22]{DX2} = -[A21]{DX1} (8)
C
C which may be used to compute improved values of the error
C {DX1} and {DX2}.
C-----------------------------------------------------------------------
C FORMAL PARAMETERS
C-----------------------------------------------------------------------
INTEGER LU_A11(25)
INTEGER LU_F11(25)
C
INTEGER LU_A12(25)
C
INTEGER LU_A21(25)
C
INTEGER LU_A22(25)
INTEGER LU_F22(25)
C
INTEGER LU_B1 (25)
INTEGER LU_X1 (25)
C
INTEGER LU_B2 (25)
INTEGER LU_X2 (25)
C
LOGICAL NEWA11
LOGICAL NEWA22
INTEGER INORM
INTEGER MAXITR
REAL*8 ERRMAX
LOGICAL NOFMSP
C-----------------------------------------------------------------------
C LOCAL VARIABLES
C-----------------------------------------------------------------------
C
CHARACTER*7 MYNAME
PARAMETER (MYNAME='JACOBI2')
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
C Data type = complex:
INTEGER IDTYPE
PARAMETER (IDTYPE = 2)
C
C Skip operations during solving (no):
INTEGER ISKIP
PARAMETER (ISKIP=0)
C
C Number of solution steps:
INTEGER NSTEPS
PARAMETER (NSTEPS=10)
C
C Constants:
COMPLEX*16 C_ONE(1)
DATA C_ONE /( 1.0D0,0.0D0)/
COMPLEX*16 C_MONE
DATA C_MONE /(-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_TWO
DATA R_TWO / 2.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
INTEGER IPRF
INTEGER IPRS
INTEGER IPRVV
C
C Matrix dimensions:
INTEGER N1
INTEGER N2
INTEGER NRHS
REAL*8 RN1
REAL*8 RN2
REAL*8 RNRHS
C
C Convergence ratios:
REAL*8 R1
REAL*8 R2
C
C Temporary files:
INTEGER LU_DX1(25,2), LU_DX1_SAVE(25)
INTEGER LU_DX2(25,2), LU_DX2_SAVE(25)
INTEGER LUA0 (25)
DATA LUA0(1)/0/
INTEGER LUS0 (25)
INTEGER LUX0 (25)
INTEGER LUD (25)
C
C CPU, Wall and I/O time:
REAL*8 TCPU(NSTEPS)
REAL*8 T1, TCPU_T1 ,TCPU_T2
REAL*8 TWALL(NSTEPS)
REAL*8 T2, TWALL_T1 ,TWALL_T2
REAL*8 TIO(NSTEPS)
REAL*8 T3, TIO_T1 ,TIO_T2
C
C Floating point operations:
REAL*8 OPS(NSTEPS), TOPS
INTEGER NCALLS(NSTEPS)
C
C Sign for adding solution:
INTEGER ISIGN
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, MAXITR, ERRMAX
C
C Initialize operation counts and timers for each step:
DO I=1,NSTEPS
OPS (I) = R_ZERO
TCPU (I) = R_ZERO
TWALL (I) = R_ZERO
TIO (I) = R_ZERO
NCALLS(I) = 0
END DO
C
C Open temporary files for residual vectors:
CALL FMSOV (N1, IDTYPE, NRHS, 'LU_DX1_1' , LU_DX1(1,1))
CALL FMSOV (N1, IDTYPE, NRHS, 'LU_DX1_2' , LU_DX1(1,2))
CALL FMSOV (N2, IDTYPE, NRHS, 'LU_DX2_1' , LU_DX2(1,1))
CALL FMSOV (N2, IDTYPE, NRHS, 'LU_DX2_2' , LU_DX2(1,2))
LUD(1) = 0
C
CALL FMSTIM (TCPU_T1, TWALL_T1, TIO_T1)
CALL TIMER (0, R_ZERO, NCALLS, OPS, TCPU, TWALL, TIO)
C
C (1) Factor [A11] into [F11]=[L11][U11]:
IF(NEWA11) THEN
WRITE(LUPR,2001)
CALL CNDAF (LU_A11, C_ONE, 1, LUS0, 0, LUA0, LU_F11,
1 LUX0, LUX0, 0)
CALL TIMER (1, R_EIGHT*RN1*RN1*RN1/R_THREE,
1 NCALLS, OPS, TCPU, TWALL, TIO)
END IF
C
C (2) Solve [F11][X1] = [B1]:
WRITE(LUPR,2002)
CALL CNDS (LU_F11, LU_B1, LU_X1, NRHS, 0)
CALL TIMER (2, R_EIGHT*RN1*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (3) Factor [A22] into [F22]=[L22][U22]:
IF(NEWA22) THEN
WRITE(LUPR,2003)
CALL CNDAF (LU_A22, C_ONE, 1, LUS0, 0, LUA0, LU_F22,
1 LUX0, LUX0, 0)
CALL TIMER (3, R_EIGHT*RN2*RN2*RN2/R_THREE,
1 NCALLS, OPS, TCPU, TWALL, TIO)
END IF
C
C (4) Solve [F22]{X2} = {B2}:
WRITE(LUPR,2004)
CALL CNDS (LU_F22, LU_B2, LU_X2, NRHS, 0)
CALL TIMER (4, R_EIGHT*RN2*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C On the first pass, use {X1} and {X2} for {DX1}(KOLD) and
C {DX2}(KOLD):
K = 2
DO I=1,25
LU_DX1_SAVE(I) = LU_DX1(I,K)
LU_DX2_SAVE(I) = LU_DX2(I,K)
LU_DX1(I,K) = LU_X1(I)
LU_DX2(I,K) = LU_X2(I)
END DO
C
C Turn off FMS printing in loop, if required:
IF(NOFMSP) THEN
CALL FMSIGT ('IPRF, IPRF)
CALL FMSIST ('IPRF, 0)
CALL FMSIGT ('IPRS, IPRS)
CALL FMSIST ('IPRS, 0)
CALL FMSIGT ('IPRVV, IPRVV)
CALL FMSIST ('IPRVV, 0)
END IF
C
C Loop over iterations:
ISIGN = 1
DO 100 NITER = 1, MAXITR
KOLD = K
K = 3 - K
ISIGN = -ISIGN
C
C (5) Multiply {DX1}(K) = [A12]{DX2}(KOLD}:
IF(.NOT.NOFMSP) WRITE(LUPR,2005) K, KOLD
CALL FMSMM (LU_DX1(1,K), LUD, 0, 'N', LU_A12,
1 LU_DX2(1,KOLD))
LU_DX1(7,K) = 0
CALL TIMER (5, R_EIGHT*RN1*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (6) Solve [F11]{DX1}(K) = {DX1}(K):
IF(.NOT.NOFMSP) WRITE(LUPR,2006) K, K
CALL CNDS (LU_F11, LU_DX1(1,K), LU_DX1(1,K), NRHS, 0)
CALL TIMER (6, R_EIGHT*RN1*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (7) Compute {X1} = {X1} + (ISIGN)*{DX1} and compute norm R1:
IF(.NOT.NOFMSP) WRITE(LUPR,2007) K
CALL FMSVAN (LU_X1, ISIGN, LU_DX1(1,K), INORM, R1, 1)
IF(INORM .GT. 1) THEN
R1 = R1/RN1
ELSE IF(INORM .EQ. 2) THEN
R1 = R1/SQRT(RN1)
END IF
CALL TIMER (7, R_TWO*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (8) Multiply {DX2}(K) = [A21]{DX1}(KOLD):
IF(.NOT.NOFMSP) WRITE(LUPR,2008) K, KOLD
CALL FMSMM (LU_DX2(1,K), LUD, 0, 'T', LU_A21,
1 LU_DX1(1,KOLD))
LU_DX2(7,K) = 0
CALL TIMER (8, R_EIGHT*RN2*RN1*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (9) Solve [F22]{DX2}(K) = {DX2}(K):
IF(.NOT.NOFMSP) WRITE(LUPR,2009) K, K
CALL CNDS (LU_F22, LU_DX2(1,K), LU_DX2(1,K), NRHS, 0)
CALL TIMER (9, R_EIGHT*RN2*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
C (10) Compute {X2} = {X2} + (ISIGN)*{DX2}(K):
IF(.NOT. NOFMSP) WRITE(LUPR,2010) K
CALL FMSVAN (LU_X2, ISIGN, LU_DX2(1,K), INORM, R2, 1)
IF(INORM .GT. 1) THEN
R2 = R2/RN2
ELSE IF(INORM .EQ. 2) THEN
R2 = R2/SQRT(RN2)
END IF
CALL TIMER (10, R_TWO*RN2*RNRHS,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C
WRITE(LUPR,2050) NITER, R1, R2
IF(NITER .EQ. 1) THEN
C Redefine LU_DX1 and LU_DX2:
DO I=1,25
LU_DX1(I,KOLD) = LU_DX1_SAVE(I)
LU_DX2(I,KOLD) = LU_DX2_SAVE(I)
END DO
END IF
C
C Test for convergence:
IF( (R1 .LT. ERRMAX) .AND.
1 (R2 .LT. ERRMAX) ) GO TO 101
100 CONTINUE
101 CONTINUE
C
C Close temporary files:
CALL FMSCV (LU_DX1(1,1))
CALL FMSCV (LU_DX1(1,2))
CALL FMSCV (LU_DX2(1,1))
CALL FMSCV (LU_DX2(1,2))
C
C Print timings:
WRITE(LUPR,2020)
CALL FMSTIM (TCPU_T2, TWALL_T2, TIO_T2)
TOPS = R_ZERO
DO 200 ISTEP = 1,NSTEPS
NCALL = NCALLS(ISTEP)
IF(NCALL .EQ. 0) THEN
C
C Operation skipped:
WRITE(LUPR,2022) ISTEP
ELSE
TOPS = TOPS + OPS(ISTEP)
C
C CPU time and Mflops:
T1 = 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)
T3 = 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,2021) ISTEP, NCALL, T1, T2, T3, MF1, MF2, RATIO
END IF
200 CONTINUE
T1 = TCPU_T2 - TCPU_T1
T2 = TWALL_T2 - TWALL_T1
T3 = TIO_T2 - TIO_T1
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
IF(NOFMSP) THEN
C Restore print settings:
CALL FMSIST ('IPRF, IPRF)
CALL FMSIST ('IPRS, IPRS)
CALL FMSIST ('IPRVV, IPRVV)
END IF
WRITE(LUPR,2023) T1, T2, T3, MF1, MF2, RATIO
CALL FMSPOP(MYNAME)
RETURN
2000 FORMAT (/
1 ' Start of Subroutine JACOBI2',/
2 ' Rank of full matrix............=',I8/
3 ' Rank of submatrix #1...........=',I8/
4 ' Rank of submatrix #2...........=',I8/
5 ' # of right-hand-sides..........=',I8/
6 ' Maximum number of iterations...=',I8/
7 ' Maximum error norm.............=',E8.1/)
2001 FORMAT (/
1 ' STEP 1. Factor [A11] into [F11]'/
2 ' ==================================')
2002 FORMAT (/
1 ' STEP 2. Solve [F11]{X1} = {B1}'/
2 ' ===================================')
2003 FORMAT (/
1 ' STEP 3. Factor [A22] into [F22]'/
2 ' ==================================')
2004 FORMAT (/
1 ' STEP 4. Solve [F22]{X2} = {B2}'/
2 ' ===================================')
2005 FORMAT (/
1 ' STEP 5. Multiply {DX1}(',i1,') = [A12]{DX2}(',i1,')'/
2 ' ==========================================')
2006 FORMAT (/
1 ' STEP 6. Solve [F11]{DX1}(',I1,') = {DX1}(',I1,')'/
2 ' ==========================================')
2007 FORMAT (/
1 ' STEP 7. Add {X1} = {X1} + (ISIGN)*{DX1}(',I1,')'/
2 ' ================================================')
2008 FORMAT (/
1 ' STEP 8. Multiply {DX2}(',I1,') = [A21]{DX1}(',I1,')'/
2 ' ==========================================')
2009 FORMAT (/
1 ' STEP 9. Solve [F22]{DX2}(',I1,') = {DX2}(',I1,')'/
2 ' ==========================================')
2010 FORMAT (/
1 ' STEP 10. Add {X2} = {X2} + (ISIGN)*{DX2}(',I1,')'/
2 ' ================================================')
2020 FORMAT (/
1 ' Summary of Subroutine JACOBI2:'//
2 ' <----------TIMES (Sec.)----------->',
2 ' <--MFLOPS-->'/
2 ' STEP RUN CPU WALL I/O',
2 ' CPU WALL',
2 ' RATIO'/
3 ' ===== === ========= ========= =========',
3 ' ==== ====',
3 ' =====')
2021 FORMAT (2I6,3F13.3,2I8,F9.3)
2022 FORMAT (I6,12X,'Skipped')
2023 FORMAT (' TOTAL',6X,3F13.3,2I8,F9.3)
2050 FORMAT (
1 ' ITERATION=',I3,
2 ', {DX1}NORM=',E12.3,
3 ', {DX2}NORM=',E12.3)
END
C=======================================================================
SUBROUTINE TIMER (ICALL, OPSINC,
1 NCALLS, OPS, TCPU, TWALL, TIO)
C=======================================================================
INTEGER ICALL
REAL*8 OPSINC
INTEGER NCALLS(*)
REAL*8 OPS(*)
REAL*8 TCPU(*)
REAL*8 TWALL(*)
REAL*8 TIO(*)
COMMON /FMSTIMES/TCPU1, TWALL1, TIO1
REAL*8 TCPU1, TWALL1, TIO1
REAL*8 TCPU2, TWALL2, TIO2
IF(ICALL .EQ. 0) THEN
CALL FMSTIM (TCPU1, TWALL1, TIO1)
ELSE
CALL FMSTIM (TCPU2, TWALL2, TIO2)
NCALLS(ICALL) = NCALLS(ICALL) + 1
OPS (ICALL) = OPS (ICALL) + OPSINC
TCPU (ICALL) = TCPU (ICALL) + (TCPU2 - TCPU1 )
TWALL (ICALL) = TWALL (ICALL) + (TWALL2 - TWALL1)
TIO (ICALL) = TIO (ICALL) + (TIO2 - TIO1 )
TCPU1 = TCPU2
TWALL1 = TWALL2
TIO1 = TIO2
END IF
RETURN
END