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 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 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 )
CALL FMSOV (N2, IDTYPE, NRHS, 'LU_X2_OLD' , LU_X2_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