CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
dgmres.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dgmres (N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSOLVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, SB, SX, RGWK, LRGW, IGWK, LIGW, RWORK, IWORK)
 
double precision function dnrm2 (N, DX, INCX)
 
subroutine dpigmr (N, R0, SR, SZ, JSCAL, MAXL, MAXLP1, KMP, NRSTS, JPRE, MATVEC, MSOLVE, NMSL, Z, V, HES, Q, LGMR, RPAR, IPAR, WK, DL, RHOL, NRMAX, B, BNRM, X, XL, ITOL, TOL, NELT, IA, JA, A, ISYM, IUNIT, IFLAG, ERR)
 
subroutine drlcal (N, KMP, LL, MAXL, V, Q, RL, SNORMW, PROD, R0NRM)
 
subroutine daxpy (N, DA, DX, INCX, DY, INCY)
 
subroutine dcopy (N, DX, INCX, DY, INCY)
 
subroutine dhels (A, LDA, N, Q, B)
 
subroutine dheqr (A, LDA, N, Q, INFO, IJOB)
 
integer function isdgmr (N, B, X, XL, NELT, IA, JA, A, ISYM, MSOLVE, NMSL, ITOL, TOL, ITMAX, ITER, ERR, IUNIT, R, Z, DZ, RWORK, IWORK, RNRM, BNRM, SB, SX, JSCAL, KMP, LGMR, MAXL, MAXLP1, V, Q, SNORMW, PROD, R0NRM, HES, JPRE)
 
subroutine dorth (VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
 
subroutine dxlcal (N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM, WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR, NELT, IA, JA, A, ISYM)
 
double precision function ddot (N, DX, INCX, DY, INCY)
 

Function/Subroutine Documentation

◆ daxpy()

subroutine daxpy (   N,
double precision  DA,
double precision, dimension(*)  DX,
  INCX,
double precision, dimension(*)  DY,
  INCY 
)
1276  use omp_lib
1277 C***BEGIN PROLOGUE DAXPY
1278 C***PURPOSE Compute a constant times a vector plus a vector.
1279 C***LIBRARY SLATEC (BLAS)
1280 C***CATEGORY D1A7
1281 C***TYPE DOUBLE PRECISION (SAXPY-S, DAXPY-D, CAXPY-C)
1282 C***KEYWORDS BLAS, LINEAR ALGEBRA, TRIAD, VECTOR
1283 C***AUTHOR Lawson, C. L., (JPL)
1284 C Hanson, R. J., (SNLA)
1285 C Kincaid, D. R., (U. of Texas)
1286 C Krogh, F. T., (JPL)
1287 C***DESCRIPTION
1288 C
1289 C B L A S Subprogram
1290 C Description of Parameters
1291 C
1292 C --Input--
1293 C N number of elements in input vector(s)
1294 C DA double precision scalar multiplier
1295 C DX double precision vector with N elements
1296 C INCX storage spacing between elements of DX
1297 C DY double precision vector with N elements
1298 C INCY storage spacing between elements of DY
1299 C
1300 C --Output--
1301 C DY double precision result (unchanged if N .LE. 0)
1302 C
1303 C Overwrite double precision DY with double precision DA*DX + DY.
1304 C For I = 0 to N-1, replace DY(LY+I*INCY) with DA*DX(LX+I*INCX) +
1305 C DY(LY+I*INCY),
1306 C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
1307 C defined in a similar way using INCY.
1308 C
1309 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
1310 C Krogh, Basic linear algebra subprograms for Fortran
1311 C usage, Algorithm No. 539, Transactions on Mathematical
1312 C Software 5, 3 (September 1979), pp. 308-323.
1313 C***ROUTINES CALLED (NONE)
1314 C***REVISION HISTORY (YYMMDD)
1315 C 791001 DATE WRITTEN
1316 C 890831 Modified array declarations. (WRB)
1317 C 890831 REVISION DATE from Version 3.2
1318 C 891214 Prologue converted to Version 4.0 format. (BAB)
1319 C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
1320 C 920501 Reformatted the REFERENCES section. (WRB)
1321 C***END PROLOGUE DAXPY
1322  DOUBLE PRECISION dx(*), dy(*), da
1323 C***FIRST EXECUTABLE STATEMENT DAXPY
1324  IF (n.LE.0 .OR. da.EQ.0.0d0) RETURN
1325  IF (incx .EQ. incy) IF (incx-1) 5,20,60
1326 C
1327 C Code for unequal or nonpositive increments.
1328 C
1329  5 ix = 1
1330  iy = 1
1331  IF (incx .LT. 0) ix = (-n+1)*incx + 1
1332  IF (incy .LT. 0) iy = (-n+1)*incy + 1
1333  DO 10 i = 1,n
1334  dy(iy) = dy(iy) + da*dx(ix)
1335  ix = ix + incx
1336  iy = iy + incy
1337  10 CONTINUE
1338  RETURN
1339 C
1340 C Code for both increments equal to 1.
1341 C
1342 C Clean-up loop so remaining vector length is a multiple of 4.
1343 C
1344  20 m = mod(n,4)
1345  IF (m .EQ. 0) GO TO 40
1346 c$omp parallel default(none)
1347 c$omp& shared(m,dy,da,dx)
1348 c$omp& private(i)
1349 c$omp do
1350  DO 30 i = 1,m
1351  dy(i) = dy(i) + da*dx(i)
1352  30 CONTINUE
1353 c$omp end do
1354 c$omp end parallel
1355  IF (n .LT. 4) RETURN
1356  40 mp1 = m + 1
1357 c$omp parallel default(none)
1358 c$omp& shared(mp1,n,dy,da,dx)
1359 c$omp& private(i)
1360 c$omp do
1361  DO 50 i = mp1,n,4
1362  dy(i) = dy(i) + da*dx(i)
1363  dy(i+1) = dy(i+1) + da*dx(i+1)
1364  dy(i+2) = dy(i+2) + da*dx(i+2)
1365  dy(i+3) = dy(i+3) + da*dx(i+3)
1366  50 CONTINUE
1367 c$omp end do
1368 c$omp end parallel
1369  RETURN
1370 C
1371 C Code for equal, positive, non-unit increments.
1372 C
1373  60 ns = n*incx
1374  DO 70 i = 1,ns,incx
1375  dy(i) = da*dx(i) + dy(i)
1376  70 CONTINUE
1377  RETURN

◆ dcopy()

subroutine dcopy (   N,
double precision, dimension(*)  DX,
  INCX,
double precision, dimension(*)  DY,
  INCY 
)
1381  use omp_lib
1382 C***BEGIN PROLOGUE DCOPY
1383 C***PURPOSE Copy a vector.
1384 C***LIBRARY SLATEC (BLAS)
1385 C***CATEGORY D1A5
1386 C***TYPE DOUBLE PRECISION (SCOPY-S, DCOPY-D, CCOPY-C, ICOPY-I)
1387 C***KEYWORDS BLAS, COPY, LINEAR ALGEBRA, VECTOR
1388 C***AUTHOR Lawson, C. L., (JPL)
1389 C Hanson, R. J., (SNLA)
1390 C Kincaid, D. R., (U. of Texas)
1391 C Krogh, F. T., (JPL)
1392 C***DESCRIPTION
1393 C
1394 C B L A S Subprogram
1395 C Description of Parameters
1396 C
1397 C --Input--
1398 C N number of elements in input vector(s)
1399 C DX double precision vector with N elements
1400 C INCX storage spacing between elements of DX
1401 C DY double precision vector with N elements
1402 C INCY storage spacing between elements of DY
1403 C
1404 C --Output--
1405 C DY copy of vector DX (unchanged if N .LE. 0)
1406 C
1407 C Copy double precision DX to double precision DY.
1408 C For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY),
1409 C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
1410 C defined in a similar way using INCY.
1411 C
1412 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
1413 C Krogh, Basic linear algebra subprograms for Fortran
1414 C usage, Algorithm No. 539, Transactions on Mathematical
1415 C Software 5, 3 (September 1979), pp. 308-323.
1416 C***ROUTINES CALLED (NONE)
1417 C***REVISION HISTORY (YYMMDD)
1418 C 791001 DATE WRITTEN
1419 C 890831 Modified array declarations. (WRB)
1420 C 890831 REVISION DATE from Version 3.2
1421 C 891214 Prologue converted to Version 4.0 format. (BAB)
1422 C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
1423 C 920501 Reformatted the REFERENCES section. (WRB)
1424 C***END PROLOGUE DCOPY
1425  DOUBLE PRECISION dx(*), dy(*)
1426 C***FIRST EXECUTABLE STATEMENT DCOPY
1427  IF (n .LE. 0) RETURN
1428  IF (incx .EQ. incy) IF (incx-1) 5,20,60
1429 C
1430 C Code for unequal or nonpositive increments.
1431 C
1432  5 ix = 1
1433  iy = 1
1434  IF (incx .LT. 0) ix = (-n+1)*incx + 1
1435  IF (incy .LT. 0) iy = (-n+1)*incy + 1
1436  DO 10 i = 1,n
1437  dy(iy) = dx(ix)
1438  ix = ix + incx
1439  iy = iy + incy
1440  10 CONTINUE
1441  RETURN
1442 C
1443 C Code for both increments equal to 1.
1444 C
1445 C Clean-up loop so remaining vector length is a multiple of 7.
1446 C
1447  20 m = mod(n,7)
1448  IF (m .EQ. 0) GO TO 40
1449  DO 30 i = 1,m
1450  dy(i) = dx(i)
1451  30 CONTINUE
1452  IF (n .LT. 7) RETURN
1453  40 mp1 = m + 1
1454 c$omp parallel default(none)
1455 c$omp& shared(mp1,n,dx,dy)
1456 c$omp& private(i)
1457 c$omp do
1458  DO 50 i = mp1,n,7
1459  dy(i) = dx(i)
1460  dy(i+1) = dx(i+1)
1461  dy(i+2) = dx(i+2)
1462  dy(i+3) = dx(i+3)
1463  dy(i+4) = dx(i+4)
1464  dy(i+5) = dx(i+5)
1465  dy(i+6) = dx(i+6)
1466  50 CONTINUE
1467 c$omp end do
1468 c$omp end parallel
1469  RETURN
1470 C
1471 C Code for equal, positive, non-unit increments.
1472 C
1473  60 ns = n*incx
1474  DO 70 i = 1,ns,incx
1475  dy(i) = dx(i)
1476  70 CONTINUE
1477  RETURN

◆ ddot()

double precision function ddot (   N,
double precision, dimension(*)  DX,
  INCX,
double precision, dimension(*)  DY,
  INCY 
)
2469  use omp_lib
2470 C***BEGIN PROLOGUE DDOT
2471 C***PURPOSE Compute the inner product of two vectors.
2472 C***LIBRARY SLATEC (BLAS)
2473 C***CATEGORY D1A4
2474 C***TYPE DOUBLE PRECISION (SDOT-S, DDOT-D, CDOTU-C)
2475 C***KEYWORDS BLAS, INNER PRODUCT, LINEAR ALGEBRA, VECTOR
2476 C***AUTHOR Lawson, C. L., (JPL)
2477 C Hanson, R. J., (SNLA)
2478 C Kincaid, D. R., (U. of Texas)
2479 C Krogh, F. T., (JPL)
2480 C***DESCRIPTION
2481 C
2482 C B L A S Subprogram
2483 C Description of Parameters
2484 C
2485 C --Input--
2486 C N number of elements in input vector(s)
2487 C DX double precision vector with N elements
2488 C INCX storage spacing between elements of DX
2489 C DY double precision vector with N elements
2490 C INCY storage spacing between elements of DY
2491 C
2492 C --Output--
2493 C DDOT double precision dot product (zero if N .LE. 0)
2494 C
2495 C Returns the dot product of double precision DX and DY.
2496 C DDOT = sum for I = 0 to N-1 of DX(LX+I*INCX) * DY(LY+I*INCY),
2497 C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
2498 C defined in a similar way using INCY.
2499 C
2500 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
2501 C Krogh, Basic linear algebra subprograms for Fortran
2502 C usage, Algorithm No. 539, Transactions on Mathematical
2503 C Software 5, 3 (September 1979), pp. 308-323.
2504 C***ROUTINES CALLED (NONE)
2505 C***REVISION HISTORY (YYMMDD)
2506 C 791001 DATE WRITTEN
2507 C 890831 Modified array declarations. (WRB)
2508 C 890831 REVISION DATE from Version 3.2
2509 C 891214 Prologue converted to Version 4.0 format. (BAB)
2510 C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
2511 C 920501 Reformatted the REFERENCES section. (WRB)
2512 C***END PROLOGUE DDOT
2513  DOUBLE PRECISION dx(*), dy(*)
2514 C***FIRST EXECUTABLE STATEMENT DDOT
2515  ddot = 0.0d0
2516  IF (n .LE. 0) RETURN
2517  IF (incx .EQ. incy) IF (incx-1) 5,20,60
2518 C
2519 C Code for unequal or nonpositive increments.
2520 C
2521  5 ix = 1
2522  iy = 1
2523  IF (incx .LT. 0) ix = (-n+1)*incx + 1
2524  IF (incy .LT. 0) iy = (-n+1)*incy + 1
2525  DO 10 i = 1,n
2526  ddot = ddot + dx(ix)*dy(iy)
2527  ix = ix + incx
2528  iy = iy + incy
2529  10 CONTINUE
2530  RETURN
2531 C
2532 C Code for both increments equal to 1.
2533 C
2534 C Clean-up loop so remaining vector length is a multiple of 5.
2535 C
2536  20 m = mod(n,5)
2537  IF (m .EQ. 0) GO TO 40
2538  DO 30 i = 1,m
2539  ddot = ddot + dx(i)*dy(i)
2540  30 CONTINUE
2541  IF (n .LT. 5) RETURN
2542  40 mp1 = m + 1
2543 c$omp parallel default(none)
2544 c$omp& shared(mp1,n,dx,dy,ddot)
2545 c$omp& private(i)
2546 c$omp do reduction(+:ddot)
2547  DO 50 i = mp1,n,5
2548  ddot = ddot + dx(i)*dy(i) + dx(i+1)*dy(i+1) + dx(i+2)*dy(i+2) +
2549  1 dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4)
2550  50 CONTINUE
2551 c$omp end do
2552 c$omp end parallel
2553  RETURN
2554 C
2555 C Code for equal, positive, non-unit increments.
2556 C
2557  60 ns = n*incx
2558  DO 70 i = 1,ns,incx
2559  ddot = ddot + dx(i)*dy(i)
2560  70 CONTINUE
2561  RETURN
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
Definition: dgmres.f:2469

◆ dgmres()

subroutine dgmres ( integer  N,
double precision, dimension(n)  B,
double precision, dimension(n)  X,
integer  NELT,
integer, dimension(nelt)  IA,
integer, dimension(nelt)  JA,
double precision, dimension(nelt)  A,
integer  ISYM,
external  MATVEC,
external  MSOLVE,
integer  ITOL,
double precision  TOL,
integer  ITMAX,
integer  ITER,
double precision  ERR,
integer  IERR,
integer  IUNIT,
double precision, dimension(n)  SB,
double precision, dimension(n)  SX,
double precision, dimension(lrgw)  RGWK,
integer  LRGW,
integer, dimension(ligw)  IGWK,
integer  LIGW,
double precision, dimension(*)  RWORK,
integer, dimension(*)  IWORK 
)
8 C***BEGIN PROLOGUE DGMRES
9 C***PURPOSE Preconditioned GMRES iterative sparse Ax=b solver.
10 C This routine uses the generalized minimum residual
11 C (GMRES) method with preconditioning to solve
12 C non-symmetric linear systems of the form: Ax = b.
13 C***LIBRARY SLATEC (SLAP)
14 C***CATEGORY D2A4, D2B4
15 C***TYPE DOUBLE PRECISION (SGMRES-S, DGMRES-D)
16 C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
17 C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
18 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
19 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
20 C Seager, Mark K., (LLNL), seager@llnl.gov
21 C Lawrence Livermore National Laboratory
22 C PO Box 808, L-60
23 C Livermore, CA 94550 (510) 423-3141
24 C***DESCRIPTION
25 C
26 C *Usage:
27 C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
28 C INTEGER ITER, IERR, IUNIT, LRGW, IGWK(LIGW), LIGW
29 C INTEGER IWORK(USER DEFINED)
30 C DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, SB(N), SX(N)
31 C DOUBLE PRECISION RGWK(LRGW), RWORK(USER DEFINED)
32 C EXTERNAL MATVEC, MSOLVE
33 C
34 C CALL DGMRES(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSOLVE,
35 C $ ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, SB, SX,
36 C $ RGWK, LRGW, IGWK, LIGW, RWORK, IWORK)
37 C
38 C *Arguments:
39 C N :IN Integer.
40 C Order of the Matrix.
41 C B :IN Double Precision B(N).
42 C Right-hand side vector.
43 C X :INOUT Double Precision X(N).
44 C On input X is your initial guess for the solution vector.
45 C On output X is the final approximate solution.
46 C NELT :IN Integer.
47 C Number of Non-Zeros stored in A.
48 C IA :IN Integer IA(NELT).
49 C JA :IN Integer JA(NELT).
50 C A :IN Double Precision A(NELT).
51 C These arrays contain the matrix data structure for A.
52 C It could take any form. See "Description", below,
53 C for more details.
54 C ISYM :IN Integer.
55 C Flag to indicate symmetric storage format.
56 C If ISYM=0, all non-zero entries of the matrix are stored.
57 C If ISYM=1, the matrix is symmetric, and only the upper
58 C or lower triangle of the matrix is stored.
59 C MATVEC :EXT External.
60 C Name of a routine which performs the matrix vector multiply
61 C Y = A*X given A and X. The name of the MATVEC routine must
62 C be declared external in the calling program. The calling
63 C sequence to MATVEC is:
64 C CALL MATVEC(N, X, Y, NELT, IA, JA, A, ISYM)
65 C where N is the number of unknowns, Y is the product A*X
66 C upon return, X is an input vector, and NELT is the number of
67 C non-zeros in the SLAP IA, JA, A storage for the matrix A.
68 C ISYM is a flag which, if non-zero, denotes that A is
69 C symmetric and only the lower or upper triangle is stored.
70 C MSOLVE :EXT External.
71 C Name of the routine which solves a linear system Mz = r for
72 C z given r with the preconditioning matrix M (M is supplied via
73 C RWORK and IWORK arrays. The name of the MSOLVE routine must
74 C be declared external in the calling program. The calling
75 C sequence to MSOLVE is:
76 C CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
77 C Where N is the number of unknowns, R is the right-hand side
78 C vector and Z is the solution upon return. NELT, IA, JA, A and
79 C ISYM are defined as above. RWORK is a double precision array
80 C that can be used to pass necessary preconditioning information
81 C and/or workspace to MSOLVE. IWORK is an integer work array
82 C for the same purpose as RWORK.
83 C ITOL :IN Integer.
84 C Flag to indicate the type of convergence criterion used.
85 C ITOL=0 Means the iteration stops when the test described
86 C below on the residual RL is satisfied. This is
87 C the "Natural Stopping Criteria" for this routine.
88 C Other values of ITOL cause extra, otherwise
89 C unnecessary, computation per iteration and are
90 C therefore much less efficient. See ISDGMR (the
91 C stop test routine) for more information.
92 C ITOL=1 Means the iteration stops when the first test
93 C described below on the residual RL is satisfied,
94 C and there is either right or no preconditioning
95 C being used.
96 C ITOL=2 Implies that the user is using left
97 C preconditioning, and the second stopping criterion
98 C below is used.
99 C ITOL=3 Means the iteration stops when the third test
100 C described below on Minv*Residual is satisfied, and
101 C there is either left or no preconditioning being
102 C used.
103 C ITOL=11 is often useful for checking and comparing
104 C different routines. For this case, the user must
105 C supply the "exact" solution or a very accurate
106 C approximation (one with an error much less than
107 C TOL) through a common block,
108 C COMMON /DSLBLK/ SOLN( )
109 C If ITOL=11, iteration stops when the 2-norm of the
110 C difference between the iterative approximation and
111 C the user-supplied solution divided by the 2-norm
112 C of the user-supplied solution is less than TOL.
113 C Note that this requires the user to set up the
114 C "COMMON /DSLBLK/ SOLN(LENGTH)" in the calling
115 C routine. The routine with this declaration should
116 C be loaded before the stop test so that the correct
117 C length is used by the loader. This procedure is
118 C not standard Fortran and may not work correctly on
119 C your system (although it has worked on every
120 C system the authors have tried). If ITOL is not 11
121 C then this common block is indeed standard Fortran.
122 C TOL :INOUT Double Precision.
123 C Convergence criterion, as described below. If TOL is set
124 C to zero on input, then a default value of 500*(the smallest
125 C positive magnitude, machine epsilon) is used.
126 C ITMAX :DUMMY Integer.
127 C Maximum number of iterations in most SLAP routines. In
128 C this routine this does not make sense. The maximum number
129 C of iterations here is given by ITMAX = MAXL*(NRMAX+1).
130 C See IGWK for definitions of MAXL and NRMAX.
131 C ITER :OUT Integer.
132 C Number of iterations required to reach convergence, or
133 C ITMAX if convergence criterion could not be achieved in
134 C ITMAX iterations.
135 C ERR :OUT Double Precision.
136 C Error estimate of error in final approximate solution, as
137 C defined by ITOL. Letting norm() denote the Euclidean
138 C norm, ERR is defined as follows..
139 C
140 C If ITOL=0, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
141 C for right or no preconditioning, and
142 C ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
143 C norm(SB*(M-inverse)*B),
144 C for left preconditioning.
145 C If ITOL=1, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
146 C since right or no preconditioning
147 C being used.
148 C If ITOL=2, then ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
149 C norm(SB*(M-inverse)*B),
150 C since left preconditioning is being
151 C used.
152 C If ITOL=3, then ERR = Max |(Minv*(B-A*X(L)))(i)/x(i)|
153 C i=1,n
154 C If ITOL=11, then ERR = norm(SB*(X(L)-SOLN))/norm(SB*SOLN).
155 C IERR :OUT Integer.
156 C Return error flag.
157 C IERR = 0 => All went well.
158 C IERR = 1 => Insufficient storage allocated for
159 C RGWK or IGWK.
160 C IERR = 2 => Routine DGMRES failed to reduce the norm
161 C of the current residual on its last call,
162 C and so the iteration has stalled. In
163 C this case, X equals the last computed
164 C approximation. The user must either
165 C increase MAXL, or choose a different
166 C initial guess.
167 C IERR =-1 => Insufficient length for RGWK array.
168 C IGWK(6) contains the required minimum
169 C length of the RGWK array.
170 C IERR =-2 => Illegal value of ITOL, or ITOL and JPRE
171 C values are inconsistent.
172 C For IERR <= 2, RGWK(1) = RHOL, which is the norm on the
173 C left-hand-side of the relevant stopping test defined
174 C below associated with the residual for the current
175 C approximation X(L).
176 C IUNIT :IN Integer.
177 C Unit number on which to write the error at each iteration,
178 C if this is desired for monitoring convergence. If unit
179 C number is 0, no writing will occur.
180 C SB :IN Double Precision SB(N).
181 C Array of length N containing scale factors for the right
182 C hand side vector B. If JSCAL.eq.0 (see below), SB need
183 C not be supplied.
184 C SX :IN Double Precision SX(N).
185 C Array of length N containing scale factors for the solution
186 C vector X. If JSCAL.eq.0 (see below), SX need not be
187 C supplied. SB and SX can be the same array in the calling
188 C program if desired.
189 C RGWK :INOUT Double Precision RGWK(LRGW).
190 C Double Precision array used for workspace by DGMRES.
191 C On return, RGWK(1) = RHOL. See IERR for definition of RHOL.
192 C LRGW :IN Integer.
193 C Length of the double precision workspace, RGWK.
194 C LRGW >= 1 + N*(MAXL+6) + MAXL*(MAXL+3).
195 C See below for definition of MAXL.
196 C For the default values, RGWK has size at least 131 + 16*N.
197 C IGWK :INOUT Integer IGWK(LIGW).
198 C The following IGWK parameters should be set by the user
199 C before calling this routine.
200 C IGWK(1) = MAXL. Maximum dimension of Krylov subspace in
201 C which X - X0 is to be found (where, X0 is the initial
202 C guess). The default value of MAXL is 10.
203 C IGWK(2) = KMP. Maximum number of previous Krylov basis
204 C vectors to which each new basis vector is made orthogonal.
205 C The default value of KMP is MAXL.
206 C IGWK(3) = JSCAL. Flag indicating whether the scaling
207 C arrays SB and SX are to be used.
208 C JSCAL = 0 => SB and SX are not used and the algorithm
209 C will perform as if all SB(I) = 1 and SX(I) = 1.
210 C JSCAL = 1 => Only SX is used, and the algorithm
211 C performs as if all SB(I) = 1.
212 C JSCAL = 2 => Only SB is used, and the algorithm
213 C performs as if all SX(I) = 1.
214 C JSCAL = 3 => Both SB and SX are used.
215 C IGWK(4) = JPRE. Flag indicating whether preconditioning
216 C is being used.
217 C JPRE = 0 => There is no preconditioning.
218 C JPRE > 0 => There is preconditioning on the right
219 C only, and the solver will call routine MSOLVE.
220 C JPRE < 0 => There is preconditioning on the left
221 C only, and the solver will call routine MSOLVE.
222 C IGWK(5) = NRMAX. Maximum number of restarts of the
223 C Krylov iteration. The default value of NRMAX = 10.
224 C if IWORK(5) = -1, then no restarts are performed (in
225 C this case, NRMAX is set to zero internally).
226 C The following IWORK parameters are diagnostic information
227 C made available to the user after this routine completes.
228 C IGWK(6) = MLWK. Required minimum length of RGWK array.
229 C IGWK(7) = NMS. The total number of calls to MSOLVE.
230 C LIGW :IN Integer.
231 C Length of the integer workspace, IGWK. LIGW >= 20.
232 C RWORK :WORK Double Precision RWORK(USER DEFINED).
233 C Double Precision array that can be used for workspace in
234 C MSOLVE.
235 C IWORK :WORK Integer IWORK(USER DEFINED).
236 C Integer array that can be used for workspace in MSOLVE.
237 C
238 C *Description:
239 C DGMRES solves a linear system A*X = B rewritten in the form:
240 C
241 C (SB*A*(M-inverse)*(SX-inverse))*(SX*M*X) = SB*B,
242 C
243 C with right preconditioning, or
244 C
245 C (SB*(M-inverse)*A*(SX-inverse))*(SX*X) = SB*(M-inverse)*B,
246 C
247 C with left preconditioning, where A is an N-by-N double precision
248 C matrix, X and B are N-vectors, SB and SX are diagonal scaling
249 C matrices, and M is a preconditioning matrix. It uses
250 C preconditioned Krylov subpace methods based on the
251 C generalized minimum residual method (GMRES). This routine
252 C optionally performs either the full orthogonalization
253 C version of the GMRES algorithm or an incomplete variant of
254 C it. Both versions use restarting of the linear iteration by
255 C default, although the user can disable this feature.
256 C
257 C The GMRES algorithm generates a sequence of approximations
258 C X(L) to the true solution of the above linear system. The
259 C convergence criteria for stopping the iteration is based on
260 C the size of the scaled norm of the residual R(L) = B -
261 C A*X(L). The actual stopping test is either:
262 C
263 C norm(SB*(B-A*X(L))) .le. TOL*norm(SB*B),
264 C
265 C for right preconditioning, or
266 C
267 C norm(SB*(M-inverse)*(B-A*X(L))) .le.
268 C TOL*norm(SB*(M-inverse)*B),
269 C
270 C for left preconditioning, where norm() denotes the Euclidean
271 C norm, and TOL is a positive scalar less than one input by
272 C the user. If TOL equals zero when DGMRES is called, then a
273 C default value of 500*(the smallest positive magnitude,
274 C machine epsilon) is used. If the scaling arrays SB and SX
275 C are used, then ideally they should be chosen so that the
276 C vectors SX*X(or SX*M*X) and SB*B have all their components
277 C approximately equal to one in magnitude. If one wants to
278 C use the same scaling in X and B, then SB and SX can be the
279 C same array in the calling program.
280 C
281 C The following is a list of the other routines and their
282 C functions used by DGMRES:
283 C DPIGMR Contains the main iteration loop for GMRES.
284 C DORTH Orthogonalizes a new vector against older basis vectors.
285 C DHEQR Computes a QR decomposition of a Hessenberg matrix.
286 C DHELS Solves a Hessenberg least-squares system, using QR
287 C factors.
288 C DRLCAL Computes the scaled residual RL.
289 C DXLCAL Computes the solution XL.
290 C ISDGMR User-replaceable stopping routine.
291 C
292 C This routine does not care what matrix data structure is
293 C used for A and M. It simply calls the MATVEC and MSOLVE
294 C routines, with the arguments as described above. The user
295 C could write any type of structure and the appropriate MATVEC
296 C and MSOLVE routines. It is assumed that A is stored in the
297 C IA, JA, A arrays in some fashion and that M (or INV(M)) is
298 C stored in IWORK and RWORK in some fashion. The SLAP
299 C routines DSDCG and DSICCG are examples of this procedure.
300 C
301 C Two examples of matrix data structures are the: 1) SLAP
302 C Triad format and 2) SLAP Column format.
303 C
304 C =================== S L A P Triad format ===================
305 C This routine requires that the matrix A be stored in the
306 C SLAP Triad format. In this format only the non-zeros are
307 C stored. They may appear in *ANY* order. The user supplies
308 C three arrays of length NELT, where NELT is the number of
309 C non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For
310 C each non-zero the user puts the row and column index of that
311 C matrix element in the IA and JA arrays. The value of the
312 C non-zero matrix element is placed in the corresponding
313 C location of the A array. This is an extremely easy data
314 C structure to generate. On the other hand it is not too
315 C efficient on vector computers for the iterative solution of
316 C linear systems. Hence, SLAP changes this input data
317 C structure to the SLAP Column format for the iteration (but
318 C does not change it back).
319 C
320 C Here is an example of the SLAP Triad storage format for a
321 C 5x5 Matrix. Recall that the entries may appear in any order.
322 C
323 C 5x5 Matrix SLAP Triad format for 5x5 matrix on left.
324 C 1 2 3 4 5 6 7 8 9 10 11
325 C |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21
326 C |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2
327 C | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1
328 C | 0 0 0 44 0|
329 C |51 0 53 0 55|
330 C
331 C =================== S L A P Column format ==================
332 C
333 C This routine requires that the matrix A be stored in the
334 C SLAP Column format. In this format the non-zeros are stored
335 C counting down columns (except for the diagonal entry, which
336 C must appear first in each "column") and are stored in the
337 C double precision array A. In other words, for each column
338 C in the matrix put the diagonal entry in A. Then put in the
339 C other non-zero elements going down the column (except the
340 C diagonal) in order. The IA array holds the row index for
341 C each non-zero. The JA array holds the offsets into the IA,
342 C A arrays for the beginning of each column. That is,
343 C IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the
344 C ICOL-th column in IA and A. IA(JA(ICOL+1)-1),
345 C A(JA(ICOL+1)-1) points to the end of the ICOL-th column.
346 C Note that we always have JA(N+1) = NELT+1, where N is the
347 C number of columns in the matrix and NELT is the number of
348 C non-zeros in the matrix.
349 C
350 C Here is an example of the SLAP Column storage format for a
351 C 5x5 Matrix (in the A and IA arrays '|' denotes the end of a
352 C column):
353 C
354 C 5x5 Matrix SLAP Column format for 5x5 matrix on left.
355 C 1 2 3 4 5 6 7 8 9 10 11
356 C |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
357 C |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3
358 C | 0 0 33 0 35| JA: 1 4 6 8 9 12
359 C | 0 0 0 44 0|
360 C |51 0 53 0 55|
361 C
362 C *Cautions:
363 C This routine will attempt to write to the Fortran logical output
364 C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that
365 C this logical unit is attached to a file or terminal before calling
366 C this routine with a non-zero value for IUNIT. This routine does
367 C not check for the validity of a non-zero IUNIT unit number.
368 C
369 C***REFERENCES 1. Peter N. Brown and A. C. Hindmarsh, Reduced Storage
370 C Matrix Methods in Stiff ODE Systems, Lawrence Liver-
371 C more National Laboratory Report UCRL-95088, Rev. 1,
372 C Livermore, California, June 1987.
373 C 2. Mark K. Seager, A SLAP for the Masses, in
374 C G. F. Carey, Ed., Parallel Supercomputing: Methods,
375 C Algorithms and Applications, Wiley, 1989, pp.135-155.
376 C***ROUTINES CALLED D1MACH, DCOPY, DNRM2, DPIGMR
377 C***REVISION HISTORY (YYMMDD)
378 C 890404 DATE WRITTEN
379 C 890404 Previous REVISION DATE
380 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
381 C 890922 Numerous changes to prologue to make closer to SLATEC
382 C standard. (FNF)
383 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
384 C 891004 Added new reference.
385 C 910411 Prologue converted to Version 4.0 format. (BAB)
386 C 910506 Corrected errors in C***ROUTINES CALLED list. (FNF)
387 C 920407 COMMON BLOCK renamed DSLBLK. (WRB)
388 C 920511 Added complete declaration section. (WRB)
389 C 920929 Corrected format of references. (FNF)
390 C 921019 Changed 500.0 to 500 to reduce SP/DP differences. (FNF)
391 C 921026 Added check for valid value of ITOL. (FNF)
392 C***END PROLOGUE DGMRES
393 C The following is for optimized compilation on LLNL/LTSS Crays.
394 CLLL. OPTIMIZE
395 C .. Scalar Arguments ..
396  DOUBLE PRECISION err, tol
397  INTEGER ierr, isym, iter, itmax, itol, iunit, ligw, lrgw, n, nelt
398 C .. Array Arguments ..
399  DOUBLE PRECISION a(nelt), b(n), rgwk(lrgw), rwork(*), sb(n),
400  + sx(n), x(n)
401  INTEGER ia(nelt), igwk(ligw), iwork(*), ja(nelt)
402 C .. Subroutine Arguments ..
403  EXTERNAL matvec, msolve
404 C .. Local Scalars ..
405  DOUBLE PRECISION bnrm, rhol, sum
406  INTEGER i, iflag, jpre, jscal, kmp, ldl, lgmr, lhes, lq, lr, lv,
407  + lw, lxl, lz, lzm1, maxl, maxlp1, nms, nmsl, nrmax, nrsts
408 C .. External Functions ..
409  DOUBLE PRECISION d1mach, dnrm2
410  EXTERNAL d1mach, dnrm2
411 C .. External Subroutines ..
412  EXTERNAL dcopy, dpigmr
413 C .. Intrinsic Functions ..
414  INTRINSIC sqrt
415 C***FIRST EXECUTABLE STATEMENT DGMRES
416  ierr = 0
417 C ------------------------------------------------------------------
418 C Load method parameters with user values or defaults.
419 C ------------------------------------------------------------------
420  maxl = igwk(1)
421  IF (maxl .EQ. 0) maxl = 10
422  IF (maxl .GT. n) maxl = n
423  kmp = igwk(2)
424  IF (kmp .EQ. 0) kmp = maxl
425  IF (kmp .GT. maxl) kmp = maxl
426  jscal = igwk(3)
427  jpre = igwk(4)
428 C Check for valid value of ITOL.
429  IF( (itol.LT.0) .OR. ((itol.GT.3).AND.(itol.NE.11)) ) GOTO 650
430 C Check for consistent values of ITOL and JPRE.
431  IF( itol.EQ.1 .AND. jpre.LT.0 ) GOTO 650
432  IF( itol.EQ.2 .AND. jpre.GE.0 ) GOTO 650
433  nrmax = igwk(5)
434  IF( nrmax.EQ.0 ) nrmax = 10
435 C If NRMAX .eq. -1, then set NRMAX = 0 to turn off restarting.
436  IF( nrmax.EQ.-1 ) nrmax = 0
437 C If input value of TOL is zero, set it to its default value.
438  IF( tol.EQ.0.0d0 ) tol = 500*d1mach(3)
439 C
440 C Initialize counters.
441  iter = 0
442  nms = 0
443  nrsts = 0
444 C ------------------------------------------------------------------
445 C Form work array segment pointers.
446 C ------------------------------------------------------------------
447  maxlp1 = maxl + 1
448  lv = 1
449  lr = lv + n*maxlp1
450  lhes = lr + n + 1
451  lq = lhes + maxl*maxlp1
452  ldl = lq + 2*maxl
453  lw = ldl + n
454  lxl = lw + n
455  lz = lxl + n
456 C
457 C Load IGWK(6) with required minimum length of the RGWK array.
458  igwk(6) = lz + n - 1
459  IF( lz+n-1.GT.lrgw ) GOTO 640
460 C ------------------------------------------------------------------
461 C Calculate scaled-preconditioned norm of RHS vector b.
462 C ------------------------------------------------------------------
463  IF (jpre .LT. 0) THEN
464  CALL msolve(n, b, rgwk(lr), nelt, ia, ja, a, isym,
465  $ rwork, iwork)
466  nms = nms + 1
467  ELSE
468  CALL dcopy(n, b, 1, rgwk(lr), 1)
469  ENDIF
470  IF( jscal.EQ.2 .OR. jscal.EQ.3 ) THEN
471  sum = 0
472  DO 10 i = 1,n
473  sum = sum + (rgwk(lr-1+i)*sb(i))**2
474  10 CONTINUE
475  bnrm = sqrt(sum)
476  ELSE
477  bnrm = dnrm2(n,rgwk(lr),1)
478  ENDIF
479 C ------------------------------------------------------------------
480 C Calculate initial residual.
481 C ------------------------------------------------------------------
482  CALL matvec(n, x, rgwk(lr), nelt, ia, ja, a, isym)
483  DO 50 i = 1,n
484  rgwk(lr-1+i) = b(i) - rgwk(lr-1+i)
485  50 CONTINUE
486 C ------------------------------------------------------------------
487 C If performing restarting, then load the residual into the
488 C correct location in the RGWK array.
489 C ------------------------------------------------------------------
490  100 CONTINUE
491  IF( nrsts.GT.nrmax ) GOTO 610
492  IF( nrsts.GT.0 ) THEN
493 C Copy the current residual to a different location in the RGWK
494 C array.
495  CALL dcopy(n, rgwk(ldl), 1, rgwk(lr), 1)
496  ENDIF
497 C ------------------------------------------------------------------
498 C Use the DPIGMR algorithm to solve the linear system A*Z = R.
499 C ------------------------------------------------------------------
500  CALL dpigmr(n, rgwk(lr), sb, sx, jscal, maxl, maxlp1, kmp,
501  $ nrsts, jpre, matvec, msolve, nmsl, rgwk(lz), rgwk(lv),
502  $ rgwk(lhes), rgwk(lq), lgmr, rwork, iwork, rgwk(lw),
503  $ rgwk(ldl), rhol, nrmax, b, bnrm, x, rgwk(lxl), itol,
504  $ tol, nelt, ia, ja, a, isym, iunit, iflag, err)
505  iter = iter + lgmr
506  nms = nms + nmsl
507 C
508 C Increment X by the current approximate solution Z of A*Z = R.
509 C
510  lzm1 = lz - 1
511  DO 110 i = 1,n
512  x(i) = x(i) + rgwk(lzm1+i)
513  110 CONTINUE
514  IF( iflag.EQ.0 ) GOTO 600
515  IF( iflag.EQ.1 ) THEN
516  nrsts = nrsts + 1
517  GOTO 100
518  ENDIF
519  IF( iflag.EQ.2 ) GOTO 620
520 C ------------------------------------------------------------------
521 C All returns are made through this section.
522 C ------------------------------------------------------------------
523 C The iteration has converged.
524 C
525  600 CONTINUE
526  igwk(7) = nms
527  rgwk(1) = rhol
528  ierr = 0
529  RETURN
530 C
531 C Max number((NRMAX+1)*MAXL) of linear iterations performed.
532  610 CONTINUE
533  igwk(7) = nms
534  rgwk(1) = rhol
535  ierr = 1
536  RETURN
537 C
538 C GMRES failed to reduce last residual in MAXL iterations.
539 C The iteration has stalled.
540  620 CONTINUE
541  igwk(7) = nms
542  rgwk(1) = rhol
543  ierr = 2
544  RETURN
545 C Error return. Insufficient length for RGWK array.
546  640 CONTINUE
547  err = tol
548  ierr = -1
549  RETURN
550 C Error return. Inconsistent ITOL and JPRE values.
551  650 CONTINUE
552  err = tol
553  ierr = -2
554  RETURN
555 C------------- LAST LINE OF DGMRES FOLLOWS ----------------------------
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
Definition: dgmres.f:1381
matvec
subroutine matvec(n, x, y, nelt, ia, ja, a, isym)
Definition: matvec.f:27
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
dpigmr
subroutine dpigmr(N, R0, SR, SZ, JSCAL, MAXL, MAXLP1, KMP, NRSTS, JPRE, MATVEC, MSOLVE, NMSL, Z, V, HES, Q, LGMR, RPAR, IPAR, WK, DL, RHOL, NRMAX, B, BNRM, X, XL, ITOL, TOL, NELT, IA, JA, A, ISYM, IUNIT, IFLAG, ERR)
Definition: dgmres.f:724
dnrm2
double precision function dnrm2(N, DX, INCX)
Definition: dgmres.f:559
msolve
subroutine msolve(n, r, z, nelt, ia, ja, a, isym, rwork, iwork)
Definition: msolve.f:22

◆ dhels()

subroutine dhels ( double precision, dimension(lda,*)  A,
integer  LDA,
integer  N,
double precision, dimension(*)  Q,
double precision, dimension(*)  B 
)
1481 C***BEGIN PROLOGUE DHELS
1482 C***SUBSIDIARY
1483 C***PURPOSE Internal routine for DGMRES.
1484 C***LIBRARY SLATEC (SLAP)
1485 C***CATEGORY D2A4, D2B4
1486 C***TYPE DOUBLE PRECISION (SHELS-S, DHELS-D)
1487 C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
1488 C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
1489 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
1490 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
1491 C Seager, Mark K., (LLNL), seager@llnl.gov
1492 C Lawrence Livermore National Laboratory
1493 C PO Box 808, L-60
1494 C Livermore, CA 94550 (510) 423-3141
1495 C***DESCRIPTION
1496 C This routine is extracted from the LINPACK routine SGESL with
1497 C changes due to the fact that A is an upper Hessenberg matrix.
1498 C
1499 C DHELS solves the least squares problem:
1500 C
1501 C MIN(B-A*X,B-A*X)
1502 C
1503 C using the factors computed by DHEQR.
1504 C
1505 C *Usage:
1506 C INTEGER LDA, N
1507 C DOUBLE PRECISION A(LDA,N), Q(2*N), B(N+1)
1508 C
1509 C CALL DHELS(A, LDA, N, Q, B)
1510 C
1511 C *Arguments:
1512 C A :IN Double Precision A(LDA,N)
1513 C The output from DHEQR which contains the upper
1514 C triangular factor R in the QR decomposition of A.
1515 C LDA :IN Integer
1516 C The leading dimension of the array A.
1517 C N :IN Integer
1518 C A is originally an (N+1) by N matrix.
1519 C Q :IN Double Precision Q(2*N)
1520 C The coefficients of the N Givens rotations
1521 C used in the QR factorization of A.
1522 C B :INOUT Double Precision B(N+1)
1523 C On input, B is the right hand side vector.
1524 C On output, B is the solution vector X.
1525 C
1526 C***SEE ALSO DGMRES
1527 C***ROUTINES CALLED DAXPY
1528 C***REVISION HISTORY (YYMMDD)
1529 C 890404 DATE WRITTEN
1530 C 890404 Previous REVISION DATE
1531 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
1532 C 890922 Numerous changes to prologue to make closer to SLATEC
1533 C standard. (FNF)
1534 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
1535 C 910411 Prologue converted to Version 4.0 format. (BAB)
1536 C 910502 Added C***FIRST EXECUTABLE STATEMENT line. (FNF)
1537 C 910506 Made subsidiary to DGMRES. (FNF)
1538 C 920511 Added complete declaration section. (WRB)
1539 C***END PROLOGUE DHELS
1540 C The following is for optimized compilation on LLNL/LTSS Crays.
1541 CLLL. OPTIMIZE
1542 C .. Scalar Arguments ..
1543  INTEGER lda, n
1544 C .. Array Arguments ..
1545  DOUBLE PRECISION a(lda,*), b(*), q(*)
1546 C .. Local Scalars ..
1547  DOUBLE PRECISION c, s, t, t1, t2
1548  INTEGER iq, k, kb, kp1
1549 C .. External Subroutines ..
1550  EXTERNAL daxpy
1551 C***FIRST EXECUTABLE STATEMENT DHELS
1552 C
1553 C Minimize(B-A*X,B-A*X). First form Q*B.
1554 C
1555  DO 20 k = 1, n
1556  kp1 = k + 1
1557  iq = 2*(k-1) + 1
1558  c = q(iq)
1559  s = q(iq+1)
1560  t1 = b(k)
1561  t2 = b(kp1)
1562  b(k) = c*t1 - s*t2
1563  b(kp1) = s*t1 + c*t2
1564  20 CONTINUE
1565 C
1566 C Now solve R*X = Q*B.
1567 C
1568  DO 40 kb = 1, n
1569  k = n + 1 - kb
1570  b(k) = b(k)/a(k,k)
1571  t = -b(k)
1572  CALL daxpy(k-1, t, a(1,k), 1, b(1), 1)
1573  40 CONTINUE
1574  RETURN
1575 C------------- LAST LINE OF DHELS FOLLOWS ----------------------------
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276

◆ dheqr()

subroutine dheqr ( double precision, dimension(lda,*)  A,
integer  LDA,
integer  N,
double precision, dimension(*)  Q,
integer  INFO,
integer  IJOB 
)
1579 C***BEGIN PROLOGUE DHEQR
1580 C***SUBSIDIARY
1581 C***PURPOSE Internal routine for DGMRES.
1582 C***LIBRARY SLATEC (SLAP)
1583 C***CATEGORY D2A4, D2B4
1584 C***TYPE DOUBLE PRECISION (SHEQR-S, DHEQR-D)
1585 C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
1586 C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
1587 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
1588 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
1589 C Seager, Mark K., (LLNL), seager@llnl.gov
1590 C Lawrence Livermore National Laboratory
1591 C PO Box 808, L-60
1592 C Livermore, CA 94550 (510) 423-3141
1593 C***DESCRIPTION
1594 C This routine performs a QR decomposition of an upper
1595 C Hessenberg matrix A using Givens rotations. There are two
1596 C options available: 1) Performing a fresh decomposition 2)
1597 C updating the QR factors by adding a row and a column to the
1598 C matrix A.
1599 C
1600 C *Usage:
1601 C INTEGER LDA, N, INFO, IJOB
1602 C DOUBLE PRECISION A(LDA,N), Q(2*N)
1603 C
1604 C CALL DHEQR(A, LDA, N, Q, INFO, IJOB)
1605 C
1606 C *Arguments:
1607 C A :INOUT Double Precision A(LDA,N)
1608 C On input, the matrix to be decomposed.
1609 C On output, the upper triangular matrix R.
1610 C The factorization can be written Q*A = R, where
1611 C Q is a product of Givens rotations and R is upper
1612 C triangular.
1613 C LDA :IN Integer
1614 C The leading dimension of the array A.
1615 C N :IN Integer
1616 C A is an (N+1) by N Hessenberg matrix.
1617 C Q :OUT Double Precision Q(2*N)
1618 C The factors c and s of each Givens rotation used
1619 C in decomposing A.
1620 C INFO :OUT Integer
1621 C = 0 normal value.
1622 C = K if A(K,K) .eq. 0.0 . This is not an error
1623 C condition for this subroutine, but it does
1624 C indicate that DHELS will divide by zero
1625 C if called.
1626 C IJOB :IN Integer
1627 C = 1 means that a fresh decomposition of the
1628 C matrix A is desired.
1629 C .ge. 2 means that the current decomposition of A
1630 C will be updated by the addition of a row
1631 C and a column.
1632 C
1633 C***SEE ALSO DGMRES
1634 C***ROUTINES CALLED (NONE)
1635 C***REVISION HISTORY (YYMMDD)
1636 C 890404 DATE WRITTEN
1637 C 890404 Previous REVISION DATE
1638 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
1639 C 890922 Numerous changes to prologue to make closer to SLATEC
1640 C standard. (FNF)
1641 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
1642 C 910411 Prologue converted to Version 4.0 format. (BAB)
1643 C 910506 Made subsidiary to DGMRES. (FNF)
1644 C 920511 Added complete declaration section. (WRB)
1645 C***END PROLOGUE DHEQR
1646 C The following is for optimized compilation on LLNL/LTSS Crays.
1647 CLLL. OPTIMIZE
1648 C .. Scalar Arguments ..
1649  INTEGER ijob, info, lda, n
1650 C .. Array Arguments ..
1651  DOUBLE PRECISION a(lda,*), q(*)
1652 C .. Local Scalars ..
1653  DOUBLE PRECISION c, s, t, t1, t2
1654  INTEGER i, iq, j, k, km1, kp1, nm1
1655 C .. Intrinsic Functions ..
1656  INTRINSIC abs, sqrt
1657 C***FIRST EXECUTABLE STATEMENT DHEQR
1658  IF (ijob .GT. 1) GO TO 70
1659 C -------------------------------------------------------------------
1660 C A new factorization is desired.
1661 C -------------------------------------------------------------------
1662 C QR decomposition without pivoting.
1663 C
1664  info = 0
1665  DO 60 k = 1, n
1666  km1 = k - 1
1667  kp1 = k + 1
1668 C
1669 C Compute K-th column of R.
1670 C First, multiply the K-th column of A by the previous
1671 C K-1 Givens rotations.
1672 C
1673  IF (km1 .LT. 1) GO TO 20
1674  DO 10 j = 1, km1
1675  i = 2*(j-1) + 1
1676  t1 = a(j,k)
1677  t2 = a(j+1,k)
1678  c = q(i)
1679  s = q(i+1)
1680  a(j,k) = c*t1 - s*t2
1681  a(j+1,k) = s*t1 + c*t2
1682  10 CONTINUE
1683 C
1684 C Compute Givens components C and S.
1685 C
1686  20 CONTINUE
1687  iq = 2*km1 + 1
1688  t1 = a(k,k)
1689  t2 = a(kp1,k)
1690  IF( t2.EQ.0.0d0 ) THEN
1691  c = 1
1692  s = 0
1693  ELSEIF( abs(t2).GE.abs(t1) ) THEN
1694  t = t1/t2
1695  s = -1.0d0/sqrt(1.0d0+t*t)
1696  c = -s*t
1697  ELSE
1698  t = t2/t1
1699  c = 1.0d0/sqrt(1.0d0+t*t)
1700  s = -c*t
1701  ENDIF
1702  q(iq) = c
1703  q(iq+1) = s
1704  a(k,k) = c*t1 - s*t2
1705  IF( a(k,k).EQ.0.0d0 ) info = k
1706  60 CONTINUE
1707  RETURN
1708 C -------------------------------------------------------------------
1709 C The old factorization of a will be updated. A row and a
1710 C column has been added to the matrix A. N by N-1 is now
1711 C the old size of the matrix.
1712 C -------------------------------------------------------------------
1713  70 CONTINUE
1714  nm1 = n - 1
1715 C -------------------------------------------------------------------
1716 C Multiply the new column by the N previous Givens rotations.
1717 C -------------------------------------------------------------------
1718  DO 100 k = 1,nm1
1719  i = 2*(k-1) + 1
1720  t1 = a(k,n)
1721  t2 = a(k+1,n)
1722  c = q(i)
1723  s = q(i+1)
1724  a(k,n) = c*t1 - s*t2
1725  a(k+1,n) = s*t1 + c*t2
1726  100 CONTINUE
1727 C -------------------------------------------------------------------
1728 C Complete update of decomposition by forming last Givens
1729 C rotation, and multiplying it times the column
1730 C vector(A(N,N),A(NP1,N)).
1731 C -------------------------------------------------------------------
1732  info = 0
1733  t1 = a(n,n)
1734  t2 = a(n+1,n)
1735  IF ( t2.EQ.0.0d0 ) THEN
1736  c = 1
1737  s = 0
1738  ELSEIF( abs(t2).GE.abs(t1) ) THEN
1739  t = t1/t2
1740  s = -1.0d0/sqrt(1.0d0+t*t)
1741  c = -s*t
1742  ELSE
1743  t = t2/t1
1744  c = 1.0d0/sqrt(1.0d0+t*t)
1745  s = -c*t
1746  ENDIF
1747  iq = 2*n - 1
1748  q(iq) = c
1749  q(iq+1) = s
1750  a(n,n) = c*t1 - s*t2
1751  IF (a(n,n) .EQ. 0.0d0) info = n
1752  RETURN
1753 C------------- LAST LINE OF DHEQR FOLLOWS ----------------------------

◆ dnrm2()

double precision function dnrm2 (   N,
double precision, dimension(*)  DX,
  INCX 
)
559 C***BEGIN PROLOGUE DNRM2
560 C***PURPOSE Compute the Euclidean length (L2 norm) of a vector.
561 C***LIBRARY SLATEC (BLAS)
562 C***CATEGORY D1A3B
563 C***TYPE DOUBLE PRECISION (SNRM2-S, DNRM2-D, SCNRM2-C)
564 C***KEYWORDS BLAS, EUCLIDEAN LENGTH, EUCLIDEAN NORM, L2,
565 C LINEAR ALGEBRA, UNITARY, VECTOR
566 C***AUTHOR Lawson, C. L., (JPL)
567 C Hanson, R. J., (SNLA)
568 C Kincaid, D. R., (U. of Texas)
569 C Krogh, F. T., (JPL)
570 C***DESCRIPTION
571 C
572 C B L A S Subprogram
573 C Description of parameters
574 C
575 C --Input--
576 C N number of elements in input vector(s)
577 C DX double precision vector with N elements
578 C INCX storage spacing between elements of DX
579 C
580 C --Output--
581 C DNRM2 double precision result (zero if N .LE. 0)
582 C
583 C Euclidean norm of the N-vector stored in DX with storage
584 C increment INCX.
585 C If N .LE. 0, return with result = 0.
586 C If N .GE. 1, then INCX must be .GE. 1
587 C
588 C Four phase method using two built-in constants that are
589 C hopefully applicable to all machines.
590 C CUTLO = maximum of SQRT(U/EPS) over all known machines.
591 C CUTHI = minimum of SQRT(V) over all known machines.
592 C where
593 C EPS = smallest no. such that EPS + 1. .GT. 1.
594 C U = smallest positive no. (underflow limit)
595 C V = largest no. (overflow limit)
596 C
597 C Brief outline of algorithm.
598 C
599 C Phase 1 scans zero components.
600 C move to phase 2 when a component is nonzero and .LE. CUTLO
601 C move to phase 3 when a component is .GT. CUTLO
602 C move to phase 4 when a component is .GE. CUTHI/M
603 C where M = N for X() real and M = 2*N for complex.
604 C
605 C Values for CUTLO and CUTHI.
606 C From the environmental parameters listed in the IMSL converter
607 C document the limiting values are as follows:
608 C CUTLO, S.P. U/EPS = 2**(-102) for Honeywell. Close seconds are
609 C Univac and DEC at 2**(-103)
610 C Thus CUTLO = 2**(-51) = 4.44089E-16
611 C CUTHI, S.P. V = 2**127 for Univac, Honeywell, and DEC.
612 C Thus CUTHI = 2**(63.5) = 1.30438E19
613 C CUTLO, D.P. U/EPS = 2**(-67) for Honeywell and DEC.
614 C Thus CUTLO = 2**(-33.5) = 8.23181D-11
615 C CUTHI, D.P. same as S.P. CUTHI = 1.30438D19
616 C DATA CUTLO, CUTHI /8.232D-11, 1.304D19/
617 C DATA CUTLO, CUTHI /4.441E-16, 1.304E19/
618 C
619 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
620 C Krogh, Basic linear algebra subprograms for Fortran
621 C usage, Algorithm No. 539, Transactions on Mathematical
622 C Software 5, 3 (September 1979), pp. 308-323.
623 C***ROUTINES CALLED (NONE)
624 C***REVISION HISTORY (YYMMDD)
625 C 791001 DATE WRITTEN
626 C 890531 Changed all specific intrinsics to generic. (WRB)
627 C 890831 Modified array declarations. (WRB)
628 C 890831 REVISION DATE from Version 3.2
629 C 891214 Prologue converted to Version 4.0 format. (BAB)
630 C 920501 Reformatted the REFERENCES section. (WRB)
631 C***END PROLOGUE DNRM2
632  INTEGER next
633  DOUBLE PRECISION dx(*), cutlo, cuthi, hitest, sum, xmax, zero,
634  + one
635  SAVE cutlo, cuthi, zero, one
636  DATA zero, one /0.0d0, 1.0d0/
637 C
638  DATA cutlo, cuthi /8.232d-11, 1.304d19/
639 C***FIRST EXECUTABLE STATEMENT DNRM2
640  IF (n .GT. 0) GO TO 10
641  dnrm2 = zero
642  GO TO 300
643 C
644  10 assign 30 to next
645  sum = zero
646  nn = n * incx
647 C
648 C BEGIN MAIN LOOP
649 C
650  i = 1
651  20 GO TO next,(30, 50, 70, 110)
652  30 IF (abs(dx(i)) .GT. cutlo) GO TO 85
653  assign 50 to next
654  xmax = zero
655 C
656 C PHASE 1. SUM IS ZERO
657 C
658  50 IF (dx(i) .EQ. zero) GO TO 200
659  IF (abs(dx(i)) .GT. cutlo) GO TO 85
660 C
661 C PREPARE FOR PHASE 2.
662 C
663  assign 70 to next
664  GO TO 105
665 C
666 C PREPARE FOR PHASE 4.
667 C
668  100 i = j
669  assign 110 to next
670  sum = (sum / dx(i)) / dx(i)
671  105 xmax = abs(dx(i))
672  GO TO 115
673 C
674 C PHASE 2. SUM IS SMALL.
675 C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
676 C
677  70 IF (abs(dx(i)) .GT. cutlo) GO TO 75
678 C
679 C COMMON CODE FOR PHASES 2 AND 4.
680 C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
681 C
682  110 IF (abs(dx(i)) .LE. xmax) GO TO 115
683  sum = one + sum * (xmax / dx(i))**2
684  xmax = abs(dx(i))
685  GO TO 200
686 C
687  115 sum = sum + (dx(i)/xmax)**2
688  GO TO 200
689 C
690 C PREPARE FOR PHASE 3.
691 C
692  75 sum = (sum * xmax) * xmax
693 C
694 C FOR REAL OR D.P. SET HITEST = CUTHI/N
695 C FOR COMPLEX SET HITEST = CUTHI/(2*N)
696 C
697  85 hitest = cuthi / n
698 C
699 C PHASE 3. SUM IS MID-RANGE. NO SCALING.
700 C
701  DO 95 j = i,nn,incx
702  IF (abs(dx(j)) .GE. hitest) GO TO 100
703  95 sum = sum + dx(j)**2
704  dnrm2 = sqrt(sum)
705  GO TO 300
706 C
707  200 CONTINUE
708  i = i + incx
709  IF (i .LE. nn) GO TO 20
710 C
711 C END OF MAIN LOOP.
712 C
713 C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
714 C
715  dnrm2 = xmax * sqrt(sum)
716  300 CONTINUE
717  RETURN
dnrm2
double precision function dnrm2(N, DX, INCX)
Definition: dgmres.f:559

◆ dorth()

subroutine dorth ( double precision, dimension(*)  VNEW,
double precision, dimension(n,*)  V,
double precision, dimension(ldhes,*)  HES,
integer  N,
integer  LL,
integer  LDHES,
integer  KMP,
double precision  SNORMW 
)
2159 C***BEGIN PROLOGUE DORTH
2160 C***SUBSIDIARY
2161 C***PURPOSE Internal routine for DGMRES.
2162 C***LIBRARY SLATEC (SLAP)
2163 C***CATEGORY D2A4, D2B4
2164 C***TYPE DOUBLE PRECISION (SORTH-S, DORTH-D)
2165 C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
2166 C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
2167 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
2168 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
2169 C Seager, Mark K., (LLNL), seager@llnl.gov
2170 C Lawrence Livermore National Laboratory
2171 C PO Box 808, L-60
2172 C Livermore, CA 94550 (510) 423-3141
2173 C***DESCRIPTION
2174 C This routine orthogonalizes the vector VNEW against the
2175 C previous KMP vectors in the V array. It uses a modified
2176 C Gram-Schmidt orthogonalization procedure with conditional
2177 C reorthogonalization.
2178 C
2179 C *Usage:
2180 C INTEGER N, LL, LDHES, KMP
2181 C DOUBLE PRECISION VNEW(N), V(N,LL), HES(LDHES,LL), SNORMW
2182 C
2183 C CALL DORTH(VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
2184 C
2185 C *Arguments:
2186 C VNEW :INOUT Double Precision VNEW(N)
2187 C On input, the vector of length N containing a scaled
2188 C product of the Jacobian and the vector V(*,LL).
2189 C On output, the new vector orthogonal to V(*,i0) to V(*,LL),
2190 C where i0 = max(1, LL-KMP+1).
2191 C V :IN Double Precision V(N,LL)
2192 C The N x LL array containing the previous LL
2193 C orthogonal vectors V(*,1) to V(*,LL).
2194 C HES :INOUT Double Precision HES(LDHES,LL)
2195 C On input, an LL x LL upper Hessenberg matrix containing,
2196 C in HES(I,K), K.lt.LL, the scaled inner products of
2197 C A*V(*,K) and V(*,i).
2198 C On return, column LL of HES is filled in with
2199 C the scaled inner products of A*V(*,LL) and V(*,i).
2200 C N :IN Integer
2201 C The order of the matrix A, and the length of VNEW.
2202 C LL :IN Integer
2203 C The current order of the matrix HES.
2204 C LDHES :IN Integer
2205 C The leading dimension of the HES array.
2206 C KMP :IN Integer
2207 C The number of previous vectors the new vector VNEW
2208 C must be made orthogonal to (KMP .le. MAXL).
2209 C SNORMW :OUT DOUBLE PRECISION
2210 C Scalar containing the l-2 norm of VNEW.
2211 C
2212 C***SEE ALSO DGMRES
2213 C***ROUTINES CALLED DAXPY, DDOT, DNRM2
2214 C***REVISION HISTORY (YYMMDD)
2215 C 890404 DATE WRITTEN
2216 C 890404 Previous REVISION DATE
2217 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
2218 C 890922 Numerous changes to prologue to make closer to SLATEC
2219 C standard. (FNF)
2220 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
2221 C 910411 Prologue converted to Version 4.0 format. (BAB)
2222 C 910506 Made subsidiary to DGMRES. (FNF)
2223 C 920511 Added complete declaration section. (WRB)
2224 C***END PROLOGUE DORTH
2225 C The following is for optimized compilation on LLNL/LTSS Crays.
2226 CLLL. OPTIMIZE
2227 C .. Scalar Arguments ..
2228  DOUBLE PRECISION snormw
2229  INTEGER kmp, ldhes, ll, n
2230 C .. Array Arguments ..
2231  DOUBLE PRECISION hes(ldhes,*), v(n,*), vnew(*)
2232 C .. Local Scalars ..
2233  DOUBLE PRECISION arg, sumdsq, tem, vnrm
2234  INTEGER i, i0
2235 C .. External Functions ..
2236  DOUBLE PRECISION ddot, dnrm2
2237  EXTERNAL ddot, dnrm2
2238 C .. External Subroutines ..
2239  EXTERNAL daxpy
2240 C .. Intrinsic Functions ..
2241  INTRINSIC max, sqrt
2242 C***FIRST EXECUTABLE STATEMENT DORTH
2243 C
2244 C Get norm of unaltered VNEW for later use.
2245 C
2246  vnrm = dnrm2(n, vnew, 1)
2247 C -------------------------------------------------------------------
2248 C Perform the modified Gram-Schmidt procedure on VNEW =A*V(LL).
2249 C Scaled inner products give new column of HES.
2250 C Projections of earlier vectors are subtracted from VNEW.
2251 C -------------------------------------------------------------------
2252  i0 = max(1,ll-kmp+1)
2253  DO 10 i = i0,ll
2254  hes(i,ll) = ddot(n, v(1,i), 1, vnew, 1)
2255  tem = -hes(i,ll)
2256  CALL daxpy(n, tem, v(1,i), 1, vnew, 1)
2257  10 CONTINUE
2258 C -------------------------------------------------------------------
2259 C Compute SNORMW = norm of VNEW. If VNEW is small compared
2260 C to its input value (in norm), then reorthogonalize VNEW to
2261 C V(*,1) through V(*,LL). Correct if relative correction
2262 C exceeds 1000*(unit roundoff). Finally, correct SNORMW using
2263 C the dot products involved.
2264 C -------------------------------------------------------------------
2265  snormw = dnrm2(n, vnew, 1)
2266  IF (vnrm + 0.001d0*snormw .NE. vnrm) RETURN
2267  sumdsq = 0
2268  DO 30 i = i0,ll
2269  tem = -ddot(n, v(1,i), 1, vnew, 1)
2270  IF (hes(i,ll) + 0.001d0*tem .EQ. hes(i,ll)) GO TO 30
2271  hes(i,ll) = hes(i,ll) - tem
2272  CALL daxpy(n, tem, v(1,i), 1, vnew, 1)
2273  sumdsq = sumdsq + tem**2
2274  30 CONTINUE
2275  IF (sumdsq .EQ. 0.0d0) RETURN
2276  arg = max(0.0d0,snormw**2 - sumdsq)
2277  snormw = sqrt(arg)
2278 C
2279  RETURN
2280 C------------- LAST LINE OF DORTH FOLLOWS ----------------------------
max
#define max(a, b)
Definition: cascade.c:32
dnrm2
double precision function dnrm2(N, DX, INCX)
Definition: dgmres.f:559
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276
ddot
double precision function ddot(N, DX, INCX, DY, INCY)
Definition: dgmres.f:2469

◆ dpigmr()

subroutine dpigmr ( integer  N,
double precision, dimension(*)  R0,
double precision, dimension(*)  SR,
double precision, dimension(*)  SZ,
integer  JSCAL,
integer  MAXL,
integer  MAXLP1,
integer  KMP,
integer  NRSTS,
integer  JPRE,
external  MATVEC,
external  MSOLVE,
integer  NMSL,
double precision, dimension(*)  Z,
double precision, dimension(n,*)  V,
double precision, dimension(maxlp1,*)  HES,
double precision, dimension(*)  Q,
integer  LGMR,
double precision, dimension(*)  RPAR,
integer, dimension(*)  IPAR,
double precision, dimension(*)  WK,
double precision, dimension(*)  DL,
double precision  RHOL,
integer  NRMAX,
double precision, dimension(*)  B,
double precision  BNRM,
double precision, dimension(*)  X,
double precision, dimension(*)  XL,
integer  ITOL,
double precision  TOL,
integer  NELT,
integer, dimension(nelt)  IA,
integer, dimension(nelt)  JA,
double precision, dimension(nelt)  A,
integer  ISYM,
integer  IUNIT,
integer  IFLAG,
double precision  ERR 
)
724 C***BEGIN PROLOGUE DPIGMR
725 C***SUBSIDIARY
726 C***PURPOSE Internal routine for DGMRES.
727 C***LIBRARY SLATEC (SLAP)
728 C***CATEGORY D2A4, D2B4
729 C***TYPE DOUBLE PRECISION (SPIGMR-S, DPIGMR-D)
730 C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
731 C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
732 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
733 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
734 C Seager, Mark K., (LLNL), seager@llnl.gov
735 C Lawrence Livermore National Laboratory
736 C PO Box 808, L-60
737 C Livermore, CA 94550 (510) 423-3141
738 C***DESCRIPTION
739 C This routine solves the linear system A * Z = R0 using a
740 C scaled preconditioned version of the generalized minimum
741 C residual method. An initial guess of Z = 0 is assumed.
742 C
743 C *Usage:
744 C INTEGER N, JSCAL, MAXL, MAXLP1, KMP, NRSTS, JPRE, NMSL, LGMR
745 C INTEGER IPAR(USER DEFINED), NRMAX, ITOL, NELT, IA(NELT), JA(NELT)
746 C INTEGER ISYM, IUNIT, IFLAG
747 C DOUBLE PRECISION R0(N), SR(N), SZ(N), Z(N), V(N,MAXLP1),
748 C $ HES(MAXLP1,MAXL), Q(2*MAXL), RPAR(USER DEFINED),
749 C $ WK(N), DL(N), RHOL, B(N), BNRM, X(N), XL(N),
750 C $ TOL, A(NELT), ERR
751 C EXTERNAL MATVEC, MSOLVE
752 C
753 C CALL DPIGMR(N, R0, SR, SZ, JSCAL, MAXL, MAXLP1, KMP,
754 C $ NRSTS, JPRE, MATVEC, MSOLVE, NMSL, Z, V, HES, Q, LGMR,
755 C $ RPAR, IPAR, WK, DL, RHOL, NRMAX, B, BNRM, X, XL,
756 C $ ITOL, TOL, NELT, IA, JA, A, ISYM, IUNIT, IFLAG, ERR)
757 C
758 C *Arguments:
759 C N :IN Integer
760 C The order of the matrix A, and the lengths
761 C of the vectors SR, SZ, R0 and Z.
762 C R0 :IN Double Precision R0(N)
763 C R0 = the right hand side of the system A*Z = R0.
764 C R0 is also used as workspace when computing
765 C the final approximation.
766 C (R0 is the same as V(*,MAXL+1) in the call to DPIGMR.)
767 C SR :IN Double Precision SR(N)
768 C SR is a vector of length N containing the non-zero
769 C elements of the diagonal scaling matrix for R0.
770 C SZ :IN Double Precision SZ(N)
771 C SZ is a vector of length N containing the non-zero
772 C elements of the diagonal scaling matrix for Z.
773 C JSCAL :IN Integer
774 C A flag indicating whether arrays SR and SZ are used.
775 C JSCAL=0 means SR and SZ are not used and the
776 C algorithm will perform as if all
777 C SR(i) = 1 and SZ(i) = 1.
778 C JSCAL=1 means only SZ is used, and the algorithm
779 C performs as if all SR(i) = 1.
780 C JSCAL=2 means only SR is used, and the algorithm
781 C performs as if all SZ(i) = 1.
782 C JSCAL=3 means both SR and SZ are used.
783 C MAXL :IN Integer
784 C The maximum allowable order of the matrix H.
785 C MAXLP1 :IN Integer
786 C MAXPL1 = MAXL + 1, used for dynamic dimensioning of HES.
787 C KMP :IN Integer
788 C The number of previous vectors the new vector VNEW
789 C must be made orthogonal to. (KMP .le. MAXL)
790 C NRSTS :IN Integer
791 C Counter for the number of restarts on the current
792 C call to DGMRES. If NRSTS .gt. 0, then the residual
793 C R0 is already scaled, and so scaling of it is
794 C not necessary.
795 C JPRE :IN Integer
796 C Preconditioner type flag.
797 C MATVEC :EXT External.
798 C Name of a routine which performs the matrix vector multiply
799 C Y = A*X given A and X. The name of the MATVEC routine must
800 C be declared external in the calling program. The calling
801 C sequence to MATVEC is:
802 C CALL MATVEC(N, X, Y, NELT, IA, JA, A, ISYM)
803 C where N is the number of unknowns, Y is the product A*X
804 C upon return, X is an input vector, and NELT is the number of
805 C non-zeros in the SLAP IA, JA, A storage for the matrix A.
806 C ISYM is a flag which, if non-zero, denotes that A is
807 C symmetric and only the lower or upper triangle is stored.
808 C MSOLVE :EXT External.
809 C Name of the routine which solves a linear system Mz = r for
810 C z given r with the preconditioning matrix M (M is supplied via
811 C RPAR and IPAR arrays. The name of the MSOLVE routine must
812 C be declared external in the calling program. The calling
813 C sequence to MSOLVE is:
814 C CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RPAR, IPAR)
815 C Where N is the number of unknowns, R is the right-hand side
816 C vector and Z is the solution upon return. NELT, IA, JA, A and
817 C ISYM are defined as below. RPAR is a double precision array
818 C that can be used to pass necessary preconditioning information
819 C and/or workspace to MSOLVE. IPAR is an integer work array
820 C for the same purpose as RPAR.
821 C NMSL :OUT Integer
822 C The number of calls to MSOLVE.
823 C Z :OUT Double Precision Z(N)
824 C The final computed approximation to the solution
825 C of the system A*Z = R0.
826 C V :OUT Double Precision V(N,MAXLP1)
827 C The N by (LGMR+1) array containing the LGMR
828 C orthogonal vectors V(*,1) to V(*,LGMR).
829 C HES :OUT Double Precision HES(MAXLP1,MAXL)
830 C The upper triangular factor of the QR decomposition
831 C of the (LGMR+1) by LGMR upper Hessenberg matrix whose
832 C entries are the scaled inner-products of A*V(*,I)
833 C and V(*,K).
834 C Q :OUT Double Precision Q(2*MAXL)
835 C A double precision array of length 2*MAXL containing the
836 C components of the Givens rotations used in the QR
837 C decomposition of HES. It is loaded in DHEQR and used in
838 C DHELS.
839 C LGMR :OUT Integer
840 C The number of iterations performed and
841 C the current order of the upper Hessenberg
842 C matrix HES.
843 C RPAR :IN Double Precision RPAR(USER DEFINED)
844 C Double Precision workspace passed directly to the MSOLVE
845 C routine.
846 C IPAR :IN Integer IPAR(USER DEFINED)
847 C Integer workspace passed directly to the MSOLVE routine.
848 C WK :IN Double Precision WK(N)
849 C A double precision work array of length N used by routines
850 C MATVEC and MSOLVE.
851 C DL :INOUT Double Precision DL(N)
852 C On input, a double precision work array of length N used for
853 C calculation of the residual norm RHO when the method is
854 C incomplete (KMP.lt.MAXL), and/or when using restarting.
855 C On output, the scaled residual vector RL. It is only loaded
856 C when performing restarts of the Krylov iteration.
857 C RHOL :OUT Double Precision
858 C A double precision scalar containing the norm of the final
859 C residual.
860 C NRMAX :IN Integer
861 C The maximum number of restarts of the Krylov iteration.
862 C NRMAX .gt. 0 means restarting is active, while
863 C NRMAX = 0 means restarting is not being used.
864 C B :IN Double Precision B(N)
865 C The right hand side of the linear system A*X = b.
866 C BNRM :IN Double Precision
867 C The scaled norm of b.
868 C X :IN Double Precision X(N)
869 C The current approximate solution as of the last
870 C restart.
871 C XL :IN Double Precision XL(N)
872 C An array of length N used to hold the approximate
873 C solution X(L) when ITOL=11.
874 C ITOL :IN Integer
875 C A flag to indicate the type of convergence criterion
876 C used. See the driver for its description.
877 C TOL :IN Double Precision
878 C The tolerance on residuals R0-A*Z in scaled norm.
879 C NELT :IN Integer
880 C The length of arrays IA, JA and A.
881 C IA :IN Integer IA(NELT)
882 C An integer array of length NELT containing matrix data.
883 C It is passed directly to the MATVEC and MSOLVE routines.
884 C JA :IN Integer JA(NELT)
885 C An integer array of length NELT containing matrix data.
886 C It is passed directly to the MATVEC and MSOLVE routines.
887 C A :IN Double Precision A(NELT)
888 C A double precision array of length NELT containing matrix
889 C data. It is passed directly to the MATVEC and MSOLVE routines.
890 C ISYM :IN Integer
891 C A flag to indicate symmetric matrix storage.
892 C If ISYM=0, all non-zero entries of the matrix are
893 C stored. If ISYM=1, the matrix is symmetric and
894 C only the upper or lower triangular part is stored.
895 C IUNIT :IN Integer
896 C The i/o unit number for writing intermediate residual
897 C norm values.
898 C IFLAG :OUT Integer
899 C An integer error flag..
900 C 0 means convergence in LGMR iterations, LGMR.le.MAXL.
901 C 1 means the convergence test did not pass in MAXL
902 C iterations, but the residual norm is .lt. norm(R0),
903 C and so Z is computed.
904 C 2 means the convergence test did not pass in MAXL
905 C iterations, residual .ge. norm(R0), and Z = 0.
906 C ERR :OUT Double Precision.
907 C Error estimate of error in final approximate solution, as
908 C defined by ITOL.
909 C
910 C *Cautions:
911 C This routine will attempt to write to the Fortran logical output
912 C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that
913 C this logical unit is attached to a file or terminal before calling
914 C this routine with a non-zero value for IUNIT. This routine does
915 C not check for the validity of a non-zero IUNIT unit number.
916 C
917 C***SEE ALSO DGMRES
918 C***ROUTINES CALLED DAXPY, DCOPY, DHELS, DHEQR, DNRM2, DORTH, DRLCAL,
919 C DSCAL, ISDGMR
920 C***REVISION HISTORY (YYMMDD)
921 C 890404 DATE WRITTEN
922 C 890404 Previous REVISION DATE
923 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
924 C 890922 Numerous changes to prologue to make closer to SLATEC
925 C standard. (FNF)
926 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
927 C 910411 Prologue converted to Version 4.0 format. (BAB)
928 C 910502 Removed MATVEC and MSOLVE from ROUTINES CALLED list. (FNF)
929 C 910506 Made subsidiary to DGMRES. (FNF)
930 C 920511 Added complete declaration section. (WRB)
931 C***END PROLOGUE DPIGMR
932 C The following is for optimized compilation on LLNL/LTSS Crays.
933 CLLL. OPTIMIZE
934 C .. Scalar Arguments ..
935  DOUBLE PRECISION bnrm, err, rhol, tol
936  INTEGER iflag, isym, itol, iunit, jpre, jscal, kmp, lgmr, maxl,
937  + maxlp1, n, nelt, nmsl, nrmax, nrsts
938 C .. Array Arguments ..
939  DOUBLE PRECISION a(nelt), b(*), dl(*), hes(maxlp1,*), q(*), r0(*),
940  + rpar(*), sr(*), sz(*), v(n,*), wk(*), x(*),
941  + xl(*), z(*)
942  INTEGER ia(nelt), ipar(*), ja(nelt)
943 C .. Subroutine Arguments ..
944  EXTERNAL matvec, msolve
945 C .. Local Scalars ..
946  DOUBLE PRECISION c, dlnrm, prod, r0nrm, rho, s, snormw, tem
947  INTEGER i, i2, info, ip1, iter, itmax, j, k, ll, llp1
948 C .. External Functions ..
949  DOUBLE PRECISION dnrm2
950  INTEGER isdgmr
951  EXTERNAL dnrm2, isdgmr
952 C .. External Subroutines ..
953  EXTERNAL daxpy, dcopy, dhels, dheqr, dorth, drlcal, dscal
954 C .. Intrinsic Functions ..
955  INTRINSIC abs
956 C***FIRST EXECUTABLE STATEMENT DPIGMR
957 C
958 C Zero out the Z array.
959 C
960  DO 5 i = 1,n
961  z(i) = 0
962  5 CONTINUE
963 C
964  iflag = 0
965  lgmr = 0
966  nmsl = 0
967 C Load ITMAX, the maximum number of iterations.
968  itmax =(nrmax+1)*maxl
969 C -------------------------------------------------------------------
970 C The initial residual is the vector R0.
971 C Apply left precon. if JPRE < 0 and this is not a restart.
972 C Apply scaling to R0 if JSCAL = 2 or 3.
973 C -------------------------------------------------------------------
974  IF ((jpre .LT. 0) .AND.(nrsts .EQ. 0)) THEN
975  CALL dcopy(n, r0, 1, wk, 1)
976  CALL msolve(n, wk, r0, nelt, ia, ja, a, isym, rpar, ipar)
977  nmsl = nmsl + 1
978  ENDIF
979  IF (((jscal.EQ.2) .OR.(jscal.EQ.3)) .AND.(nrsts.EQ.0)) THEN
980  DO 10 i = 1,n
981  v(i,1) = r0(i)*sr(i)
982  10 CONTINUE
983  ELSE
984  DO 20 i = 1,n
985  v(i,1) = r0(i)
986  20 CONTINUE
987  ENDIF
988  r0nrm = dnrm2(n, v, 1)
989  iter = nrsts*maxl
990 C
991 C Call stopping routine ISDGMR.
992 C
993  IF (isdgmr(n, b, x, xl, nelt, ia, ja, a, isym, msolve,
994  $ nmsl, itol, tol, itmax, iter, err, iunit, v(1,1), z, wk,
995  $ rpar, ipar, r0nrm, bnrm, sr, sz, jscal,
996  $ kmp, lgmr, maxl, maxlp1, v, q, snormw, prod, r0nrm,
997  $ hes, jpre) .NE. 0) RETURN
998  tem = 1.0d0/r0nrm
999  CALL dscal(n, tem, v(1,1), 1)
1000 C
1001 C Zero out the HES array.
1002 C
1003  DO 50 j = 1,maxl
1004  DO 40 i = 1,maxlp1
1005  hes(i,j) = 0
1006  40 CONTINUE
1007  50 CONTINUE
1008 C -------------------------------------------------------------------
1009 C Main loop to compute the vectors V(*,2) to V(*,MAXL).
1010 C The running product PROD is needed for the convergence test.
1011 C -------------------------------------------------------------------
1012  prod = 1
1013  DO 90 ll = 1,maxl
1014  lgmr = ll
1015 C -------------------------------------------------------------------
1016 C Unscale the current V(LL) and store in WK. Call routine
1017 C MSOLVE to compute(M-inverse)*WK, where M is the
1018 C preconditioner matrix. Save the answer in Z. Call routine
1019 C MATVEC to compute VNEW = A*Z, where A is the the system
1020 C matrix. save the answer in V(LL+1). Scale V(LL+1). Call
1021 C routine DORTH to orthogonalize the new vector VNEW =
1022 C V(*,LL+1). Call routine DHEQR to update the factors of HES.
1023 C -------------------------------------------------------------------
1024  IF ((jscal .EQ. 1) .OR.(jscal .EQ. 3)) THEN
1025  DO 60 i = 1,n
1026  wk(i) = v(i,ll)/sz(i)
1027  60 CONTINUE
1028  ELSE
1029  CALL dcopy(n, v(1,ll), 1, wk, 1)
1030  ENDIF
1031  IF (jpre .GT. 0) THEN
1032  CALL msolve(n, wk, z, nelt, ia, ja, a, isym, rpar, ipar)
1033  nmsl = nmsl + 1
1034  CALL matvec(n, z, v(1,ll+1), nelt, ia, ja, a, isym)
1035  ELSE
1036  CALL matvec(n, wk, v(1,ll+1), nelt, ia, ja, a, isym)
1037  ENDIF
1038  IF (jpre .LT. 0) THEN
1039  CALL dcopy(n, v(1,ll+1), 1, wk, 1)
1040  CALL msolve(n,wk,v(1,ll+1),nelt,ia,ja,a,isym,rpar,ipar)
1041  nmsl = nmsl + 1
1042  ENDIF
1043  IF ((jscal .EQ. 2) .OR.(jscal .EQ. 3)) THEN
1044  DO 65 i = 1,n
1045  v(i,ll+1) = v(i,ll+1)*sr(i)
1046  65 CONTINUE
1047  ENDIF
1048  CALL dorth(v(1,ll+1), v, hes, n, ll, maxlp1, kmp, snormw)
1049  hes(ll+1,ll) = snormw
1050  CALL dheqr(hes, maxlp1, ll, q, info, ll)
1051  IF (info .EQ. ll) GO TO 120
1052 C -------------------------------------------------------------------
1053 C Update RHO, the estimate of the norm of the residual R0-A*ZL.
1054 C If KMP < MAXL, then the vectors V(*,1),...,V(*,LL+1) are not
1055 C necessarily orthogonal for LL > KMP. The vector DL must then
1056 C be computed, and its norm used in the calculation of RHO.
1057 C -------------------------------------------------------------------
1058  prod = prod*q(2*ll)
1059  rho = abs(prod*r0nrm)
1060  IF ((ll.GT.kmp) .AND.(kmp.LT.maxl)) THEN
1061  IF (ll .EQ. kmp+1) THEN
1062  CALL dcopy(n, v(1,1), 1, dl, 1)
1063  DO 75 i = 1,kmp
1064  ip1 = i + 1
1065  i2 = i*2
1066  s = q(i2)
1067  c = q(i2-1)
1068  DO 70 k = 1,n
1069  dl(k) = s*dl(k) + c*v(k,ip1)
1070  70 CONTINUE
1071  75 CONTINUE
1072  ENDIF
1073  s = q(2*ll)
1074  c = q(2*ll-1)/snormw
1075  llp1 = ll + 1
1076  DO 80 k = 1,n
1077  dl(k) = s*dl(k) + c*v(k,llp1)
1078  80 CONTINUE
1079  dlnrm = dnrm2(n, dl, 1)
1080  rho = rho*dlnrm
1081  ENDIF
1082  rhol = rho
1083 C -------------------------------------------------------------------
1084 C Test for convergence. If passed, compute approximation ZL.
1085 C If failed and LL < MAXL, then continue iterating.
1086 C -------------------------------------------------------------------
1087  iter = nrsts*maxl + lgmr
1088  IF (isdgmr(n, b, x, xl, nelt, ia, ja, a, isym, msolve,
1089  $ nmsl, itol, tol, itmax, iter, err, iunit, dl, z, wk,
1090  $ rpar, ipar, rhol, bnrm, sr, sz, jscal,
1091  $ kmp, lgmr, maxl, maxlp1, v, q, snormw, prod, r0nrm,
1092  $ hes, jpre) .NE. 0) GO TO 200
1093  IF (ll .EQ. maxl) GO TO 100
1094 C -------------------------------------------------------------------
1095 C Rescale so that the norm of V(1,LL+1) is one.
1096 C -------------------------------------------------------------------
1097  tem = 1.0d0/snormw
1098  CALL dscal(n, tem, v(1,ll+1), 1)
1099  90 CONTINUE
1100  100 CONTINUE
1101  IF (rho .LT. r0nrm) GO TO 150
1102  120 CONTINUE
1103  iflag = 2
1104 C
1105 C Load approximate solution with zero.
1106 C
1107  DO 130 i = 1,n
1108  z(i) = 0
1109  130 CONTINUE
1110  RETURN
1111  150 iflag = 1
1112 C
1113 C Tolerance not met, but residual norm reduced.
1114 C
1115  IF (nrmax .GT. 0) THEN
1116 C
1117 C If performing restarting (NRMAX > 0) calculate the residual
1118 C vector RL and store it in the DL array. If the incomplete
1119 C version is being used (KMP < MAXL) then DL has already been
1120 C calculated up to a scaling factor. Use DRLCAL to calculate
1121 C the scaled residual vector.
1122 C
1123  CALL drlcal(n, kmp, maxl, maxl, v, q, dl, snormw, prod,
1124  $ r0nrm)
1125  ENDIF
1126 C -------------------------------------------------------------------
1127 C Compute the approximation ZL to the solution. Since the
1128 C vector Z was used as workspace, and the initial guess
1129 C of the linear iteration is zero, Z must be reset to zero.
1130 C -------------------------------------------------------------------
1131  200 CONTINUE
1132  ll = lgmr
1133  llp1 = ll + 1
1134  DO 210 k = 1,llp1
1135  r0(k) = 0
1136  210 CONTINUE
1137  r0(1) = r0nrm
1138  CALL dhels(hes, maxlp1, ll, q, r0)
1139  DO 220 k = 1,n
1140  z(k) = 0
1141  220 CONTINUE
1142  DO 230 i = 1,ll
1143  CALL daxpy(n, r0(i), v(1,i), 1, z, 1)
1144  230 CONTINUE
1145  IF ((jscal .EQ. 1) .OR.(jscal .EQ. 3)) THEN
1146  DO 240 i = 1,n
1147  z(i) = z(i)/sz(i)
1148  240 CONTINUE
1149  ENDIF
1150  IF (jpre .GT. 0) THEN
1151  CALL dcopy(n, z, 1, wk, 1)
1152  CALL msolve(n, wk, z, nelt, ia, ja, a, isym, rpar, ipar)
1153  nmsl = nmsl + 1
1154  ENDIF
1155  RETURN
1156 C------------- LAST LINE OF DPIGMR FOLLOWS ----------------------------
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
Definition: dgmres.f:1381
matvec
subroutine matvec(n, x, y, nelt, ia, ja, a, isym)
Definition: matvec.f:27
dorth
subroutine dorth(VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
Definition: dgmres.f:2159
dnrm2
double precision function dnrm2(N, DX, INCX)
Definition: dgmres.f:559
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276
dscal
subroutine dscal(n, da, dx, incx)
Definition: dgesv.f:1746
isdgmr
integer function isdgmr(N, B, X, XL, NELT, IA, JA, A, ISYM, MSOLVE, NMSL, ITOL, TOL, ITMAX, ITER, ERR, IUNIT, R, Z, DZ, RWORK, IWORK, RNRM, BNRM, SB, SX, JSCAL, KMP, LGMR, MAXL, MAXLP1, V, Q, SNORMW, PROD, R0NRM, HES, JPRE)
Definition: dgmres.f:1760
dheqr
subroutine dheqr(A, LDA, N, Q, INFO, IJOB)
Definition: dgmres.f:1579
dhels
subroutine dhels(A, LDA, N, Q, B)
Definition: dgmres.f:1481
msolve
subroutine msolve(n, r, z, nelt, ia, ja, a, isym, rwork, iwork)
Definition: msolve.f:22
drlcal
subroutine drlcal(N, KMP, LL, MAXL, V, Q, RL, SNORMW, PROD, R0NRM)
Definition: dgmres.f:1161

◆ drlcal()

subroutine drlcal ( integer  N,
integer  KMP,
integer  LL,
integer  MAXL,
double precision, dimension(n,*)  V,
double precision, dimension(*)  Q,
double precision, dimension(n)  RL,
double precision  SNORMW,
double precision  PROD,
double precision  R0NRM 
)
1161 C***BEGIN PROLOGUE DRLCAL
1162 C***SUBSIDIARY
1163 C***PURPOSE Internal routine for DGMRES.
1164 C***LIBRARY SLATEC (SLAP)
1165 C***CATEGORY D2A4, D2B4
1166 C***TYPE DOUBLE PRECISION (SRLCAL-S, DRLCAL-D)
1167 C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
1168 C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
1169 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
1170 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
1171 C Seager, Mark K., (LLNL), seager@llnl.gov
1172 C Lawrence Livermore National Laboratory
1173 C PO Box 808, L-60
1174 C Livermore, CA 94550 (510) 423-3141
1175 C***DESCRIPTION
1176 C This routine calculates the scaled residual RL from the
1177 C V(I)'s.
1178 C *Usage:
1179 C INTEGER N, KMP, LL, MAXL
1180 C DOUBLE PRECISION V(N,LL), Q(2*MAXL), RL(N), SNORMW, PROD, R0NORM
1181 C
1182 C CALL DRLCAL(N, KMP, LL, MAXL, V, Q, RL, SNORMW, PROD, R0NRM)
1183 C
1184 C *Arguments:
1185 C N :IN Integer
1186 C The order of the matrix A, and the lengths
1187 C of the vectors SR, SZ, R0 and Z.
1188 C KMP :IN Integer
1189 C The number of previous V vectors the new vector VNEW
1190 C must be made orthogonal to. (KMP .le. MAXL)
1191 C LL :IN Integer
1192 C The current dimension of the Krylov subspace.
1193 C MAXL :IN Integer
1194 C The maximum dimension of the Krylov subspace.
1195 C V :IN Double Precision V(N,LL)
1196 C The N x LL array containing the orthogonal vectors
1197 C V(*,1) to V(*,LL).
1198 C Q :IN Double Precision Q(2*MAXL)
1199 C A double precision array of length 2*MAXL containing the
1200 C components of the Givens rotations used in the QR
1201 C decomposition of HES. It is loaded in DHEQR and used in
1202 C DHELS.
1203 C RL :OUT Double Precision RL(N)
1204 C The residual vector RL. This is either SB*(B-A*XL) if
1205 C not preconditioning or preconditioning on the right,
1206 C or SB*(M-inverse)*(B-A*XL) if preconditioning on the
1207 C left.
1208 C SNORMW :IN Double Precision
1209 C Scale factor.
1210 C PROD :IN Double Precision
1211 C The product s1*s2*...*sl = the product of the sines of the
1212 C Givens rotations used in the QR factorization of
1213 C the Hessenberg matrix HES.
1214 C R0NRM :IN Double Precision
1215 C The scaled norm of initial residual R0.
1216 C
1217 C***SEE ALSO DGMRES
1218 C***ROUTINES CALLED DCOPY, DSCAL
1219 C***REVISION HISTORY (YYMMDD)
1220 C 890404 DATE WRITTEN
1221 C 890404 Previous REVISION DATE
1222 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
1223 C 890922 Numerous changes to prologue to make closer to SLATEC
1224 C standard. (FNF)
1225 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
1226 C 910411 Prologue converted to Version 4.0 format. (BAB)
1227 C 910506 Made subsidiary to DGMRES. (FNF)
1228 C 920511 Added complete declaration section. (WRB)
1229 C***END PROLOGUE DRLCAL
1230 C The following is for optimized compilation on LLNL/LTSS Crays.
1231 CLLL. OPTIMIZE
1232 C .. Scalar Arguments ..
1233  DOUBLE PRECISION prod, r0nrm, snormw
1234  INTEGER kmp, ll, maxl, n
1235 C .. Array Arguments ..
1236  DOUBLE PRECISION q(*), rl(n), v(n,*)
1237 C .. Local Scalars ..
1238  DOUBLE PRECISION c, s, tem
1239  INTEGER i, i2, ip1, k, llm1, llp1
1240 C .. External Subroutines ..
1241  EXTERNAL dcopy, dscal
1242 C***FIRST EXECUTABLE STATEMENT DRLCAL
1243  IF (kmp .EQ. maxl) THEN
1244 C
1245 C calculate RL. Start by copying V(*,1) into RL.
1246 C
1247  CALL dcopy(n, v(1,1), 1, rl, 1)
1248  llm1 = ll - 1
1249  DO 20 i = 1,llm1
1250  ip1 = i + 1
1251  i2 = i*2
1252  s = q(i2)
1253  c = q(i2-1)
1254  DO 10 k = 1,n
1255  rl(k) = s*rl(k) + c*v(k,ip1)
1256  10 CONTINUE
1257  20 CONTINUE
1258  s = q(2*ll)
1259  c = q(2*ll-1)/snormw
1260  llp1 = ll + 1
1261  DO 30 k = 1,n
1262  rl(k) = s*rl(k) + c*v(k,llp1)
1263  30 CONTINUE
1264  ENDIF
1265 C
1266 C When KMP < MAXL, RL vector already partially calculated.
1267 C Scale RL by R0NRM*PROD to obtain the residual RL.
1268 C
1269  tem = r0nrm*prod
1270  CALL dscal(n, tem, rl, 1)
1271  RETURN
1272 C------------- LAST LINE OF DRLCAL FOLLOWS ----------------------------
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
Definition: dgmres.f:1381
dscal
subroutine dscal(n, da, dx, incx)
Definition: dgesv.f:1746

◆ dxlcal()

subroutine dxlcal ( integer  N,
integer  LGMR,
double precision, dimension(n)  X,
double precision, dimension(n)  XL,
double precision, dimension(n)  ZL,
double precision, dimension(maxlp1,*)  HES,
integer  MAXLP1,
double precision, dimension(*)  Q,
double precision, dimension(n,*)  V,
double precision  R0NRM,
double precision, dimension(n)  WK,
double precision, dimension(*)  SZ,
integer  JSCAL,
integer  JPRE,
external  MSOLVE,
integer  NMSL,
double precision, dimension(*)  RPAR,
integer, dimension(*)  IPAR,
integer  NELT,
integer, dimension(nelt)  IA,
integer, dimension(nelt)  JA,
double precision, dimension(nelt)  A,
integer  ISYM 
)
2286 C***BEGIN PROLOGUE DXLCAL
2287 C***SUBSIDIARY
2288 C***PURPOSE Internal routine for DGMRES.
2289 C***LIBRARY SLATEC (SLAP)
2290 C***CATEGORY D2A4, D2B4
2291 C***TYPE DOUBLE PRECISION (SXLCAL-S, DXLCAL-D)
2292 C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
2293 C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
2294 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
2295 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
2296 C Seager, Mark K., (LLNL), seager@llnl.gov
2297 C Lawrence Livermore National Laboratory
2298 C PO Box 808, L-60
2299 C Livermore, CA 94550 (510) 423-3141
2300 C***DESCRIPTION
2301 C This routine computes the solution XL, the current DGMRES
2302 C iterate, given the V(I)'s and the QR factorization of the
2303 C Hessenberg matrix HES. This routine is only called when
2304 C ITOL=11.
2305 C
2306 C *Usage:
2307 C INTEGER N, LGMR, MAXLP1, JSCAL, JPRE, NMSL, IPAR(USER DEFINED)
2308 C INTEGER NELT, IA(NELT), JA(NELT), ISYM
2309 C DOUBLE PRECISION X(N), XL(N), ZL(N), HES(MAXLP1,MAXL), Q(2*MAXL),
2310 C $ V(N,MAXLP1), R0NRM, WK(N), SZ(N),
2311 C $ RPAR(USER DEFINED), A(NELT)
2312 C EXTERNAL MSOLVE
2313 C
2314 C CALL DXLCAL(N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM,
2315 C $ WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR,
2316 C $ NELT, IA, JA, A, ISYM)
2317 C
2318 C *Arguments:
2319 C N :IN Integer
2320 C The order of the matrix A, and the lengths
2321 C of the vectors SR, SZ, R0 and Z.
2322 C LGMR :IN Integer
2323 C The number of iterations performed and
2324 C the current order of the upper Hessenberg
2325 C matrix HES.
2326 C X :IN Double Precision X(N)
2327 C The current approximate solution as of the last restart.
2328 C XL :OUT Double Precision XL(N)
2329 C An array of length N used to hold the approximate
2330 C solution X(L).
2331 C Warning: XL and ZL are the same array in the calling routine.
2332 C ZL :IN Double Precision ZL(N)
2333 C An array of length N used to hold the approximate
2334 C solution Z(L).
2335 C HES :IN Double Precision HES(MAXLP1,MAXL)
2336 C The upper triangular factor of the QR decomposition
2337 C of the (LGMR+1) by LGMR upper Hessenberg matrix whose
2338 C entries are the scaled inner-products of A*V(*,i) and V(*,k).
2339 C MAXLP1 :IN Integer
2340 C MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
2341 C MAXL is the maximum allowable order of the matrix HES.
2342 C Q :IN Double Precision Q(2*MAXL)
2343 C A double precision array of length 2*MAXL containing the
2344 C components of the Givens rotations used in the QR
2345 C decomposition of HES. It is loaded in DHEQR.
2346 C V :IN Double Precision V(N,MAXLP1)
2347 C The N by(LGMR+1) array containing the LGMR
2348 C orthogonal vectors V(*,1) to V(*,LGMR).
2349 C R0NRM :IN Double Precision
2350 C The scaled norm of the initial residual for the
2351 C current call to DPIGMR.
2352 C WK :IN Double Precision WK(N)
2353 C A double precision work array of length N.
2354 C SZ :IN Double Precision SZ(N)
2355 C A vector of length N containing the non-zero
2356 C elements of the diagonal scaling matrix for Z.
2357 C JSCAL :IN Integer
2358 C A flag indicating whether arrays SR and SZ are used.
2359 C JSCAL=0 means SR and SZ are not used and the
2360 C algorithm will perform as if all
2361 C SR(i) = 1 and SZ(i) = 1.
2362 C JSCAL=1 means only SZ is used, and the algorithm
2363 C performs as if all SR(i) = 1.
2364 C JSCAL=2 means only SR is used, and the algorithm
2365 C performs as if all SZ(i) = 1.
2366 C JSCAL=3 means both SR and SZ are used.
2367 C JPRE :IN Integer
2368 C The preconditioner type flag.
2369 C MSOLVE :EXT External.
2370 C Name of the routine which solves a linear system Mz = r for
2371 C z given r with the preconditioning matrix M (M is supplied via
2372 C RPAR and IPAR arrays. The name of the MSOLVE routine must
2373 C be declared external in the calling program. The calling
2374 C sequence to MSOLVE is:
2375 C CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RPAR, IPAR)
2376 C Where N is the number of unknowns, R is the right-hand side
2377 C vector and Z is the solution upon return. NELT, IA, JA, A and
2378 C ISYM are defined as below. RPAR is a double precision array
2379 C that can be used to pass necessary preconditioning information
2380 C and/or workspace to MSOLVE. IPAR is an integer work array
2381 C for the same purpose as RPAR.
2382 C NMSL :IN Integer
2383 C The number of calls to MSOLVE.
2384 C RPAR :IN Double Precision RPAR(USER DEFINED)
2385 C Double Precision workspace passed directly to the MSOLVE
2386 C routine.
2387 C IPAR :IN Integer IPAR(USER DEFINED)
2388 C Integer workspace passed directly to the MSOLVE routine.
2389 C NELT :IN Integer
2390 C The length of arrays IA, JA and A.
2391 C IA :IN Integer IA(NELT)
2392 C An integer array of length NELT containing matrix data.
2393 C It is passed directly to the MATVEC and MSOLVE routines.
2394 C JA :IN Integer JA(NELT)
2395 C An integer array of length NELT containing matrix data.
2396 C It is passed directly to the MATVEC and MSOLVE routines.
2397 C A :IN Double Precision A(NELT)
2398 C A double precision array of length NELT containing matrix
2399 C data.
2400 C It is passed directly to the MATVEC and MSOLVE routines.
2401 C ISYM :IN Integer
2402 C A flag to indicate symmetric matrix storage.
2403 C If ISYM=0, all non-zero entries of the matrix are
2404 C stored. If ISYM=1, the matrix is symmetric and
2405 C only the upper or lower triangular part is stored.
2406 C
2407 C***SEE ALSO DGMRES
2408 C***ROUTINES CALLED DAXPY, DCOPY, DHELS
2409 C***REVISION HISTORY (YYMMDD)
2410 C 890404 DATE WRITTEN
2411 C 890404 Previous REVISION DATE
2412 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
2413 C 890922 Numerous changes to prologue to make closer to SLATEC
2414 C standard. (FNF)
2415 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
2416 C 910411 Prologue converted to Version 4.0 format. (BAB)
2417 C 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF)
2418 C 910506 Made subsidiary to DGMRES. (FNF)
2419 C 920511 Added complete declaration section. (WRB)
2420 C***END PROLOGUE DXLCAL
2421 C The following is for optimized compilation on LLNL/LTSS Crays.
2422 CLLL. OPTIMIZE
2423 C .. Scalar Arguments ..
2424  DOUBLE PRECISION r0nrm
2425  INTEGER isym, jpre, jscal, lgmr, maxlp1, n, nelt, nmsl
2426 C .. Array Arguments ..
2427  DOUBLE PRECISION a(nelt), hes(maxlp1,*), q(*), rpar(*), sz(*),
2428  + v(n,*), wk(n), x(n), xl(n), zl(n)
2429  INTEGER ia(nelt), ipar(*), ja(nelt)
2430 C .. Subroutine Arguments ..
2431  EXTERNAL msolve
2432 C .. Local Scalars ..
2433  INTEGER i, k, ll, llp1
2434 C .. External Subroutines ..
2435  EXTERNAL daxpy, dcopy, dhels
2436 C***FIRST EXECUTABLE STATEMENT DXLCAL
2437  ll = lgmr
2438  llp1 = ll + 1
2439  DO 10 k = 1,llp1
2440  wk(k) = 0
2441  10 CONTINUE
2442  wk(1) = r0nrm
2443  CALL dhels(hes, maxlp1, ll, q, wk)
2444  DO 20 k = 1,n
2445  zl(k) = 0
2446  20 CONTINUE
2447  DO 30 i = 1,ll
2448  CALL daxpy(n, wk(i), v(1,i), 1, zl, 1)
2449  30 CONTINUE
2450  IF ((jscal .EQ. 1) .OR.(jscal .EQ. 3)) THEN
2451  DO 40 k = 1,n
2452  zl(k) = zl(k)/sz(k)
2453  40 CONTINUE
2454  ENDIF
2455  IF (jpre .GT. 0) THEN
2456  CALL dcopy(n, zl, 1, wk, 1)
2457  CALL msolve(n, wk, zl, nelt, ia, ja, a, isym, rpar, ipar)
2458  nmsl = nmsl + 1
2459  ENDIF
2460 C calculate XL from X and ZL.
2461  DO 50 k = 1,n
2462  xl(k) = x(k) + zl(k)
2463  50 CONTINUE
2464  RETURN
2465 C------------- LAST LINE OF DXLCAL FOLLOWS ----------------------------
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
Definition: dgmres.f:1381
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
Definition: dgmres.f:1276
dhels
subroutine dhels(A, LDA, N, Q, B)
Definition: dgmres.f:1481
msolve
subroutine msolve(n, r, z, nelt, ia, ja, a, isym, rwork, iwork)
Definition: msolve.f:22

◆ isdgmr()

integer function isdgmr ( integer  N,
double precision, dimension(*)  B,
double precision, dimension(*)  X,
double precision, dimension(*)  XL,
integer  NELT,
integer, dimension(*)  IA,
integer, dimension(*)  JA,
double precision, dimension(*)  A,
integer  ISYM,
external  MSOLVE,
integer  NMSL,
integer  ITOL,
double precision  TOL,
integer  ITMAX,
integer  ITER,
double precision  ERR,
integer  IUNIT,
double precision, dimension(*)  R,
double precision, dimension(*)  Z,
double precision, dimension(*)  DZ,
double precision, dimension(*)  RWORK,
integer, dimension(*)  IWORK,
double precision  RNRM,
double precision  BNRM,
double precision, dimension(*)  SB,
double precision, dimension(*)  SX,
integer  JSCAL,
integer  KMP,
integer  LGMR,
integer  MAXL,
integer  MAXLP1,
double precision, dimension(n,*)  V,
double precision, dimension(*)  Q,
double precision  SNORMW,
double precision  PROD,
double precision  R0NRM,
double precision, dimension(maxlp1, maxl)  HES,
integer  JPRE 
)
1760 C***BEGIN PROLOGUE ISDGMR
1761 C***SUBSIDIARY
1762 C***PURPOSE Generalized Minimum Residual Stop Test.
1763 C This routine calculates the stop test for the Generalized
1764 C Minimum RESidual (GMRES) iteration scheme. It returns a
1765 C non-zero if the error estimate (the type of which is
1766 C determined by ITOL) is less than the user specified
1767 C tolerance TOL.
1768 C***LIBRARY SLATEC (SLAP)
1769 C***CATEGORY D2A4, D2B4
1770 C***TYPE DOUBLE PRECISION (ISSGMR-S, ISDGMR-D)
1771 C***KEYWORDS GMRES, LINEAR SYSTEM, SLAP, SPARSE, STOP TEST
1772 C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
1773 C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
1774 C Seager, Mark K., (LLNL), seager@llnl.gov
1775 C Lawrence Livermore National Laboratory
1776 C PO Box 808, L-60
1777 C Livermore, CA 94550 (510) 423-3141
1778 C***DESCRIPTION
1779 C
1780 C *Usage:
1781 C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, NMSL, ITOL
1782 C INTEGER ITMAX, ITER, IUNIT, IWORK(USER DEFINED), JSCAL
1783 C INTEGER KMP, LGMR, MAXL, MAXLP1, JPRE
1784 C DOUBLE PRECISION B(N), X(N), XL(MAXL), A(NELT), TOL, ERR,
1785 C $ R(N), Z(N), DZ(N), RWORK(USER DEFINED),
1786 C $ RNRM, BNRM, SB(N), SX(N), V(N,MAXLP1),
1787 C $ Q(2*MAXL), SNORMW, PROD, R0NRM,
1788 C $ HES(MAXLP1,MAXL)
1789 C EXTERNAL MSOLVE
1790 C
1791 C IF (ISDGMR(N, B, X, XL, NELT, IA, JA, A, ISYM, MSOLVE,
1792 C $ NMSL, ITOL, TOL, ITMAX, ITER, ERR, IUNIT, R, Z, DZ,
1793 C $ RWORK, IWORK, RNRM, BNRM, SB, SX, JSCAL,
1794 C $ KMP, LGMR, MAXL, MAXLP1, V, Q, SNORMW, PROD, R0NRM,
1795 C $ HES, JPRE) .NE. 0) THEN ITERATION DONE
1796 C
1797 C *Arguments:
1798 C N :IN Integer.
1799 C Order of the Matrix.
1800 C B :IN Double Precision B(N).
1801 C Right-hand-side vector.
1802 C X :IN Double Precision X(N).
1803 C Approximate solution vector as of the last restart.
1804 C XL :OUT Double Precision XL(N)
1805 C An array of length N used to hold the approximate
1806 C solution as of the current iteration. Only computed by
1807 C this routine when ITOL=11.
1808 C NELT :IN Integer.
1809 C Number of Non-Zeros stored in A.
1810 C IA :IN Integer IA(NELT).
1811 C JA :IN Integer JA(NELT).
1812 C A :IN Double Precision A(NELT).
1813 C These arrays contain the matrix data structure for A.
1814 C It could take any form. See "Description", in the DGMRES,
1815 C DSLUGM and DSDGMR routines for more details.
1816 C ISYM :IN Integer.
1817 C Flag to indicate symmetric storage format.
1818 C If ISYM=0, all non-zero entries of the matrix are stored.
1819 C If ISYM=1, the matrix is symmetric, and only the upper
1820 C or lower triangle of the matrix is stored.
1821 C MSOLVE :EXT External.
1822 C Name of a routine which solves a linear system Mz = r for z
1823 C given r with the preconditioning matrix M (M is supplied via
1824 C RWORK and IWORK arrays. The name of the MSOLVE routine must
1825 C be declared external in the calling program. The calling
1826 C sequence to MSOLVE is:
1827 C CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
1828 C Where N is the number of unknowns, R is the right-hand side
1829 C vector and Z is the solution upon return. NELT, IA, JA, A and
1830 C ISYM are defined as above. RWORK is a double precision array
1831 C that can be used to pass necessary preconditioning information
1832 C and/or workspace to MSOLVE. IWORK is an integer work array
1833 C for the same purpose as RWORK.
1834 C NMSL :INOUT Integer.
1835 C A counter for the number of calls to MSOLVE.
1836 C ITOL :IN Integer.
1837 C Flag to indicate the type of convergence criterion used.
1838 C ITOL=0 Means the iteration stops when the test described
1839 C below on the residual RL is satisfied. This is
1840 C the "Natural Stopping Criteria" for this routine.
1841 C Other values of ITOL cause extra, otherwise
1842 C unnecessary, computation per iteration and are
1843 C therefore much less efficient.
1844 C ITOL=1 Means the iteration stops when the first test
1845 C described below on the residual RL is satisfied,
1846 C and there is either right or no preconditioning
1847 C being used.
1848 C ITOL=2 Implies that the user is using left
1849 C preconditioning, and the second stopping criterion
1850 C below is used.
1851 C ITOL=3 Means the iteration stops when the third test
1852 C described below on Minv*Residual is satisfied, and
1853 C there is either left or no preconditioning begin
1854 C used.
1855 C ITOL=11 is often useful for checking and comparing
1856 C different routines. For this case, the user must
1857 C supply the "exact" solution or a very accurate
1858 C approximation (one with an error much less than
1859 C TOL) through a common block,
1860 C COMMON /DSLBLK/ SOLN( )
1861 C If ITOL=11, iteration stops when the 2-norm of the
1862 C difference between the iterative approximation and
1863 C the user-supplied solution divided by the 2-norm
1864 C of the user-supplied solution is less than TOL.
1865 C Note that this requires the user to set up the
1866 C "COMMON /DSLBLK/ SOLN(LENGTH)" in the calling
1867 C routine. The routine with this declaration should
1868 C be loaded before the stop test so that the correct
1869 C length is used by the loader. This procedure is
1870 C not standard Fortran and may not work correctly on
1871 C your system (although it has worked on every
1872 C system the authors have tried). If ITOL is not 11
1873 C then this common block is indeed standard Fortran.
1874 C TOL :IN Double Precision.
1875 C Convergence criterion, as described above.
1876 C ITMAX :IN Integer.
1877 C Maximum number of iterations.
1878 C ITER :IN Integer.
1879 C The iteration for which to check for convergence.
1880 C ERR :OUT Double Precision.
1881 C Error estimate of error in final approximate solution, as
1882 C defined by ITOL. Letting norm() denote the Euclidean
1883 C norm, ERR is defined as follows..
1884 C
1885 C If ITOL=0, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
1886 C for right or no preconditioning, and
1887 C ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
1888 C norm(SB*(M-inverse)*B),
1889 C for left preconditioning.
1890 C If ITOL=1, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
1891 C since right or no preconditioning
1892 C being used.
1893 C If ITOL=2, then ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
1894 C norm(SB*(M-inverse)*B),
1895 C since left preconditioning is being
1896 C used.
1897 C If ITOL=3, then ERR = Max |(Minv*(B-A*X(L)))(i)/x(i)|
1898 C i=1,n
1899 C If ITOL=11, then ERR = norm(SB*(X(L)-SOLN))/norm(SB*SOLN).
1900 C IUNIT :IN Integer.
1901 C Unit number on which to write the error at each iteration,
1902 C if this is desired for monitoring convergence. If unit
1903 C number is 0, no writing will occur.
1904 C R :INOUT Double Precision R(N).
1905 C Work array used in calling routine. It contains
1906 C information necessary to compute the residual RL = B-A*XL.
1907 C Z :WORK Double Precision Z(N).
1908 C Workspace used to hold the pseudo-residual M z = r.
1909 C DZ :WORK Double Precision DZ(N).
1910 C Workspace used to hold temporary vector(s).
1911 C RWORK :WORK Double Precision RWORK(USER DEFINED).
1912 C Double Precision array that can be used by MSOLVE.
1913 C IWORK :WORK Integer IWORK(USER DEFINED).
1914 C Integer array that can be used by MSOLVE.
1915 C RNRM :IN Double Precision.
1916 C Norm of the current residual. Type of norm depends on ITOL.
1917 C BNRM :IN Double Precision.
1918 C Norm of the right hand side. Type of norm depends on ITOL.
1919 C SB :IN Double Precision SB(N).
1920 C Scaling vector for B.
1921 C SX :IN Double Precision SX(N).
1922 C Scaling vector for X.
1923 C JSCAL :IN Integer.
1924 C Flag indicating if scaling arrays SB and SX are being
1925 C used in the calling routine DPIGMR.
1926 C JSCAL=0 means SB and SX are not used and the
1927 C algorithm will perform as if all
1928 C SB(i) = 1 and SX(i) = 1.
1929 C JSCAL=1 means only SX is used, and the algorithm
1930 C performs as if all SB(i) = 1.
1931 C JSCAL=2 means only SB is used, and the algorithm
1932 C performs as if all SX(i) = 1.
1933 C JSCAL=3 means both SB and SX are used.
1934 C KMP :IN Integer
1935 C The number of previous vectors the new vector VNEW
1936 C must be made orthogonal to. (KMP .le. MAXL)
1937 C LGMR :IN Integer
1938 C The number of GMRES iterations performed on the current call
1939 C to DPIGMR (i.e., # iterations since the last restart) and
1940 C the current order of the upper Hessenberg
1941 C matrix HES.
1942 C MAXL :IN Integer
1943 C The maximum allowable order of the matrix H.
1944 C MAXLP1 :IN Integer
1945 C MAXPL1 = MAXL + 1, used for dynamic dimensioning of HES.
1946 C V :IN Double Precision V(N,MAXLP1)
1947 C The N by (LGMR+1) array containing the LGMR
1948 C orthogonal vectors V(*,1) to V(*,LGMR).
1949 C Q :IN Double Precision Q(2*MAXL)
1950 C A double precision array of length 2*MAXL containing the
1951 C components of the Givens rotations used in the QR
1952 C decomposition of HES.
1953 C SNORMW :IN Double Precision
1954 C A scalar containing the scaled norm of VNEW before it
1955 C is renormalized in DPIGMR.
1956 C PROD :IN Double Precision
1957 C The product s1*s2*...*sl = the product of the sines of the
1958 C Givens rotations used in the QR factorization of the
1959 C Hessenberg matrix HES.
1960 C R0NRM :IN Double Precision
1961 C The scaled norm of initial residual R0.
1962 C HES :IN Double Precision HES(MAXLP1,MAXL)
1963 C The upper triangular factor of the QR decomposition
1964 C of the (LGMR+1) by LGMR upper Hessenberg matrix whose
1965 C entries are the scaled inner-products of A*V(*,I)
1966 C and V(*,K).
1967 C JPRE :IN Integer
1968 C Preconditioner type flag.
1969 C (See description of IGWK(4) in DGMRES.)
1970 C
1971 C *Description
1972 C When using the GMRES solver, the preferred value for ITOL
1973 C is 0. This is due to the fact that when ITOL=0 the norm of
1974 C the residual required in the stopping test is obtained for
1975 C free, since this value is already calculated in the GMRES
1976 C algorithm. The variable RNRM contains the appropriate
1977 C norm, which is equal to norm(SB*(RL - A*XL)) when right or
1978 C no preconditioning is being performed, and equal to
1979 C norm(SB*Minv*(RL - A*XL)) when using left preconditioning.
1980 C Here, norm() is the Euclidean norm. Nonzero values of ITOL
1981 C require additional work to calculate the actual scaled
1982 C residual or its scaled/preconditioned form, and/or the
1983 C approximate solution XL. Hence, these values of ITOL will
1984 C not be as efficient as ITOL=0.
1985 C
1986 C *Cautions:
1987 C This routine will attempt to write to the Fortran logical output
1988 C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that
1989 C this logical unit is attached to a file or terminal before calling
1990 C this routine with a non-zero value for IUNIT. This routine does
1991 C not check for the validity of a non-zero IUNIT unit number.
1992 C
1993 C This routine does not verify that ITOL has a valid value.
1994 C The calling routine should make such a test before calling
1995 C ISDGMR, as is done in DGMRES.
1996 C
1997 C***SEE ALSO DGMRES
1998 C***ROUTINES CALLED D1MACH, DCOPY, DNRM2, DRLCAL, DSCAL, DXLCAL
1999 C***COMMON BLOCKS DSLBLK
2000 C***REVISION HISTORY (YYMMDD)
2001 C 890404 DATE WRITTEN
2002 C 890404 Previous REVISION DATE
2003 C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
2004 C 890922 Numerous changes to prologue to make closer to SLATEC
2005 C standard. (FNF)
2006 C 890929 Numerous changes to reduce SP/DP differences. (FNF)
2007 C 910411 Prologue converted to Version 4.0 format. (BAB)
2008 C 910502 Corrected conversion errors, etc. (FNF)
2009 C 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF)
2010 C 910506 Made subsidiary to DGMRES. (FNF)
2011 C 920407 COMMON BLOCK renamed DSLBLK. (WRB)
2012 C 920511 Added complete declaration section. (WRB)
2013 C 921026 Corrected D to E in output format. (FNF)
2014 C 921113 Corrected C***CATEGORY line. (FNF)
2015 C***END PROLOGUE ISDGMR
2016 C .. Scalar Arguments ..
2017  DOUBLE PRECISION bnrm, err, prod, r0nrm, rnrm, snormw, tol
2018  INTEGER isym, iter, itmax, itol, iunit, jpre, jscal, kmp, lgmr,
2019  + maxl, maxlp1, n, nelt, nmsl
2020 C .. Array Arguments ..
2021  DOUBLE PRECISION a(*), b(*), dz(*), hes(maxlp1, maxl), q(*), r(*),
2022  + rwork(*), sb(*), sx(*), v(n,*), x(*), xl(*), z(*)
2023  INTEGER ia(*), iwork(*), ja(*)
2024 C .. Subroutine Arguments ..
2025  EXTERNAL msolve
2026 C .. Arrays in Common ..
2027  DOUBLE PRECISION soln(1)
2028 C .. Local Scalars ..
2029  DOUBLE PRECISION dxnrm, fuzz, rat, ratmax, solnrm, tem
2030  INTEGER i, ielmax
2031 C .. External Functions ..
2032  DOUBLE PRECISION d1mach, dnrm2
2033  EXTERNAL d1mach, dnrm2
2034 C .. External Subroutines ..
2035  EXTERNAL dcopy, drlcal, dscal, dxlcal
2036 C .. Intrinsic Functions ..
2037  INTRINSIC abs, max, sqrt
2038 C .. Common blocks ..
2039  COMMON /dslblk/ soln
2040 C .. Save statement ..
2041  SAVE solnrm
2042 C***FIRST EXECUTABLE STATEMENT ISDGMR
2043  isdgmr = 0
2044  IF ( itol.EQ.0 ) THEN
2045 C
2046 C Use input from DPIGMR to determine if stop conditions are met.
2047 C
2048  err = rnrm/bnrm
2049  ENDIF
2050  IF ( (itol.GT.0) .AND. (itol.LE.3) ) THEN
2051 C
2052 C Use DRLCAL to calculate the scaled residual vector.
2053 C Store answer in R.
2054 C
2055  IF ( lgmr.NE.0 ) CALL drlcal(n, kmp, lgmr, maxl, v, q, r,
2056  $ snormw, prod, r0nrm)
2057  IF ( itol.LE.2 ) THEN
2058 C err = ||Residual||/||RightHandSide||(2-Norms).
2059  err = dnrm2(n, r, 1)/bnrm
2060 C
2061 C Unscale R by R0NRM*PROD when KMP < MAXL.
2062 C
2063  IF ( (kmp.LT.maxl) .AND. (lgmr.NE.0) ) THEN
2064  tem = 1.0d0/(r0nrm*prod)
2065  CALL dscal(n, tem, r, 1)
2066  ENDIF
2067  ELSEIF ( itol.EQ.3 ) THEN
2068 C err = Max |(Minv*Residual)(i)/x(i)|
2069 C When JPRE .lt. 0, R already contains Minv*Residual.
2070  IF ( jpre.GT.0 ) THEN
2071  CALL msolve(n, r, dz, nelt, ia, ja, a, isym, rwork,
2072  $ iwork)
2073  nmsl = nmsl + 1
2074  ENDIF
2075 C
2076 C Unscale R by R0NRM*PROD when KMP < MAXL.
2077 C
2078  IF ( (kmp.LT.maxl) .AND. (lgmr.NE.0) ) THEN
2079  tem = 1.0d0/(r0nrm*prod)
2080  CALL dscal(n, tem, r, 1)
2081  ENDIF
2082 C
2083  fuzz = d1mach(1)
2084  ielmax = 1
2085  ratmax = abs(dz(1))/max(abs(x(1)),fuzz)
2086  DO 25 i = 2, n
2087  rat = abs(dz(i))/max(abs(x(i)),fuzz)
2088  IF( rat.GT.ratmax ) THEN
2089  ielmax = i
2090  ratmax = rat
2091  ENDIF
2092  25 CONTINUE
2093  err = ratmax
2094  IF( ratmax.LE.tol ) isdgmr = 1
2095  IF( iunit.GT.0 ) WRITE(iunit,1020) iter, ielmax, ratmax
2096  RETURN
2097  ENDIF
2098  ENDIF
2099  IF ( itol.EQ.11 ) THEN
2100 C
2101 C Use DXLCAL to calculate the approximate solution XL.
2102 C
2103  IF ( (lgmr.NE.0) .AND. (iter.GT.0) ) THEN
2104  CALL dxlcal(n, lgmr, x, xl, xl, hes, maxlp1, q, v, r0nrm,
2105  $ dz, sx, jscal, jpre, msolve, nmsl, rwork, iwork,
2106  $ nelt, ia, ja, a, isym)
2107  ELSEIF ( iter.EQ.0 ) THEN
2108 C Copy X to XL to check if initial guess is good enough.
2109  CALL dcopy(n, x, 1, xl, 1)
2110  ELSE
2111 C Return since this is the first call to DPIGMR on a restart.
2112  RETURN
2113  ENDIF
2114 C
2115  IF ((jscal .EQ. 0) .OR.(jscal .EQ. 2)) THEN
2116 C err = ||x-TrueSolution||/||TrueSolution||(2-Norms).
2117  IF ( iter.EQ.0 ) solnrm = dnrm2(n, soln, 1)
2118  DO 30 i = 1, n
2119  dz(i) = xl(i) - soln(i)
2120  30 CONTINUE
2121  err = dnrm2(n, dz, 1)/solnrm
2122  ELSE
2123  IF (iter .EQ. 0) THEN
2124  solnrm = 0
2125  DO 40 i = 1,n
2126  solnrm = solnrm + (sx(i)*soln(i))**2
2127  40 CONTINUE
2128  solnrm = sqrt(solnrm)
2129  ENDIF
2130  dxnrm = 0
2131  DO 50 i = 1,n
2132  dxnrm = dxnrm + (sx(i)*(xl(i)-soln(i)))**2
2133  50 CONTINUE
2134  dxnrm = sqrt(dxnrm)
2135 C err = ||SX*(x-TrueSolution)||/||SX*TrueSolution|| (2-Norms).
2136  err = dxnrm/solnrm
2137  ENDIF
2138  ENDIF
2139 C
2140  IF( iunit.NE.0 ) THEN
2141  IF( iter.EQ.0 ) THEN
2142  WRITE(iunit,1000) n, itol, maxl, kmp
2143  ENDIF
2144  WRITE(iunit,1010) iter, rnrm/bnrm, err
2145  ENDIF
2146  IF ( err.LE.tol ) isdgmr = 1
2147 C
2148  RETURN
2149  1000 FORMAT(' Generalized Minimum Residual(',i3,i3,') for ',
2150  $ 'N, ITOL = ',i5, i5,
2151  $ /' ITER',' Natural Err Est',' Error Estimate')
2152  1010 FORMAT(1x,i4,1x,d16.7,1x,d16.7)
2153  1020 FORMAT(1x,' ITER = ',i5, ' IELMAX = ',i5,
2154  $ ' |R(IELMAX)/X(IELMAX)| = ',d12.5)
2155 C------------- LAST LINE OF ISDGMR FOLLOWS ----------------------------
max
#define max(a, b)
Definition: cascade.c:32
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
Definition: dgmres.f:1381
d1mach
double precision function d1mach(I)
Definition: ddeabm.f:2012
dnrm2
double precision function dnrm2(N, DX, INCX)
Definition: dgmres.f:559
dxlcal
subroutine dxlcal(N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM, WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR, NELT, IA, JA, A, ISYM)
Definition: dgmres.f:2286
dscal
subroutine dscal(n, da, dx, incx)
Definition: dgesv.f:1746
isdgmr
integer function isdgmr(N, B, X, XL, NELT, IA, JA, A, ISYM, MSOLVE, NMSL, ITOL, TOL, ITMAX, ITER, ERR, IUNIT, R, Z, DZ, RWORK, IWORK, RNRM, BNRM, SB, SX, JSCAL, KMP, LGMR, MAXL, MAXLP1, V, Q, SNORMW, PROD, R0NRM, HES, JPRE)
Definition: dgmres.f:1760
msolve
subroutine msolve(n, r, z, nelt, ia, ja, a, isym, rwork, iwork)
Definition: msolve.f:22
drlcal
subroutine drlcal(N, KMP, LL, MAXL, V, Q, RL, SNORMW, PROD, R0NRM)
Definition: dgmres.f:1161
Hosted by OpenAircraft.com, (Michigan UAV, LLC)