← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/integ/prja.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/integ/prja.f
 
1
C/MEMBR ADD NAME=PRJA,SSI=0
 
2
      subroutine prja (neq, y, yh, nyh, ewt, ftem, savf, wm, iwm,
 
3
     1   f, jac)
 
4
clll. optimize
 
5
      external f, jac
 
6
      integer neq, nyh, iwm
 
7
      integer iownd, iowns,
 
8
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
 
9
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
 
10
      integer iownd2, iowns2, jtyp, mused, mxordn, mxords
 
11
      integer i, i1, i2, ier, ii, j, j1, jj, lenp,
 
12
     1   mba, mband, meb1, meband, ml, ml3, mu
 
13
      double precision y, yh, ewt, ftem, savf, wm
 
14
      double precision rownd, rowns,
 
15
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
 
16
      double precision rownd2, rowns2, pdnorm
 
17
      double precision con, fac, hl0, r, r0, srur, yi, yj, yjj,
 
18
     1   vmnorm, fnorm, bnorm
 
19
      dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
 
20
     1   wm(*), iwm(*)
 
21
      integer         iero
 
22
      common /ierode/ iero
 
23
      common /ls0001/ rownd, rowns(209),
 
24
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
 
25
     3   iownd(14), iowns(6),
 
26
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
 
27
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
 
28
      common /lsa001/ rownd2, rowns2(20), pdnorm,
 
29
     1   iownd2(3), iowns2(2), jtyp, mused, mxordn, mxords
 
30
c-----------------------------------------------------------------------
 
31
c%purpose
 
32
c prja is called by stoda to compute and process the matrix
 
33
c p = i - h*el(1)*j , where j is an approximation to the jacobian.
 
34
c here j is computed by the user-supplied routine jac if
 
35
c miter = 1 or 4 or by finite differencing if miter = 2 or 5.
 
36
c j, scaled by -h*el(1), is stored in wm.  then the norm of j (the
 
37
c matrix norm consistent with the weighted max-norm on vectors given
 
38
c by vmnorm) is computed, and j is overwritten by p.  p is then
 
39
c subjected to lu decomposition in preparation for later solution
 
40
c of linear systems with p as coefficient matrix. this is done
 
41
c by dgefa if miter = 1 or 2, and by dgbfa if miter = 4 or 5.
 
42
c
 
43
c%additional  parameters
 
44
c in addition to variables described previously, communication
 
45
c with prja uses the following..
 
46
c y     = array containing predicted values on entry.
 
47
c ftem  = work array of length n (acor in stoda).
 
48
c savf  = array containing f evaluated at predicted y.
 
49
c wm    = real work space for matrices.  on output it contains the
 
50
c         lu decomposition of p.
 
51
c         storage of matrix elements starts at wm(3).
 
52
c         wm also contains the following matrix-related data..
 
53
c         wm(1) = sqrt(uround), used in numerical jacobian increments.
 
54
c iwm   = integer work space containing pivot information, starting at
 
55
c         iwm(21).   iwm also contains the band parameters
 
56
c         ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
 
57
c el0   = el(1) (input).
 
58
c pdnorm= norm of jacobian matrix. (output).
 
59
c ierpj = output error flag,  = 0 if no trouble, .gt. 0 if
 
60
c         p matrix found to be singular.
 
61
c jcur  = output flag = 1 to indicate that the jacobian matrix
 
62
c         (or approximation) is now current.
 
63
c this routine also uses the common variables el0, h, tn, uround,
 
64
c miter, n, nfe, and nje.
 
65
c!
 
66
c-----------------------------------------------------------------------
 
67
      nje = nje + 1
 
68
      ierpj = 0
 
69
      jcur = 1
 
70
      hl0 = h*el0
 
71
      go to (100, 200, 300, 400, 500), miter
 
72
c if miter = 1, call jac and multiply by scalar. -----------------------
 
73
 100  lenp = n*n
 
74
      do 110 i = 1,lenp
 
75
 110    wm(i+2) = 0.0d+0
 
76
      call jac (neq, tn, y, 0, 0, wm(3), n)
 
77
      if(iero.gt.0) return
 
78
      con = -hl0
 
79
      do 120 i = 1,lenp
 
80
 120    wm(i+2) = wm(i+2)*con
 
81
      go to 240
 
82
c if miter = 2, make n calls to f to approximate j. --------------------
 
83
 200  fac = vmnorm (n, savf, ewt)
 
84
      r0 = 1000.0d+0*abs(h)*uround*dble(n)*fac
 
85
      if (r0 .eq. 0.0d+0) r0 = 1.0d+0
 
86
      srur = wm(1)
 
87
      j1 = 2
 
88
      do 230 j = 1,n
 
89
        yj = y(j)
 
90
        r = max(srur*abs(yj),r0/ewt(j))
 
91
        y(j) = y(j) + r
 
92
        fac = -hl0/r
 
93
        call f (neq, tn, y, ftem)
 
94
      if(iero.gt.0) return
 
95
        do 220 i = 1,n
 
96
 220      wm(i+j1) = (ftem(i) - savf(i))*fac
 
97
        y(j) = yj
 
98
        j1 = j1 + n
 
99
 230    continue
 
100
      nfe = nfe + n
 
101
 240  continue
 
102
c compute norm of jacobian. --------------------------------------------
 
103
      pdnorm = fnorm (n, wm(3), ewt)/abs(hl0)
 
104
c add identity matrix. -------------------------------------------------
 
105
      j = 3
 
106
      do 250 i = 1,n
 
107
        wm(j) = wm(j) + 1.0d+0
 
108
 250    j = j + (n + 1)
 
109
c do lu decomposition on p. --------------------------------------------
 
110
      call dgefa (wm(3), n, n, iwm(21), ier)
 
111
      if (ier .ne. 0) ierpj = 1
 
112
      return
 
113
c dummy block only, since miter is never 3 in this routine. ------------
 
114
 300  return
 
115
c if miter = 4, call jac and multiply by scalar. -----------------------
 
116
 400  ml = iwm(1)
 
117
      mu = iwm(2)
 
118
      ml3 = ml + 3
 
119
      mband = ml + mu + 1
 
120
      meband = mband + ml
 
121
      lenp = meband*n
 
122
      do 410 i = 1,lenp
 
123
 410    wm(i+2) = 0.0d+0
 
124
      call jac (neq, tn, y, ml, mu, wm(ml3), meband)
 
125
      if(iero.gt.0) return
 
126
      con = -hl0
 
127
      do 420 i = 1,lenp
 
128
 420    wm(i+2) = wm(i+2)*con
 
129
      go to 570
 
130
c if miter = 5, make mband calls to f to approximate j. ----------------
 
131
 500  ml = iwm(1)
 
132
      mu = iwm(2)
 
133
      mband = ml + mu + 1
 
134
      mba = min(mband,n)
 
135
      meband = mband + ml
 
136
      meb1 = meband - 1
 
137
      srur = wm(1)
 
138
      fac = vmnorm (n, savf, ewt)
 
139
      r0 = 1000.0d+0*abs(h)*uround*dble(n)*fac
 
140
      if (r0 .eq. 0.0d+0) r0 = 1.0d+0
 
141
      do 560 j = 1,mba
 
142
        do 530 i = j,n,mband
 
143
          yi = y(i)
 
144
          r = max(srur*abs(yi),r0/ewt(i))
 
145
 530      y(i) = y(i) + r
 
146
        call f (neq, tn, y, ftem)
 
147
      if(iero.gt.0) return
 
148
        do 550 jj = j,n,mband
 
149
          y(jj) = yh(jj,1)
 
150
          yjj = y(jj)
 
151
          r = max(srur*abs(yjj),r0/ewt(jj))
 
152
          fac = -hl0/r
 
153
          i1 = max(jj-mu,1)
 
154
          i2 = min(jj+ml,n)
 
155
          ii = jj*meb1 - ml + 2
 
156
          do 540 i = i1,i2
 
157
 540        wm(ii+i) = (ftem(i) - savf(i))*fac
 
158
 550      continue
 
159
 560    continue
 
160
      nfe = nfe + mba
 
161
 570  continue
 
162
c compute norm of jacobian. --------------------------------------------
 
163
      pdnorm = bnorm (n, wm(3), meband, ml, mu, ewt)/abs(hl0)
 
164
c add identity matrix. -------------------------------------------------
 
165
      ii = mband + 2
 
166
      do 580 i = 1,n
 
167
        wm(ii) = wm(ii) + 1.0d+0
 
168
 580    ii = ii + meband
 
169
c do lu decomposition of p. --------------------------------------------
 
170
      call dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier)
 
171
      if (ier .ne. 0) ierpj = 1
 
172
      return
 
173
c----------------------- end of subroutine prja ------------------------
 
174
      end