GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
datv.f
Go to the documentation of this file.
1 C Work performed under the auspices of the U.S. Department of Energy
2 C by Lawrence Livermore National Laboratory under contract number
3 C W-7405-Eng-48.
4 C
5  SUBROUTINE datv (NEQ, Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, RES,
6  * IRES, PSOL, Z, VTEM, WP, IWP, CJ, EPLIN, IER, NRE, NPSL,
7  * RPAR,IPAR)
8 C
9 C***BEGIN PROLOGUE DATV
10 C***DATE WRITTEN 890101 (YYMMDD)
11 C***REVISION DATE 900926 (YYMMDD)
12 C
13 C
14 C-----------------------------------------------------------------------
15 C***DESCRIPTION
16 C
17 C This routine computes the product
18 C
19 C Z = (D-inverse)*(P-inverse)*(dF/dY)*(D*V),
20 C
21 C where F(Y) = G(T, Y, CJ*(Y-A)), CJ is a scalar proportional to 1/H,
22 C and A involves the past history of Y. The quantity CJ*(Y-A) is
23 C an approximation to the first derivative of Y and is stored
24 C in the array YPRIME. Note that dF/dY = dG/dY + CJ*dG/dYPRIME.
25 C
26 C D is a diagonal scaling matrix, and P is the left preconditioning
27 C matrix. V is assumed to have L2 norm equal to 1.
28 C The product is stored in Z and is computed by means of a
29 C difference quotient, a call to RES, and one call to PSOL.
30 C
31 C On entry
32 C
33 C NEQ = Problem size, passed to RES and PSOL.
34 C
35 C Y = Array containing current dependent variable vector.
36 C
37 C YPRIME = Array containing current first derivative of y.
38 C
39 C SAVR = Array containing current value of G(T,Y,YPRIME).
40 C
41 C V = Real array of length NEQ (can be the same array as Z).
42 C
43 C WGHT = Array of length NEQ containing scale factors.
44 C 1/WGHT(I) are the diagonal elements of the matrix D.
45 C
46 C YPTEM = Work array of length NEQ.
47 C
48 C VTEM = Work array of length NEQ used to store the
49 C unscaled version of V.
50 C
51 C WP = Real work array used by preconditioner PSOL.
52 C
53 C IWP = Integer work array used by preconditioner PSOL.
54 C
55 C CJ = Scalar proportional to current value of
56 C 1/(step size H).
57 C
58 C
59 C On return
60 C
61 C Z = Array of length NEQ containing desired scaled
62 C matrix-vector product.
63 C
64 C IRES = Error flag from RES.
65 C
66 C IER = Error flag from PSOL.
67 C
68 C NRE = The number of calls to RES.
69 C
70 C NPSL = The number of calls to PSOL.
71 C
72 C-----------------------------------------------------------------------
73 C***ROUTINES CALLED
74 C RES, PSOL
75 C
76 C***END PROLOGUE DATV
77 C
78  INTEGER NEQ, IRES, IWP, IER, NRE, NPSL, IPAR
79  DOUBLE PRECISION Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, Z, VTEM,
80  1 wp, cj, rpar
81  dimension y(*), yprime(*), savr(*), v(*), wght(*), yptem(*),
82  1 z(*), vtem(*), wp(*), iwp(*), rpar(*), ipar(*)
83  INTEGER I
84  DOUBLE PRECISION EPLIN
85  EXTERNAL res, psol
86 C
87  ires = 0
88 C-----------------------------------------------------------------------
89 C Set VTEM = D * V.
90 C-----------------------------------------------------------------------
91  DO 10 i = 1,neq
92  10 vtem(i) = v(i)/wght(i)
93  ier = 0
94 C-----------------------------------------------------------------------
95 C Store Y in Z and increment Z by VTEM.
96 C Store YPRIME in YPTEM and increment YPTEM by VTEM*CJ.
97 C-----------------------------------------------------------------------
98  DO 20 i = 1,neq
99  yptem(i) = yprime(i) + vtem(i)*cj
100  20 z(i) = y(i) + vtem(i)
101 C-----------------------------------------------------------------------
102 C Call RES with incremented Y, YPRIME arguments
103 C stored in Z, YPTEM. VTEM is overwritten with new residual.
104 C-----------------------------------------------------------------------
105  CONTINUE
106  CALL res(tn,z,yptem,cj,vtem,ires,rpar,ipar)
107  nre = nre + 1
108  IF (ires .LT. 0) RETURN
109 C-----------------------------------------------------------------------
110 C Set Z = (dF/dY) * VBAR using difference quotient.
111 C (VBAR is old value of VTEM before calling RES)
112 C-----------------------------------------------------------------------
113  DO 70 i = 1,neq
114  70 z(i) = vtem(i) - savr(i)
115 C-----------------------------------------------------------------------
116 C Apply inverse of left preconditioner to Z.
117 C-----------------------------------------------------------------------
118  CALL psol (neq, tn, y, yprime, savr, yptem, cj, wght, wp, iwp,
119  1 z, eplin, ier, rpar, ipar)
120  npsl = npsl + 1
121  IF (ier .NE. 0) RETURN
122 C-----------------------------------------------------------------------
123 C Apply D-inverse to Z and return.
124 C-----------------------------------------------------------------------
125  DO 90 i = 1,neq
126  90 z(i) = z(i)*wght(i)
127  RETURN
128 C
129 C------END OF SUBROUTINE DATV-------------------------------------------
130  END
datv
subroutine datv(NEQ, Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, RES, IRES, PSOL, Z, VTEM, WP, IWP, CJ, EPLIN, IER, NRE, NPSL, RPAR, IPAR)
Definition: datv.f:8