| Title: | Digital Extended Math Library |
| Notice: | Kit locations: 9.last (UNIX), 10.last (VMS) |
| Moderator: | RTL::CHAO FGREN |
| Created: | Mon Apr 30 1990 |
| Last Modified: | Tue Jun 03 1997 |
| Last Successful Update: | Fri Jun 06 1997 |
| Number of topics: | 324 |
| Total number of notes: | 1402 |
Hi,
I've got a combustion/CFD benchmark from GM, and am being asked what
can be done with it to optimize performance for SMP by the end of the
month. Consider, if you will (and have some time), the attached
subroutine.
Both the conjugate and bi-conjugate entry points are used. This code
and its callees consume circa 90% of the CPU time. If the inline
diagonal scaling preconditioner is used, KAP does an OK job of
parallelizing, but performance on an 8400 5/440 is alleged to be only
a bit better than R8000.
I have a couple of tacks that I can pursue on my own, but I'm wondering
about the prospects for plugging in a DXML solver. As near as I can
make out, the coefficient matrix has a bandwidth of 7, with the main
diagonal stored in AP (length 190488) and the smaller diagonals stored
in AMINUS and APLUS. So, something would need to be done to get the
coefficients into a form that DXML can deal with.
My questions are:
- does a DXML solution look feasible?
- do we have any idea what the parallel scaling would be like?
Thanks,
Chuck
SUBROUTINE SOLVER ! Entries <SOLVEB/SOLVEC>
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C+ +
C+ Bi-conjugate-gradient solver <SOLVEB> +
C+ Conjugate-gradient solver <SOLVEC> +
C+ +
C+---------------------------------------------------------------------+
C+ +
C+ Preconditioner: +
C+ +
C+ METHOD = 1 Diagonal scaling (DS ) +
C+ METHOD = 2 Symmetric successive over-relaxation (SSOR) +
C+ METHOD = 3 Incomplete Chelosky LU decomposition (ICLU) +
C+ +
C+ Convergence criteria: +
C+ +
C+ ICON = 1 for absolute convergence +
C+ ICON =-1 for relative convergence +
C+ ICON = 2 for both +
C+ ICON =-2 for either +
C+ +
C+---------------------------------------------------------------------+
C+ => (WRK2.7) must NOT be used in this routine (ref. CALCP) +
C+ => convergence tested via one-norm of the residual, rather +
C+ than the two-norm (Euclidean norm); which one is better? +
C+ => (XLU) should be equivalenced to something +
C+---------------------------------------------------------------------+
C+ Called by : CALCU , CALCE , CALCTE , CALCED , CALCY +
C+ CALCP +
C+ Calls : PSSOR , PICLU +
C+ : BSGI01 , BSGI02 , CSGI01 , CSGI02 +
C+---------------------------------------------------------------------+
C+ Modified on : Feb-22-93 +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
INCLUDE 'cd00.h'
INCLUDE 'cd10.h'
INCLUDE 'parm.h'
INCLUDE 'cntl.h'
INCLUDE 'geom.h'
INCLUDE 'data.h'
INCLUDE 'work.h'
DIMENSION VAR(-NBM:NCM),SU(0:NCM),AP(0:NCM)
DIMENSION AMINUS(0:NCMAX,3),APLUS(0:NCMAX,3)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
EQUIVALENCE (WRK2(0,1),AMINUS(0,1))
EQUIVALENCE (WRK2(0,4),APLUS (0,1))
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
DIMENSION R (NCMAX),P (0:NCMAX),Q (NCMAX),Z (0:NCMAX)
DIMENSION RC(NCMAX),PC(0:NCMAX),QC(NCMAX),ZC(0:NCMAX)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
EQUIVALENCE (WRK3(1,1),R (1))
EQUIVALENCE (WRK3(1,2),RC(1))
EQUIVALENCE (WRK3(1,3),P (1))
EQUIVALENCE (WRK3(1,4),PC(1))
EQUIVALENCE (WRK3(1,5),Q (1))
EQUIVALENCE (WRK3(1,6),QC(1))
EQUIVALENCE (WRK3(1,7),Z (1))
EQUIVALENCE (WRK3(1,8),ZC(1))
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
DIMENSION XLU(0:NCMAX,3,2),XLP(NCMAX)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
EQUIVALENCE (WRK3(1,9),XLP(1))
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
CHARACTER ABSREL*5,CHID*11
DATA GREAT,SMALL /1.0E+50,1.0E-50/
C
C
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C @@ BI-CONJUGATE-GRADIENT SOLVER @@
C @@ BI-CONJUGATE-GRADIENT SOLVER @@
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C
ENTRY SOLVEB(IPR,NBM,NCM,VAR,SU,AP,RES0,KTER,LUPRTO,CHID)
C
C
C ======================================================================
C (A) Pre-liminaries:
C ======================================================================
C
C
METHOD = ISOL(IPR)
KTER = 0
VAR(0) = 0.0
P (0) = 0.0
PC (0) = 0.0
Z (0) = 0.0
ZC (0) = 0.0
C
C
C ======================================================================
C (B) Compute initial residuals; establish convergence tolerance:
C ======================================================================
C
C
DO 110 IC = 1,NC
N11 = MAX(NEIGHB(IC,1,1),0)
N12 = MAX(NEIGHB(IC,1,2),0)
N21 = MAX(NEIGHB(IC,2,1),0)
N22 = MAX(NEIGHB(IC,2,2),0)
N31 = MAX(NEIGHB(IC,3,1),0)
N32 = MAX(NEIGHB(IC,3,2),0)
L11 = MIJ(LDTHRF(IC,1,1))
L12 = MIJ(LDTHRF(IC,1,2))
L21 = MIJ(LDTHRF(IC,2,1))
L22 = MIJ(LDTHRF(IC,2,2))
L31 = MIJ(LDTHRF(IC,3,1))
L32 = MIJ(LDTHRF(IC,3,2))
F11 = FIJ(LDTHRF(IC,1,1))
F12 = FIJ(LDTHRF(IC,1,2))
F21 = FIJ(LDTHRF(IC,2,1))
F22 = FIJ(LDTHRF(IC,2,2))
F31 = FIJ(LDTHRF(IC,3,1))
F32 = FIJ(LDTHRF(IC,3,2))
R (IC) = SU(IC)-AP(IC)*VAR(IC)
+ +AMINUS(IC,1)*VAR(N11)+APLUS(IC,1)*VAR(N12)
+ +AMINUS(IC,2)*VAR(N21)+APLUS(IC,2)*VAR(N22)
+ +AMINUS(IC,3)*VAR(N31)+APLUS(IC,3)*VAR(N32)
RC(IC) = SU(IC)-AP(IC)*VAR(IC)
+ +(APLUS(N11,L11)*F11+AMINUS(N11,L11)*(1.-F11))*VAR(N11)
+ +(APLUS(N12,L12)*F12+AMINUS(N12,L12)*(1.-F12))*VAR(N12)
+ +(APLUS(N21,L21)*F21+AMINUS(N21,L21)*(1.-F21))*VAR(N21)
+ +(APLUS(N22,L22)*F22+AMINUS(N22,L22)*(1.-F22))*VAR(N22)
+ +(APLUS(N31,L31)*F31+AMINUS(N31,L31)*(1.-F31))*VAR(N31)
+ +(APLUS(N32,L32)*F32+AMINUS(N32,L32)*(1.-F32))*VAR(N32)
110 CONTINUE
IF (NSGIS.GT.0) CALL BSGI01(IPR,NBM,NCM,VAR,R,RC,LUPRTO)
RES0 = 0.0
DO 140 IC = 1,NC
FAC0 = 1.0
IF (AP(IC).GE.GREAT) FAC0 = SMALL
RES0 = RES0+ABS(R(IC))*FAC0
140 CONTINUE
RES0 = RES0/SNORM(IPR)
ABSCON = ACON(IPR)
RELCON = RCON(IPR)*RES0
RSM = 0.0
IF (ICON(IPR).EQ. 1) RSM = ABSCON
IF (ICON(IPR).EQ.-1) RSM = RELCON
IF (ICON(IPR).EQ. 2) RSM = MIN(ABSCON,RELCON)
IF (ICON(IPR).EQ.-2) RSM = MAX(ABSCON,RELCON)
ABSREL = '(???)'
IF (RSM.EQ.ABSCON) ABSREL = '(ABS)'
IF (RSM.EQ.RELCON) ABSREL = '(REL)'
C
C
C ======================================================================
C (C) Precondition the initial residuals:
C ======================================================================
C
C
IF (METHOD.EQ.1) THEN
DO 210 IC = 1,NC
XLP(IC) = 1.0/AP(IC)
Z (IC) = XLP(IC)*R (IC)
ZC (IC) = XLP(IC)*RC(IC)
210 CONTINUE
ELSEIF (METHOD.EQ.2) THEN
CALL PSSOR(IPR,NCM,R ,Z ,AMINUS,APLUS,AP,XLP,-1,LUPRTO)
CALL PSSOR(IPR,NCM,R ,Z ,AMINUS,APLUS,AP,XLP, 0,LUPRTO)
CALL PSSOR(IPR,NCM,RC,ZC,AMINUS,APLUS,AP,XLP, 1,LUPRTO)
ELSEIF (METHOD.EQ.3) THEN
CALL PICLU(IPR,NCM,R ,Z ,AMINUS,APLUS,AP,XLU,XLP,-1,LUPRTO)
CALL PICLU(IPR,NCM,R ,Z ,AMINUS,APLUS,AP,XLU,XLP, 0,LUPRTO)
CALL PICLU(IPR,NCM,RC,ZC,AMINUS,APLUS,AP,XLU,XLP, 1,LUPRTO)
ENDIF
C
C
C ======================================================================
C (D) Define initial search directions:
C ======================================================================
C
C
RCZNEW = 0.0
DO 310 IC = 1,NC
P (IC) = Z (IC)
PC(IC) = ZC(IC)
RCZNEW = RCZNEW+RC(IC)*Z(IC)
310 CONTINUE
C
C
C ======================================================================
C (E) Iteration loop starts here:
C ======================================================================
C
C
410 KTER = KTER+1
IF (LSTRU.EQ.1 .OR. LSTRU.EQ.2) THEN
DO 415 IC = 1,NC
N11 = MAX(NEIGHB(IC,1,1),0)
N12 = MAX(NEIGHB(IC,1,2),0)
N21 = MAX(NEIGHB(IC,2,1),0)
N22 = MAX(NEIGHB(IC,2,2),0)
N31 = MAX(NEIGHB(IC,3,1),0)
N32 = MAX(NEIGHB(IC,3,2),0)
Q (IC) = AP(IC)*P(IC)
+ -AMINUS(IC,1)*P(N11)-APLUS(IC,1)*P(N12)
+ -AMINUS(IC,2)*P(N21)-APLUS(IC,2)*P(N22)
+ -AMINUS(IC,3)*P(N31)-APLUS(IC,3)*P(N32)
QC(IC) = AP(IC)*PC(IC)
+ -APLUS(N11,1)*PC(N11)-AMINUS(N12,1)*PC(N12)
+ -APLUS(N21,2)*PC(N21)-AMINUS(N22,2)*PC(N22)
+ -APLUS(N31,3)*PC(N31)-AMINUS(N32,3)*PC(N32)
415 CONTINUE
ELSE
DO 420 IC = 1,NC
N11 = MAX(NEIGHB(IC,1,1),0)
N12 = MAX(NEIGHB(IC,1,2),0)
N21 = MAX(NEIGHB(IC,2,1),0)
N22 = MAX(NEIGHB(IC,2,2),0)
N31 = MAX(NEIGHB(IC,3,1),0)
N32 = MAX(NEIGHB(IC,3,2),0)
L11 = MIJ(LDTHRF(IC,1,1))
L12 = MIJ(LDTHRF(IC,1,2))
L21 = MIJ(LDTHRF(IC,2,1))
L22 = MIJ(LDTHRF(IC,2,2))
L31 = MIJ(LDTHRF(IC,3,1))
L32 = MIJ(LDTHRF(IC,3,2))
F11 = FIJ(LDTHRF(IC,1,1))
F12 = FIJ(LDTHRF(IC,1,2))
F21 = FIJ(LDTHRF(IC,2,1))
F22 = FIJ(LDTHRF(IC,2,2))
F31 = FIJ(LDTHRF(IC,3,1))
F32 = FIJ(LDTHRF(IC,3,2))
Q (IC) = AP(IC)*P(IC)
+ -AMINUS(IC,1)*P(N11)-APLUS(IC,1)*P(N12)
+ -AMINUS(IC,2)*P(N21)-APLUS(IC,2)*P(N22)
+ -AMINUS(IC,3)*P(N31)-APLUS(IC,3)*P(N32)
QC(IC) = AP(IC)*PC(IC)
+ -(APLUS(N11,L11)*F11+AMINUS(N11,L11)*(1.-F11))*PC(N11)
+ -(APLUS(N12,L12)*F12+AMINUS(N12,L12)*(1.-F12))*PC(N12)
+ -(APLUS(N21,L21)*F21+AMINUS(N21,L21)*(1.-F21))*PC(N21)
+ -(APLUS(N22,L22)*F22+AMINUS(N22,L22)*(1.-F22))*PC(N22)
+ -(APLUS(N31,L31)*F31+AMINUS(N31,L31)*(1.-F31))*PC(N31)
+ -(APLUS(N32,L32)*F32+AMINUS(N32,L32)*(1.-F32))*PC(N32)
420 CONTINUE
ENDIF
IF (NSGIS.GT.0) CALL BSGI02(IPR,NBM,NCM,Q,QC,P,PC,LUPRTO)
PCTQ = 0.0
DO 450 IC = 1,NC
450 PCTQ = PCTQ+PC(IC)*Q(IC)
IF (ABS(PCTQ).LE.SMALL) THEN
WRITE(LUPRTO,460)
460 FORMAT(' = M =>',' Solver returns the old solution')
GO TO 900
ENDIF
C
C
C ======================================================================
C (F) Update solution & residuals:
C ======================================================================
C
C
ALPHA = RCZNEW/PCTQ
DO 510 IC = 1,NC
VAR(IC) = VAR(IC)+ALPHA*P (IC)
R (IC) = R (IC)-ALPHA*Q (IC)
RC (IC) = RC (IC)-ALPHA*QC(IC)
510 CONTINUE
C
C
C ======================================================================
C (G) Precondition the updated residuals:
C ======================================================================
C
C
IF (METHOD.EQ.1) THEN
DO 610 IC = 1,NC
Z (IC) = XLP(IC)*R (IC)
ZC(IC) = XLP(IC)*RC(IC)
610 CONTINUE
ELSEIF (METHOD.EQ.2) THEN
CALL PSSOR(IPR,NCM,R ,Z ,AMINUS,APLUS,AP,XLP,0,LUPRTO)
CALL PSSOR(IPR,NCM,RC,ZC,AMINUS,APLUS,AP,XLP,1,LUPRTO)
ELSEIF (METHOD.EQ.3) THEN
CALL PICLU(IPR,NCM,R ,Z ,AMINUS,APLUS,AP,XLU,XLP,0,LUPRTO)
CALL PICLU(IPR,NCM,RC,ZC,AMINUS,APLUS,AP,XLU,XLP,1,LUPRTO)
ENDIF
RCZOLD = RCZNEW
RCZNEW = 0.0
DO 660 IC = 1,NC
660 RCZNEW = RCZNEW+RC(IC)*Z(IC)
C
C
C ======================================================================
C (H) Update the search directions:
C ======================================================================
C
C
BETA = RCZNEW/RCZOLD
DO 710 IC = 1,NC
P (IC) = Z (IC)+BETA*P (IC)
PC(IC) = ZC(IC)+BETA*PC(IC)
710 CONTINUE
C
C
C ======================================================================
C (I) Convergence test:
C ======================================================================
C
C
RESK = 0.0
DO 810 IC = 1,NC
FACK = 1.0
IF (AP(IC).GE.GREAT) FACK = SMALL
RESK = RESK+ABS(R(IC))*FACK
810 CONTINUE
RESK = RESK/SNORM(IPR)
C
C --- didn't converge within the iteration limit
C
IF (KTER.EQ.MCON(IPR)) THEN
WRITE(LUPRTO,820) CHID,KTER,RESK,RES0,SNORM(IPR),ABSREL
820 FORMAT(' = W => ',A11,'; iter = ',I4,'; res = ',1P,E9.3
+ ,'; res0 = ',E9.3,'; snorm = ',E9.3,'; ',A5)
GO TO 900
ENDIF
IF (RESK.GT.RSM) GO TO 410
C
C --- process converged
C
WRITE(LUPRTO,830) CHID,KTER,RESK,RES0,SNORM(IPR),ABSREL
830 FORMAT(' ',A11,'; iter = ',I4,'; res = ',1P,E9.3
+ ,'; res0 = ',E9.3,'; snorm = ',E9.3,'; ',A5)
900 CONTINUE
C
RETURN
C
C
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C @@ CONJUGATE-GRADIENT SOLVER @@
C @@ CONJUGATE-GRADIENT SOLVER @@
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C
ENTRY SOLVEC(IPR,NBM,NCM,VAR,SU,AP,RES0,KTER,LUPRTO,CHID)
C
C
C ======================================================================
C (A) Pre-liminaries:
C ======================================================================
C
C
METHOD = ISOL(IPR)
KTER = 0
VAR(0) = 0.0
P (0) = 0.0
Z (0) = 0.0
C
C
C ======================================================================
C (B) Compute initial residual; establish convergence tolerance:
C ======================================================================
C
C
DO 1110 IC = 1,NC
N11 = NEIGHB(IC,1,1)
N12 = NEIGHB(IC,1,2)
N21 = NEIGHB(IC,2,1)
N22 = NEIGHB(IC,2,2)
N31 = NEIGHB(IC,3,1)
N32 = NEIGHB(IC,3,2)
Comment ----------------------------------------------------------------
C Although Nxx is negative on a boundary, APLUS/AMINUS was set
C explicitly to zero in ASMATS/ASMATP; therefore, the residual
C R is computed correctly
Comment ----------------------------------------------------------------
R(IC) = AMINUS(IC,1)*VAR(N11)+APLUS(IC,1)*VAR(N12)
+ +AMINUS(IC,2)*VAR(N21)+APLUS(IC,2)*VAR(N22)
+ +AMINUS(IC,3)*VAR(N31)+APLUS(IC,3)*VAR(N32)
+ +SU(IC)-AP(IC)*VAR(IC)
1110 CONTINUE
IF (NSGIS.GT.0) CALL CSGI01(IPR,NBM,NCM,VAR,R,LUPRTO)
RES0 = 0.0
DO 1140 IC = 1,NC
1140 RES0 = RES0+ABS(R(IC))
RES0 = RES0/SNORM(IPR)
ABSCON = ACON(IPR)
RELCON = RCON(IPR)*RES0
RSM = 0.0
IF (ICON(IPR).EQ. 1) RSM = ABSCON
IF (ICON(IPR).EQ.-1) RSM = RELCON
IF (ICON(IPR).EQ. 2) RSM = MIN(ABSCON,RELCON)
IF (ICON(IPR).EQ.-2) RSM = MAX(ABSCON,RELCON)
ABSREL = '(???)'
IF (RSM.EQ.ABSCON) ABSREL = '(ABS)'
IF (RSM.EQ.RELCON) ABSREL = '(REL)'
C
C
C ======================================================================
C (C) Precondition the initial residual:
C ======================================================================
C
C
IF (METHOD.EQ.1) THEN
DO 1210 IC = 1,NC
XLP(IC) = 1.0/AP(IC)
Z (IC) = XLP(IC)*R(IC)
1210 CONTINUE
ELSEIF (METHOD.EQ.2) THEN
CALL PSSOR(IPR,NCM,R,Z,AMINUS,APLUS,AP,XLP,-1,LUPRTO)
CALL PSSOR(IPR,NCM,R,Z,AMINUS,APLUS,AP,XLP, 0,LUPRTO)
ELSEIF (METHOD.EQ.3) THEN
CALL PICLU(IPR,NCM,R,Z,AMINUS,APLUS,AP,XLU,XLP,-1,LUPRTO)
CALL PICLU(IPR,NCM,R,Z,AMINUS,APLUS,AP,XLU,XLP, 0,LUPRTO)
ENDIF
C
C
C ======================================================================
C (D) Define initial search direction:
C ======================================================================
C
C
RZNEW = 0.0
DO 1310 IC = 1,NC
P(IC) = Z(IC)
RZNEW = RZNEW+R(IC)*Z(IC)
1310 CONTINUE
C
C
C ======================================================================
C (E) Iteration loop starts here:
C ======================================================================
C
C
1410 KTER = KTER+1
DO 1420 IC = 1,NC
N11 = NEIGHB(IC,1,1)
N12 = NEIGHB(IC,1,2)
N21 = NEIGHB(IC,2,1)
N22 = NEIGHB(IC,2,2)
N31 = NEIGHB(IC,3,1)
N32 = NEIGHB(IC,3,2)
Comment ---------------------------- to prevent out-of-bound message ---
C N11 = MAX(NEIGHB(IC,1,1),0)
C N12 = MAX(NEIGHB(IC,1,2),0)
C N21 = MAX(NEIGHB(IC,2,1),0)
C N22 = MAX(NEIGHB(IC,2,2),0)
C N31 = MAX(NEIGHB(IC,3,1),0)
C N32 = MAX(NEIGHB(IC,3,2),0)
Comment ---------------------------- to prevent out-of-bound message ---
Q(IC) = AP(IC)*P(IC)
+ -AMINUS(IC,1)*P(N11)-APLUS(IC,1)*P(N12)
+ -AMINUS(IC,2)*P(N21)-APLUS(IC,2)*P(N22)
+ -AMINUS(IC,3)*P(N31)-APLUS(IC,3)*P(N32)
1420 CONTINUE
IF (NSGIS.GT.0) CALL CSGI02(IPR,NBM,NCM,Q,P,LUPRTO)
PTQ = 0.0
DO 1450 IC = 1,NC
1450 PTQ = PTQ+P(IC)*Q(IC)
IF (ABS(PTQ).LE.SMALL) THEN
WRITE(LUPRTO,1460)
1460 FORMAT(' = M =>',' Solver returns the old solution')
GO TO 1900
ENDIF
C
C
C ======================================================================
C (F) Update solution & residual:
C ======================================================================
C
C
ALPHA = RZNEW/PTQ
DO 1510 IC = 1,NC
VAR(IC) = VAR(IC)+ALPHA*P(IC)
R (IC) = R (IC)-ALPHA*Q(IC)
1510 CONTINUE
C
C
C ======================================================================
C (G) Precondition the updated residual:
C ======================================================================
C
C
IF (METHOD.EQ.1) THEN
DO 1610 IC = 1,NC
1610 Z(IC) = XLP(IC)*R(IC)
ELSEIF (METHOD.EQ.2) THEN
CALL PSSOR(IPR,NCM,R,Z,AMINUS,APLUS,AP,XLP,0,LUPRTO)
ELSEIF (METHOD.EQ.3) THEN
CALL PICLU(IPR,NCM,R,Z,AMINUS,APLUS,AP,XLU,XLP,0,LUPRTO)
ENDIF
RZOLD = RZNEW
RZNEW = 0.0
DO 1660 IC = 1,NC
1660 RZNEW = RZNEW+R(IC)*Z(IC)
C
C
C ======================================================================
C (H) Update the search direction:
C ======================================================================
C
C
BETA = RZNEW/RZOLD
DO 1710 IC = 1,NC
1710 P(IC) = Z(IC)+BETA*P(IC)
C
C
C ======================================================================
C (I) Convergence test:
C ======================================================================
C
C
RESK = 0.0
DO 1810 IC = 1,NC
1810 RESK = RESK+ABS(R(IC))
RESK = RESK/SNORM(IPR)
C
C --- didn't converge within the iteration limit
C
IF (KTER.EQ.MCON(IPR)) THEN
WRITE(LUPRTO,1820) CHID,KTER,RESK,RES0,SNORM(IPR),ABSREL
1820 FORMAT(' = W => ',A11,'; iter = ',I4,'; res = ',1P,E9.3
+ ,'; res0 = ',E9.3,'; snorm = ',E9.3,'; ',A5)
GO TO 1900
ENDIF
IF (RESK.GT.RSM) GO TO 1410
C
C --- process converged
C
WRITE(LUPRTO,1830) CHID,KTER,RESK,RES0,SNORM(IPR),ABSREL
1830 FORMAT(' ',A11,'; iter = ',I4,'; res = ',1P,E9.3
+ ,'; res0 = ',E9.3,'; snorm = ',E9.3,'; ',A5)
1900 CONTINUE
C
RETURN
C
C
END
| T.R | Title | User | Personal Name | Date | Lines |
|---|