Next: Installation Guide Up: NSPCG User's Guide Previous: Sample Usage

Using NSPCG in Matrix Format-Free Mode

When the NSPCG subroutine is called, the user is required to use one of the allowed data structures. However, it is possible to call the acceleration routines directly in which case the user supplies customized routines for performing certain matrix operations. In this manner, an iterative algorithm can be designed which is suitable to a particular application. Also, storage demands may preclude copying the matrix into one of NSPCG's allowed matrix formats. In this case, it may be necessary to use the data format available from the application code itself when writing the subroutines for matrix operations. The subroutines which the user must supply to perform certain matrix operations are denoted SUBA, SUBAT, SUBQL, SUBQR, SUBQLT, SUBQRT, SUBQ, and SUBADP. Not all of these subroutines are needed for every accelerator. The user should consult the calling sequences for the accelerators to determine which subroutines need to be supplied.

SUBA and SUBAT are routines for performing y=Ax and y=A^Tx, respectively. Their calling sequences are

       SUBA (COEF,JCOEF,WFAC,JWFAC,N,X,Y)

and

       SUBAT (COEF,JCOEF,WFAC,JWFAC,N,X,Y)

where

COEF: is a real vector. [real, input]
JCOEF: is an integer vector. [integer, input]
WFAC: is a real vector. [real, input]
JWFAC: is an integer vector. [integer, input]
N: is the order of the linear system. [integer, input]
X: is a real vector of length N containing the input vector (the vector to be multiplied by) on input. It should be unchanged on output. [real, input]
Y: is a real vector of length N containing the output vector (resultant product vector) on output. [real, output]

The vectors COEF, JCOEF, WFAC, and JWFAC are simply passed by the acceleration routines to the matrix-vector product routines SUBA and SUBAT. Hence, the user may use these vectors for any purpose inside SUBA or SUBAT. For example, they can be used to represent the operator A or workspace needed to compute Ax or A^Tx. If any of these vectors is not needed, it is not necessary to declare them. SUBQL, SUBQR, SUBQLT, and SUBQRT are routines for performing y=Q_L^-1x, y=Q_R^-1x, y=Q_L^-Tx, and y=Q_R^-Tx, respectively. These routines assume that the preconditioning operator Q can be written as Q=Q_LQ_R and that the preconditioned matrix is Q_L^-1AQ_R^-1. Note that Q_R=I for left preconditioning and Q_L=I for right preconditioning. These four subroutines have the same calling sequences as SUBA and SUBAT. As with these two subroutines, the parameters COEF, JCOEF, WFAC, and JWFAC can be used for any purpose by the user since they are only passed along by the accelerator. Also, X is the input vector (the right-hand-side vector in the system Q_Ly=x) and Y is the output vector (the solution to the system). The NSPCG library contains a routine called COPY with the argument list given above which simply performs a copy of the vector X into vector Y. Hence if the user wishes to specify that Q_R=I (i.e., only a left preconditioner is desired), then COPY should be used for SUBQR and SUBQRT. SUBQ is used only in the SOR routine and computes an SOR iteration. Thus SUBQ computes

u⁽ⁿ⁺¹⁾ = Gu⁽ⁿ⁾ + k

where

G = I - Q^-1A

and

k = Q^-1b

Hence SUBQ computes

$\begin{displaymath}u^{(n+1)} = \left(\frac{1}{\omega}D-C_L \right)^{-1}\left[\left(\frac{1}{\omega}-1 \right)Du^{(n)} + C_Uu^{(n)}+b \right]\end{displaymath}$

where

A = D-C_L-C_U

Its calling sequence is
SUBQ (COEF,JCOEF,WFAC,JWFAC,N,U,RHS,UNEW)
The interpretation of the parameters COEF, JCOEF, WFAC, JWFAC, and N is the same as before. U is a real vector of length N which contains the current solution vector u⁽ⁿ⁾ on input, and it should not be changed on output. UNEW is a real vector of length N which should contain on output the new solution vector u⁽ⁿ⁺¹⁾ after one SOR iteration. RHS is a real vector of length N containing the right-hand-side b of the linear system on input; it does not need to be preserved on output and so can be used for workspace if needed. SUBADP is used only in the adaptive SRCG and SRSI routines and is used to calculate certain quantities needed for the adaptive procedures for $\omega$ . Its calling sequence is
SUBADP (COEF,JCOEF,WFAC,JWFAC,N,P,R,PDP,PLDUP)
The interpretation of the parameters COEF, JCOEF, WFAC, JWFAC, and N is the same as before. P is a real vector of length N which contains on input a direction vector from the conjugate gradient or Chebyshev algorithms. It should not be changed on output. R is a real vector of length N from the accelerator routine which can be used for workspace. PDP is a real constant which contains on output the quantity (p,Dp) where A=D-C_L-C_U. PLDUP is a real constant which contains on output the quantity (p,C_LD^-1C_Up). Note that the user is free to define D as the diagonal of A or as the block diagonal part of A, as long as the SUBQL, SUBQR, SUBQLT, and SUBQRT routines are consistently defined. If the user does not wish to adapt on $\omega$ , the nonadaptive CG and SI routines should be called instead of the adaptive ones since they do not require the user to supply SUBADP. If the SOR, SRCG, or SRSI accelerators are called, the user may need the value of $\omega$ used in the accelerator for the lower level matrix operation routines. The user should include the following two lines of code in each subroutine requiring $\omega$ :
LOGICAL OMGADP COMMON / ITCOM5 / OMGADP, OMEGA, ALPHAB, BETAB, FFF, SPECR
The variable OMEGA contains the value of $\omega$ . The calling sequences for the acceleration routines with the corresponding $\langle$ accel $\rangle$ entries on the left are as follows:

CG
CGW (SUBA,SUBQL, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
SI
SIW (SUBA,SUBQL, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
SOR
SORW (SUBA,SUBQ, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
SRCG
SRCGW (SUBA,SUBQL,SUBADP, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
SRSI
SRSIW (SUBA,SUBQL,SUBADP, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
BASIC
BASICW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
ME
MEW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
CGNR
CGNRW (SUBA,SUBAT,SUBQL,SUBQLT,SUBQR,SUBQRT, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
LSQR
LSQRW (SUBA,SUBAT,SUBQL,SUBQLT,SUBQR,SUBQRT, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
ODIR
ODIRW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
OMIN
OMINW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
ORES
ORESW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
IOM
IOMW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
GMRES
GMRESW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
USYMLQ
USLQW (SUBA,SUBAT,SUBQL,SUBQLT,SUBQR,SUBQRT, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
USYMQR
USQRW (SUBA,SUBAT,SUBQL,SUBQLT,SUBQR,SUBQRT, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
LANDIR
LDIRW (SUBA,SUBAT,SUBQL,SUBQLT,SUBQR,SUBQRT, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
LANMIN
LMINW (SUBA,SUBAT,SUBQL,SUBQLT,SUBQR,SUBQRT, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)
LANRES
LRESW (SUBA,SUBAT,SUBQL,SUBQLT,SUBQR,SUBQRT, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)

BCGS
BCGSW (SUBA,SUBQL,SUBQR, COEF,JCOEF,WFAC,JWFAC,N,U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER)

where

SUBA
is the name of a matrix-vector multiplication routine which must be declared EXTERNAL in the user's calling routine. It computes y=Ax. [name]
SUBAT
is the name of a transpose matrix-vector multiplication routine which must be declared EXTERNAL in the user's calling routine. It computes y=A^Tx. [name]
SUBQL
is the name of the left-preconditioning routine which must be declared EXTERNAL in the user's calling routine. It solves Q_Ly=x for y given x. [name]
SUBQR
is the name of the right-preconditioning routine which must be declared EXTERNAL in the user's calling routine. It solves Q_Ry=x for y given x. [name]
SUBQLT
is the name of the transpose left-preconditioning routine which must be declared EXTERNAL in the user's calling routine. It solves Q_L^Ty=x for y given x. [name]
SUBQRT
is the name of the transpose right-preconditioning routine which must be declared EXTERNAL in the user's calling routine. It solves Q_R^Ty=x for y given x. [name]
SUBQ
is the name of the SOR iteration routine which must be declared EXTERNAL in the user's calling routine. It computes an SOR iteration. [name]
SUBADP
is the name of the SSOR adaptive procedure routine which must be declared EXTERNAL in the user's calling routine. It computes parameters necessary for adapting on $\omega$ for the SRCGW and SRSIW acceleration routines. [name]
COEF
is a real vector which is passed to the lower level matrix routines. [real, input]
JCOEF
is an integer vector which is passed to the lower level matrix routines. [integer, input]
WFAC
is a real vector which is passed to the lower level matrix routines. [real, input]
JWFAC
is an integer vector which is passed to the lower level matrix routines. [integer, input]
N
is the order of the linear system. [integer, input]
U
is a real vector of length N containing the current estimate of the solution to the linear system on input. The user must supply an initial guess, which can be all zeros if no information is known. On output, U contains an improved estimate of the solution vector if convergence is achieved. [real, input/output]
UBAR
is a vector of length N containing the exact solution to the linear system, if known. This vector is used only if the exact stopping test is used. Otherwise, it may be dimensioned to be of length one. [real, input]
RHS
is a vector of length N containing the right-hand-side of the matrix problem. [real, input]
WKSP
is a real vector of length NW used for workspace. [real, input]
NW
is an integer variable indicating the length of the WKSP vector on input and the amount used on output. If insufficient real workspace is provided, NW is set on output to the amount of workspace needed at the point execution terminated. [integer, input/output]
IPARM
is an integer vector of length 25 giving various integer parameters which are described in Section 7. Only parameters 1 through 11 and 22 are relevant when the package is called at this level. [integer, input/output]
RPARM
is a real vector of length 16 giving various real parameters which are described in Section 7. Only parameters 1 through 12 are relevant when the package is called at this level. [real, input/output]
IER
is an error flag. A nonzero value on output indicates a fatal or warning error was detected during the iteration. See Section13 on error conditions for the meaning of its assigned values. [integer, output]
Some examples of the use of NSPCG in matrix format-free mode follow:

Example 1:

In this example, NSPCG was used to solve the linear system Ax=bwhich resulted from discretizing the problem

$\begin{displaymath}\left\{\begin{array}{cc} u_{xx} + 2u_{yy} = 0 & \mbox{on S $... ... \\ u = 1 + xy & \mbox{on boundary of S} \end{array} \right. \end{displaymath}$

using the standard 5-point central finite difference formula with a mesh size of $h = \frac{1}{11}$ . The finite difference stencil at node (i,j) is thus

6u_i,j-u_i-1,j-u_i+1,j-2u_i,j+1-2u_i,j-1 = 0

This resulted in a symmetric and positive definite matrix of order 100 with five nonzero diagonals. The matrix coefficients were not stored in an array, but a subroutine called MULT was written to compute the matrix-vector product, y=Ax. The iterative method used was Richardson's method with conjugate gradient acceleration. Note that for Richardson's method, Q=I, so COPY was used for the SUBQL parameter in the CGW parameter list. The computer used to run the program was a Cyber 170/750 and the compiler used was FTN5. The program to set up the problem and the output resulting from this call to CGW is given below:

PROGRAM MAIN (OUTPUT,TAPE6=OUTPUT) C C ... ARRAY DECLARATIONS. C REAL RHS(100), U(100), WKSP(600), UBAR(1), RPARM(30) INTEGER IPARM(30) EXTERNAL MULT, COPY C NW = 600 C C ... GENERATE RHS. C NX = 10 NY = 10 N = NX*NY H = 1.0/FLOAT(NX + 1) DO 10 I = 1,N RHS(I) = 0.0 10 CONTINUE K = 0 DO 30 J = 1,NY Y = FLOAT(J)*H DO 25 I = 1,NX X = FLOAT(I)*H K = K + 1 IF (J .EQ. 1) RHS(K) = RHS(K) + 2.0 IF (J .EQ. NY) RHS(K) = RHS(K) + 2.0*(1.0 + X) IF (I .EQ. 1) RHS(K) = RHS(K) + 1.0 IF (I .EQ. NX) RHS(K) = RHS(K) + 1.0 + Y 25 CONTINUE 30 CONTINUE CALL DFAULT (IPARM,RPARM) C C ... NOW, RESET SOME DEFAULT VALUES. C IPARM(3) = 3 RPARM(1) = 1.0E-8 C C ... GENERATE AN INITIAL GUESS FOR U AND CALL NSPCG. C CALL VFILL (N,U,0.0) C CALL CGW (MULT,COPY,COEF,JCOEF,WFAC,JWFAC,N, A U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER) STOP END SUBROUTINE MULT (COEF,JCOEF,WFAC,JWFAC,N,X,Y) REAL X(N), Y(N) NX = 10 C C ... COMPUTE PRODUCT AS IF FIRST SUPERDIAGONAL AND FIRST C SUBDIAGONAL WERE FULL. C Y(1) = 6.0*X(1) - X(2) - 2.0*X(NX+1) DO 10 I = 2,NX Y(I) = 6.0*X(I) - X(I+1) - X(I-1) - 2.0*X(I+NX) 10 CONTINUE DO 15 I = NX+1,N-NX Y(I) = 6.0*X(I) - X(I+1) - X(I-1) - 2.0*(X(I+NX) + X(I-NX)) 15 CONTINUE DO 20 I = N-NX+1,N-1 Y(I) = 6.0*X(I) - X(I+1) - X(I-1) - 2.0*X(I-NX) 20 CONTINUE Y(N) = 6.0*X(N) - X(N-1) - 2.0*X(N-NX) C C ... MAKE CORRECTIONS TO Y VECTOR FOR ZEROS IN FIRST SUPERDIAGONAL C AND FIRST SUBDIAGONAL. C DO 25 I = NX,N-NX,NX Y(I) = Y(I) + X(I+1) 25 CONTINUE DO 30 I = NX+1,N-NX+1,NX Y(I) = Y(I) + X(I-1) 30 CONTINUE RETURN END
INITIAL ITERATIVE PARAMETERS GENERAL AND ACCELERATION PARAMETERS IPARM( 1) = 2 (NTEST ) IPARM( 2) = 100 (ITMAX ) IPARM( 3) = 3 (LEVEL ) IPARM( 4) = 6 (NOUT ) IPARM( 5) = 0 (IDGTS ) IPARM( 6) = 1 (MAXADP) IPARM( 7) = 1 (MINADP) IPARM( 8) = 1 (IOMGAD) IPARM( 9) = 5 (NS1 ) IPARM(10) = 100000 (NS2 ) IPARM(11) = 0 (NS3 ) RPARM( 1) = .10000000E-07 (ZETA ) RPARM( 2) = .20000000E+01 (EMAX ) RPARM( 3) = .10000000E+01 (EMIN ) RPARM( 4) = .75000000E+00 (FF ) RPARM( 5) = .75000000E+00 (FFF ) RPARM( 6) = .00000000E+00 (TIMIT ) RPARM( 7) = .00000000E+00 (DIGIT1) RPARM( 8) = .00000000E+00 (DIGIT2) RPARM( 9) = .10000000E+01 (OMEGA ) RPARM(10) = .00000000E+00 (ALPHAB) RPARM(11) = .25000000E+00 (BETAB ) RPARM(12) = .00000000E+00 (SPECR )
CG INTERMEDIATE OUTPUT AFTER EACH ITERATION ITERATION CONVERGENCE EMAX EMIN N S TEST 0 0 .13900E+02 .20000E+01 .10000E+01 1 1 .56024E+00 .37349E+01 .37349E+01 2 2 .49004E+00 .61940E+01 .19544E+01 3 3 .51510E+00 .74467E+01 .12374E+01 4 4 .52558E+00 .86621E+01 .86417E+00 5 5 .63725E+00 .97284E+01 .61624E+00 6 6 .89475E+00 .10417E+02 .39239E+00 7 7 .66075E+00 .10698E+02 .32758E+00 8 8 .68211E+00 .10905E+02 .28757E+00 9 9 .48761E+00 .11011E+02 .26502E+00 10 10 .24106E+00 .11096E+02 .25229E+00 11 11 .15764E+00 .11146E+02 .24748E+00 12 12 .94345E-01 .11245E+02 .24514E+00 13 13 .49491E-01 .11486E+02 .24418E+00 14 14 .29955E-01 .11621E+02 .24383E+00 15 15 .24829E-01 .11659E+02 .24363E+00 16 16 .23087E-01 .11676E+02 .24346E+00 17 17 .15022E-01 .11699E+02 .24323E+00 18 18 .83114E-02 .11718E+02 .24312E+00 19 19 .43446E-02 .11734E+02 .24308E+00 20 20 .31278E-02 .11743E+02 .24307E+00 21 21 .18751E-02 .11751E+02 .24306E+00 22 22 .16765E-02 .11754E+02 .24305E+00 23 23 .93634E-03 .11755E+02 .24304E+00 24 24 .47839E-03 .11756E+02 .24304E+00 25 25 .28610E-03 .11756E+02 .24304E+00 26 26 .16126E-03 .11757E+02 .24304E+00 27 27 .85505E-04 .11757E+02 .24304E+00 28 28 .56106E-04 .11757E+02 .24304E+00 29 29 .28246E-04 .11757E+02 .24304E+00 30 30 .16505E-04 .11757E+02 .24304E+00 31 31 .96348E-05 .11757E+02 .24304E+00 32 32 .62848E-05 .11757E+02 .24304E+00 33 33 .33737E-05 .11757E+02 .24304E+00 34 34 .12847E-05 .11757E+02 .24304E+00 35 35 .63580E-06 .11757E+02 .24304E+00 36 36 .25529E-06 .11757E+02 .24304E+00 37 37 .12834E-06 .11757E+02 .24304E+00 38 38 .64785E-07 .11757E+02 .24304E+00 39 39 .24616E-07 .11757E+02 .24304E+00 40 40 .10040E-07 .11757E+02 .24304E+00 41 41 .41463E-08 .11757E+02 .24304E+00 CG HAS CONVERGED IN 41 ITERATIONS
FINAL ITERATIVE PARAMETERS GENERAL AND ACCELERATION PARAMETERS IPARM( 1) = 2 (NTEST ) IPARM( 2) = 41 (ITMAX ) IPARM( 3) = 3 (LEVEL ) IPARM( 4) = 6 (NOUT ) IPARM( 5) = 0 (IDGTS ) IPARM( 6) = 1 (MAXADP) IPARM( 7) = 1 (MINADP) IPARM( 8) = 1 (IOMGAD) IPARM( 9) = 5 (NS1 ) IPARM(10) = 100000 (NS2 ) IPARM(11) = 0 (NS3 ) RPARM( 1) = .10000000E-07 (ZETA ) RPARM( 2) = .11756958E+02 (EMAX ) RPARM( 3) = .24304216E+00 (EMIN ) RPARM( 4) = .75000000E+00 (FF ) RPARM( 5) = .75000000E+00 (FFF ) RPARM( 6) = .59400000E+00 (TIMIT ) RPARM( 7) = .83823401E+01 (DIGIT1) RPARM( 8) = .90374486E+01 (DIGIT2) RPARM( 9) = .10000000E+01 (OMEGA ) RPARM(10) = .00000000E+00 (ALPHAB) RPARM(11) = .25000000E+00 (BETAB ) RPARM(12) = .00000000E+00 (SPECR )

Example 2:

In this example, the same problem given in Example 1 was solved using the SOR method with the SORW module. In addition to the module called MULT to compute the matrix-vector product, a routine called SRPASS was written to compute an SOR iteration. Note that the ITCOM5 common statement had to be included in the SRPASS routine so that the value of $\omega$ could be used. The program to set up the problem and the output resulting from this call to SORW is given below:

PROGRAM MAIN (OUTPUT,TAPE6=OUTPUT) C C ... ARRAY DECLARATIONS. C REAL RHS(100), U(100), WKSP(600), UBAR(1), RPARM(30) INTEGER IPARM(30) EXTERNAL MULT, SRPASS C NW = 600 C C ... GENERATE RHS. C NX = 10 NY = 10 N = NX*NY H = 1.0/FLOAT(NX + 1) DO 10 I = 1,N RHS(I) = 0.0 10 CONTINUE K = 0 DO 30 J = 1,NY Y = FLOAT(J)*H DO 25 I = 1,NX X = FLOAT(I)*H K = K + 1 IF (J .EQ. 1) RHS(K) = RHS(K) + 2.0 IF (J .EQ. NY) RHS(K) = RHS(K) + 2.0*(1.0 + X) IF (I .EQ. 1) RHS(K) = RHS(K) + 1.0 IF (I .EQ. NX) RHS(K) = RHS(K) + 1.0 + Y 25 CONTINUE 30 CONTINUE CALL DFAULT (IPARM,RPARM) C C ... NOW, RESET SOME DEFAULT VALUES. C IPARM(3) = 3 RPARM(1) = 1.0E-8 C C ... GENERATE AN INITIAL GUESS FOR U AND CALL NSPCG. C CALL VFILL (N,U,0.0) C CALL SORW (MULT,SRPASS,COEF,JCOEF,WFAC,JWFAC,N, A U,UBAR,RHS,WKSP,NW,IPARM,RPARM,IER) STOP END SUBROUTINE MULT (COEF,JCOEF,WFAC,JWFAC,N,X,Y) REAL X(N), Y(N) NX = 10 C C ... COMPUTE PRODUCT AS IF FIRST SUPERDIAGONAL AND FIRST C SUBDIAGONAL WERE FULL. C Y(1) = 6.0*X(1) - X(2) - 2.0*X(NX+1) DO 10 I = 2,NX Y(I) = 6.0*X(I) - X(I+1) - X(I-1) - 2.0*X(I+NX) 10 CONTINUE DO 15 I = NX+1,N-NX Y(I) = 6.0*X(I) - X(I+1) - X(I-1) - 2.0*(X(I+NX) + X(I-NX)) 15 CONTINUE DO 20 I = N-NX+1,N-1 Y(I) = 6.0*X(I) - X(I+1) - X(I-1) - 2.0*X(I-NX) 20 CONTINUE Y(N) = 6.0*X(N) - X(N-1) - 2.0*X(N-NX) C C ... MAKE CORRECTIONS TO Y VECTOR FOR ZEROS IN FIRST SUPERDIAGONAL C AND FIRST SUBDIAGONAL. C DO 25 I = NX,N-NX,NX Y(I) = Y(I) + X(I+1) 25 CONTINUE DO 30 I = NX+1,N-NX+1,NX Y(I) = Y(I) + X(I-1) 30 CONTINUE RETURN END SUBROUTINE SRPASS (COEF,JCOEF,WFAC,JWFAC,N,U,RHS,UNEW) C C ... SRPASS DOES AN SOR ITERATION. C C UNEW = INV((1/W)*D + L)*(((1/W)-1)*D*UN + RHS - U*UN) C C ... PARAMETERS -- C C N ORDER OF SYSTEM C U CURRENT SOLUTION VECTOR C RHS RIGHT HAND SIDE C UNEW UPDATED SOLUTION VECTOR C C ... SPECIFICATIONS FOR PARAMETERS C DIMENSION U(1), UNEW(1), RHS(1) LOGICAL OMGADP COMMON / ITCOM5 / OMGADP, OMEGA, ALPHAB, BETAB, FFF, SPECR C C ... TEMP = ((1/W)-1)*D*UN + RHS - U*UN C UNEW IS USED FOR TEMP. C NX = 10 CON = 6.0*(1.0/OMEGA - 1.0) DO 10 I = 1,N UNEW(I) = CON*U(I) + RHS(I) 10 CONTINUE DO 15 I = 1,N-1 UNEW(I) = UNEW(I) + U(I+1) 15 CONTINUE DO 20 I = 1,N-NX UNEW(I) = UNEW(I) + 2.0*U(I+NX) 20 CONTINUE DO 25 I = NX,N-NX,NX UNEW(I) = UNEW(I) - U(I+1) 25 CONTINUE C C ... UNEW = INV((1/W)*D + L)*TEMP C CON = OMEGA/6.0 DO 40 J = 1,NX IBGN = (J-1)*NX + 1 IEND = J*NX UNEW(IBGN) = CON*UNEW(IBGN) DO 30 I = IBGN+1,IEND UNEW(I) = CON*(UNEW(I) + UNEW(I-1)) 30 CONTINUE IF (J .EQ. NX) GO TO 40 DO 35 I = IBGN,IEND UNEW(I+NX) = UNEW(I+NX) + 2.0*UNEW(I) 35 CONTINUE 40 CONTINUE RETURN END
INITIAL ITERATIVE PARAMETERS GENERAL AND ACCELERATION PARAMETERS IPARM( 1) = 2 (NTEST ) IPARM( 2) = 100 (ITMAX ) IPARM( 3) = 3 (LEVEL ) IPARM( 4) = 6 (NOUT ) IPARM( 5) = 0 (IDGTS ) IPARM( 6) = 1 (MAXADP) IPARM( 7) = 1 (MINADP) IPARM( 8) = 1 (IOMGAD) IPARM( 9) = 5 (NS1 ) IPARM(10) = 100000 (NS2 ) IPARM(11) = 0 (NS3 ) RPARM( 1) = .10000000E-07 (ZETA ) RPARM( 2) = .20000000E+01 (EMAX ) RPARM( 3) = .10000000E+01 (EMIN ) RPARM( 4) = .75000000E+00 (FF ) RPARM( 5) = .75000000E+00 (FFF ) RPARM( 6) = .00000000E+00 (TIMIT ) RPARM( 7) = .00000000E+00 (DIGIT1) RPARM( 8) = .00000000E+00 (DIGIT2) RPARM( 9) = .10000000E+01 (OMEGA ) RPARM(10) = .00000000E+00 (ALPHAB) RPARM(11) = .25000000E+00 (BETAB ) RPARM(12) = .00000000E+00 (SPECR )
SOR INTERMEDIATE OUTPUT AFTER EACH ITERATION NUMBER OF CONVERGENCE EMAX OMEGA SPECTRAL ITERATIONS TEST RADIUS 0 0 .10000E+04 .20000E+01 .10000E+01 .00000E+00 1 0 .10000E+04 .20000E+01 .10000E+01 .53872E+00 2 0 .72510E+00 .20000E+01 .10000E+01 .70482E+00 3 0 .70622E+00 .20000E+01 .10000E+01 .78917E+00 4 0 .70165E+00 .91538E+00 .10000E+01 .83791E+00 5 1 .70165E+00 .91538E+00 .14259E+01 .18722E+01 6 1 .70165E+00 .91538E+00 .14259E+01 .81405E+00 7 1 .70165E+00 .91538E+00 .14259E+01 .82807E+00 8 1 .73605E+00 .91538E+00 .14259E+01 .83682E+00 9 1 .63941E+00 .91538E+00 .14259E+01 .84186E+00 10 1 .33339E+00 .91538E+00 .14259E+01 .84339E+00 11 1 .27588E+00 .91538E+00 .14259E+01 .84086E+00 12 1 .22191E+00 .91538E+00 .14259E+01 .83484E+00 13 1 .17539E+00 .91538E+00 .14259E+01 .82715E+00 14 1 .13764E+00 .91538E+00 .14259E+01 .81950E+00 15 1 .96272E-01 .91538E+00 .14259E+01 .81267E+00 16 1 .75684E-01 .91538E+00 .14259E+01 .80756E+00 17 1 .59578E-01 .91538E+00 .14259E+01 .80356E+00 18 1 .46828E-01 .91538E+00 .14259E+01 .80005E+00 19 1 .36758E-01 .91538E+00 .14259E+01 .79698E+00
46 2 .33557E-05 .95950E+00 .15604E+01 .59139E+00 47 2 .19459E-05 .95950E+00 .15604E+01 .58663E+00 48 2 .11480E-05 .95950E+00 .15604E+01 .58801E+00 49 2 .69725E-06 .95950E+00 .15604E+01 .59582E+00 50 2 .43310E-06 .95950E+00 .15604E+01 .60581E+00 51 2 .26430E-06 .95950E+00 .15604E+01 .60755E+00 52 2 .14596E-06 .95950E+00 .15604E+01 .58457E+00 53 2 .83253E-07 .95950E+00 .15604E+01 .57860E+00 54 2 .51143E-07 .95950E+00 .15604E+01 .59313E+00 55 2 .30889E-07 .95950E+00 .15604E+01 .59749E+00 56 2 .17899E-07 .95950E+00 .15604E+01 .59011E+00 57 2 .10123E-07 .95950E+00 .15604E+01 .57979E+00 58 2 .55955E-08 .95950E+00 .15604E+01 .56811E+00 SOR HAS CONVERGED IN 58 ITERATIONS
FINAL ITERATIVE PARAMETERS GENERAL AND ACCELERATION PARAMETERS IPARM( 1) = 2 (NTEST ) IPARM( 2) = 58 (ITMAX ) IPARM( 3) = 3 (LEVEL ) IPARM( 4) = 6 (NOUT ) IPARM( 5) = 0 (IDGTS ) IPARM( 6) = 1 (MAXADP) IPARM( 7) = 1 (MINADP) IPARM( 8) = 1 (IOMGAD) IPARM( 9) = 5 (NS1 ) IPARM(10) = 100000 (NS2 ) IPARM(11) = 0 (NS3 ) RPARM( 1) = .10000000E-07 (ZETA ) RPARM( 2) = .95950463E+00 (EMAX ) RPARM( 3) = .10000000E+01 (EMIN ) RPARM( 4) = .75000000E+00 (FF ) RPARM( 5) = .75000000E+00 (FFF ) RPARM( 6) = .62900000E+00 (TIMIT ) RPARM( 7) = .82521604E+01 (DIGIT1) RPARM( 8) = .87703729E+01 (DIGIT2) RPARM( 9) = .15604363E+01 (OMEGA ) RPARM(10) = .00000000E+00 (ALPHAB) RPARM(11) = .25000000E+00 (BETAB ) RPARM(12) = .56811199E+00 (SPECR )

Next: Installation Guide Up: NSPCG User's Guide Previous: Sample Usage