* F04QAF Example Program Text * Mark 14 Revised. NAG Copyright 1989. * .. Parameters .. INTEGER MMAX, NMAX, LIWORK, LRWORK PARAMETER (MMAX=100,NMAX=99,LIWORK=1,LRWORK=1) INTEGER NOUT PARAMETER (NOUT=6) * .. Scalars in Common .. INTEGER NCOLS, NROWS * .. Local Scalars .. DOUBLE PRECISION ACOND, ANORM, ARNORM, ATOL, BTOL, C, CONLIM, + DAMP, H, RNORM, XNORM INTEGER I, I1, IFAIL, INFORM, ITN, ITNLIM, K, M, MSGLVL, + N * .. Local Arrays .. DOUBLE PRECISION B(MMAX), RWORK(LRWORK), SE(NMAX), WORK(NMAX,2), + X(NMAX) INTEGER IWORK(LIWORK) * .. External Subroutines .. EXTERNAL APROD, F04QAF, X04ABF * .. Common blocks .. COMMON /USER/NROWS, NCOLS * .. Executable Statements .. WRITE (NOUT,*) 'F04QAF Example Program Results' WRITE (NOUT,*) CALL X04ABF(1,NOUT) NROWS = 4 NCOLS = 4 H = 0.1D0 N = NCOLS*NROWS - 4 M = N + 1 IF (NROWS.LT.3 .OR. NCOLS.LT.3 .OR. M.GT.MMAX) THEN WRITE (NOUT,99998) 'NROWS or NCOLS is out of range: NROWS = ', + NROWS, ' NCOLS = ', NCOLS ELSE * * Initialize RHS and other quantities required by F04QAF. * Convergence will be sooner if we do not regard A as exact, * so ATOL is not set to zero. * DO 20 I = 1, N B(I) = 0.0D0 20 CONTINUE C = -H**2 I1 = NROWS DO 60 K = 3, NCOLS DO 40 I = I1, I1 + NROWS - 3 B(I) = C 40 CONTINUE I1 = I1 + NROWS 60 CONTINUE B(M) = 1.0D0/H DAMP = 0.0D0 ATOL = 1.0D-5 BTOL = 1.0D-4 CONLIM = 1.0D0/ATOL ITNLIM = 100 * * Set MSGLVL to 2 to get output at each iteration * MSGLVL = 1 IFAIL = 1 * CALL F04QAF(M,N,B,X,SE,APROD,DAMP,ATOL,BTOL,CONLIM,ITNLIM, + MSGLVL,ITN,ANORM,ACOND,RNORM,ARNORM,XNORM,WORK, + RWORK,LRWORK,IWORK,LIWORK,INFORM,IFAIL) * WRITE (NOUT,*) IF (IFAIL.NE.0) THEN WRITE (NOUT,99999) 'F04QAF fails. IFAIL =', IFAIL ELSE WRITE (NOUT,*) 'Solution returned by F04QAF' WRITE (NOUT,99997) (X(I),I=1,N) WRITE (NOUT,*) WRITE (NOUT,99996) 'Norm of the residual = ', RNORM END IF END IF STOP * 99999 FORMAT (1X,A,I3) 99998 FORMAT (1X,A,I6,A,I6) 99997 FORMAT (1X,5F9.3) 99996 FORMAT (1X,A,1P,D12.2) END * SUBROUTINE APROD(MODE,M,N,X,Y,RWORK,LRWORK,IWORK,LIWORK) * APROD returns * Y = Y + A*X when MODE = 1 * X = X + ( A**T )*Y when MODE = 2 * for a given X and Y. * .. Scalar Arguments .. INTEGER LIWORK, LRWORK, M, MODE, N * .. Array Arguments .. DOUBLE PRECISION RWORK(LRWORK), X(N), Y(M) INTEGER IWORK(LIWORK) * .. Scalars in Common .. INTEGER NCOLS, NROWS * .. Local Scalars .. INTEGER J, J1, J2 * .. External Subroutines .. EXTERNAL ATIMES * .. Common blocks .. COMMON /USER/NROWS, NCOLS * .. Executable Statements .. IF (MODE.NE.2) THEN CALL ATIMES(NROWS,NCOLS,N,X,Y) DO 20 J = 1, NROWS - 2 Y(M) = Y(M) + X(J) 20 CONTINUE DO 40 J = 1, NCOLS - 2 Y(M) = Y(M) + X(J*NROWS-1) + X(J*NROWS+NROWS-2) 40 CONTINUE DO 60 J = M - NROWS + 2, N Y(M) = Y(M) + X(J) 60 CONTINUE ELSE CALL ATIMES(NROWS,NCOLS,N,Y,X) DO 80 J = 1, NROWS - 2 X(J) = X(J) + Y(M) 80 CONTINUE DO 100 J = 1, NCOLS - 2 J1 = J*NROWS - 1 J2 = J1 + NROWS - 1 X(J1) = X(J1) + Y(M) X(J2) = X(J2) + Y(M) 100 CONTINUE DO 120 J = M - NROWS + 2, N X(J) = X(J) + Y(M) 120 CONTINUE END IF RETURN END * SUBROUTINE ATIMES(NROWS,NCOLS,N,X,Y) * ATIMES is called by routine APROD and returns * Y = Y + ANN*X, * where ANN is the N by N symmetric part of the matrix A. * .. Scalar Arguments .. INTEGER N, NCOLS, NROWS * .. Array Arguments .. DOUBLE PRECISION X(N), Y(N) * .. Local Scalars .. INTEGER I, I1, I2, I3, IL, J * .. Executable Statements .. DO 20 J = 1, NROWS - 2 Y(J) = Y(J) + X(J) - X(J+NROWS-1) 20 CONTINUE DO 60 J = 1, NCOLS - 2 I = J*NROWS - 1 Y(I) = Y(I) + X(I) - X(I+1) I1 = I + 1 IL = I1 + NROWS - 3 DO 40 I = I1, IL I2 = I - NROWS IF (J.EQ.1) I2 = I2 + 1 I3 = I + NROWS IF (J.EQ.NCOLS-2) I3 = I3 - 1 Y(I) = Y(I) - X(I2) - X(I-1) + 4.0D0*X(I) - X(I+1) - X(I3) 40 CONTINUE I = IL + 1 Y(I) = Y(I) - X(I-1) + X(I) 60 CONTINUE DO 80 J = N - NROWS + 3, N Y(J) = Y(J) - X(J-NROWS+1) + X(J) 80 CONTINUE RETURN END