GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
std::string dimension(void) const
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:5