1 subroutine ddasj (neq, ldd, x, y, yprime, delta, h, wt, e, 2 & wm, iwork, rwork, ddasf, info, ires) 3c Copyright (c) 2006, Math a la Carte, Inc. 4c>> 2010-08-26 ddasj Krogh Changed declaration of info to info(*). 5c>> 2008-08-26 ddasj Hanson add argument of leading dimension to ddasf 6c>> 2001-12-12 ddasj Krogh Changed code for reverse communication 7c>> 2001-11-23 ddasj Krogh Changed many names per library conventions. 8c>> 2001-11-04 ddasj Krogh Fixes for F77 and conversion to single & C 9c>> 2001-11-01 ddasj Hanson Provide code to Math a la Carte. 10c--D replaces "?": ?dasj, ?dasf, ?daslx, ?dasdb, ?swap, ?gbfa, ?gefa 11c***BEGIN PROLOGUE DDASJ 12c***SUBSIDIARY 13c***PURPOSE Compute the iteration matrix for DDASLX and form the 14c LU-decomposition. 15c***LIBRARY SLATEC (DDASLX) 16c***AUTHOR Petzold, Linda R., (LLNL) 17c***DESCRIPTION 18c----------------------------------------------------------------------- 19 20c THIS ROUTINE COMPUTES THE ITERATION MATRIX PD=DG/DY+CJ*DG/DYPRIME 21c (WHERE G(X,Y,YPRIME)=0, AND CJ IS CONTAINED IN RWORK(LCJ)). 22c HERE PD IS COMPUTED BY THE USER-SUPPLIED ROUTINE DDASF IF |INFO(5)| 23c IS 2, 4, 9, OR 13 AND IS COMPUTED BY FINITE DIFFERENCES IF |INFO(5)| 24c IS 1, 3, 8, OR 12. IF INFO(5) < 0, THEN COMPUTATIONS OF J ARE DONE 25c BY USING REVERSE COMMUNICTION IF THESE COMPUTATIONS ARE DONE BY THE 26c THE USER. WHEN |INFO(5) = 5 OR 6, THE USER IS DOING ALL 27c COMPUTATIONS ASSOCIATED WITH J AND THE ASSOCIATED LINEAR ALGEBRA. 28c THE PARAMETERS HAVE THE FOLLOWING MEANINGS. 29c NEQ = NUMBER OF EQUATIONS 30c LDD = FIRST (ROW) DIMENSION OF THE MATRIX. 31c X = CURRENT VALUE OF THE INDEPENDENT VARIABLE. 32c Y = ARRAY CONTAINING PREDICTED VALUES 33c YPRIME = ARRAY CONTAINING PREDICTED DERIVATIVES 34c DELTA = RESIDUAL EVALUATED AT (X,Y,YPRIME) 35c (USED ONLY FOR BAND MATRICES AND FOR REVERSE 36c COMMUNICATION.) 37c H = CURRENT STEPSIZE IN INTEGRATION 38c WT = VECTOR OF WEIGHTS FOR COMPUTING NORMS 39c E = WORK SPACE (TEMPORARY) OF LENGTH NEQ 40c WM = REAL WORK SPACE FOR MATRICES. ON 41c OUTPUT IT CONTAINS THE LU DECOMPOSITION 42c OF THE ITERATION MATRIX. 43c IWORK = INTEGER WORK SPACE CONTAINING MATRIX INFORMATION 44c DDASF = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE 45c TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME) 46c IRES = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES 47c IN DDASF, AND LESS THAN ZERO OTHERWISE. (IF IRES 48c IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED) 49c----------------------------------------------------------------------- 50c***ROUTINES CALLED DGBFA, DGEFA 51c***REVISION HISTORY (YYMMDD) 52c 830315 DATE WRITTEN 53c 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch) 54c 901010 Modified three MAX calls to be all on one line. (FNF) 55c 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format. 56c 901026 Added explicit declarations for all variables and minor 57c cosmetic changes to prologue. (FNF) 58c 901101 Corrected PURPOSE. (FNF) 59c***END PROLOGUE DDASJ 60c 61 integer ldd, neq, iwork(*), ires, info(*) 62 double precision x, y(*), yprime(*), delta(*), h, wt(*), 63 & e(*), wm(*), rwork(*) 64 external ddasf 65c 66 external dgbfa, dgefa 67c 68 integer i, i1, i2, ii, lmata, ipsave, isave, j, k, l, 69 & mba, mband, meb1, meband, msave, 70 & n, nrow 71 integer ier, locate 72 double precision del, delinv, squr, ypsave, ysave 73c 74c 75c POINTERS INTO IWORK 76 integer lml, lmu, lires, ldelt, lwm, lmxord, lk, lkold, lmat, 77 & lcnstr, lns, lnstl, lnst, lnre, lnje, letf, lctf, lnpd, 78 & ljcalc, lphase, revloc, mxstep, le, lwt, lphi, ntemp, lipvt 79 parameter (lml=1, lmu=2, lires=3, ldelt=4, lwm=5, lmxord=6, lk=7, 80 & lkold=8, lmat=9, lcnstr=10, lns=11, lnstl=12, lnst=13, 81 & lnre=14, lnje=15, letf=16, lctf=17, lnpd=18, ljcalc=19, 82 & lphase=20, revloc=21, mxstep=22, le=23, lwt=24, lphi=25, 83 & ntemp=26, lipvt=31) 84c 85c POINTERS INTO RWORK 86 integer lcj, ltstop, lhmax, lh, ltn, lcjold, lhold, lnjac, 87 & lround, lhmin, lalpha, lbeta, lgamma, lpsi, lsigma, ldelta 88 parameter (lcj=1, ltstop=2, lhmax=3, lh=4, ltn=5, lcjold=6, 89 & lhold=7, lnjac=8, lround=9, lhmin=10, lalpha=11, lbeta=17, 90 & lgamma=23, lpsi=29, lsigma=35, ldelta=46) 91c 92c POINTERS INTO INFO 93 integer itol, iout, istop, imat, idb, imaxh, ih0, iord, icnstr, 94 & inityp, ixstep 95 parameter (itol=2, iout=3, istop=4, imat=5, idb=6, imaxh=7, 96 & ih0=8, iord=9, icnstr=10, inityp=11, ixstep= 12) 97 98 save 99c***FIRST EXECUTABLE STATEMENT DDASJ 100 if (iwork(revloc) .ne. 0) then 101 ires = iwork(lires) 102 locate = mod(iwork(revloc), 8) 103 iwork(revloc) = iwork(revloc) / 8 104 go to (90, 50, 210, 170), locate 105 end if 106 107c The first entry drops through here. 108 lmata = abs(iwork(lmat)) 109c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 110 go to (30,10,140,120,250,280,30,10,140,120,30,10,140,120), lmata 111 go to 30 112c 113c 114c Dense full user-supplied matrix 115 10 CONTINUE 116 do 20 i=1, iwork(lnpd) 117 wm(i) = 0.0d0 118 20 continue 119 ires = 2 120 if (info(idb) .ne. 0) call ddasdb(2, neq, x, y, yprime, 121 & info, rwork, iwork, ires, rwork, rwork) 122 if (iwork(lmat) .ge. 0) then 123 call ddasf(x, y, yprime, e, wm, LDD, rwork(lcj), ires, 124 $ rwork, iwork) 125 go to 90 126 end if 127c Reverse communication to compute the partials. 128 iwork(revloc) = 1 129 iwork(lires) = ires 130c This location will be the same as WM when the user responds. 131 go to 240 132c 133c 134c DENSE FINITE-DIFFERENCE-GENERATED MATRIX 135 30 continue 136 ires = 1 137 nrow = 0 138 squr = sqrt(rwork(lround)) 139 i = 0 140c Loop over the columns to generate the approximate Jacobian. 141 40 i = i + 1 142 if (i .gt. neq) go to 80 143 del = squr*max(abs(y(i)),max(abs(h*yprime(i)),abs(wt(i)))) 144 del = sign(del,h*yprime(i)) 145 del = (y(i)+del) - y(i) 146 ysave = y(i) 147 ypsave = yprime(i) 148 y(i) = y(i) + del 149 yprime(i) = yprime(i) + rwork(lcj)*del 150 if (info(idb) .ne. 0) call ddasdb(2, neq, x, y, yprime, 151 & info, rwork, iwork, ires, rwork, rwork) 152 if (iwork(lmat) .ge. 0) then 153 call ddasf(x, y, yprime, e, wm, LDD, rwork(lcj), ires, 154 & rwork, iwork) 155 go to 60 156 else 157 iwork(revloc) = 2 158 iwork(lires) = ires 159c This is placed in the start of DELTA(*), in units of RWORK(*). 160 call dswap (neq, e, 1, delta, 1) 161 go to 240 162 end if 163c REVERSE ENTRY 2: 164 50 continue 165 call dswap (neq, e, 1, delta, 1) 166 60 continue 167 if (info(idb) .ne. 0) call ddasdb(3, neq, x, y, yprime, 168 & info, rwork, iwork, ires, rwork, rwork) 169 if (ires .lt. 0) then 170 if (ires .ge. -2) go to 240 171 info(idb) = -ires 172 ires = 1 173 end if 174 delinv = 1.0d0/del 175 do 70 l=1, neq 176 70 wm(nrow+l) = (e(l)-delta(l))*delinv 177 nrow = nrow + neq 178 y(i) = ysave 179 yprime(i) = ypsave 180 go to 40 181 80 continue 182c 183c 184c DO DENSE-MATRIX LU DECOMPOSITION ON PD 185c REVERSE ENTRY 1: 186 90 continue 187 if (info(idb) .ne. 0) call ddasdb(3, neq, x, y, yprime, 188 & info, rwork, iwork, ires, rwork, rwork) 189 if (ires .lt. 0) then 190 if (ires .ge. -2) go to 240 191 info(idb) = -ires 192 ires = 0 193 end if 194 if (lmata .gt. 4) then 195 ires = 3 196 if (lmata .le. 10) go to 260 197 go to 290 198 else 199 call dgefa (wm, ldd, neq, iwork(lipvt), ier) 200 if (ier .eq. 0) ires = 0 201 end if 202 go to 240 203c 204c BANDED USER-SUPPLIED MATRIX 205 120 do 130 i=1, iwork(lnpd) 206 wm(i) = 0.0d0 207 130 continue 208 ires = 2 209 if (info(idb) .ne. 0) call ddasdb(2, neq, x, y, yprime, 210 & info, rwork, iwork, ires, rwork, rwork) 211 if (iwork(lmat) .lt. 0) then 212 iwork(revloc) = 3 213 iwork(lires) = ires 214 go to 240 215 end if 216 call ddasf (x, y, yprime, e, wm, LDD, rwork(lcj), ires, 217 & rwork, iwork) 218 meband = 2*iwork(lml) + iwork(lmu) + 1 219 go to 210 220 221c 222c 223c BANDED FINITE-DIFFERENCE-GENERATED MATRIX 224 140 mband = iwork(lml) + iwork(lmu) + 1 225 mba = min(mband,neq) 226 meband = mband + iwork(lml) 227 meb1 = meband - 1 228 msave = (neq/mband) + 1 229 isave = iwork(ntemp) - 1 230 ipsave = isave + msave 231 ires = 1 232 squr = sqrt(rwork(lround)) 233 j = 0 234 150 continue 235 j = j + 1 236 if (j .gt. mba) go to 220 237 do 160 n=j, neq, mband 238 k = (n-j)/mband + 1 239 wm(isave+k) = y(n) 240 wm(ipsave+k) = yprime(n) 241 del = squr*max(abs(y(n)),max(abs(h*yprime(n)),abs(wt(n)))) 242 del = sign(del,h*yprime(n)) 243 del = (y(n)+del) - y(n) 244 y(n) = y(n) + del 245 160 yprime(n) = yprime(n) + rwork(lcj)*del 246 if (info(idb) .ne. 0) call ddasdb(2, neq, x, y, yprime, 247 & info, rwork, iwork, ires, rwork, rwork) 248 if (iwork(lmat) .ge. 0) then 249 call ddasf(x, y, yprime, e, wm, LDD, rwork(lcj), ires, 250 & rwork, iwork) 251 if (ires .lt. 0) then 252 if (ires .ge. -2) go to 240 253 info(idb) = -ires 254 ires = 0 255 end if 256 go to 180 257 else 258 iwork(revloc) = 4 259 ires = 1 260 iwork(lires) = ires 261c This is placed in the start of DELTA(*), in units of RWORK(*). 262 call dswap (neq, e, 1, delta, 1) 263 go to 240 264 end if 265c REVERSE ENTRY 4: 266 170 continue 267 call dswap (neq, e, 1, delta, 1) 268 180 continue 269 if (info(idb) .ne. 0) call ddasdb(3, neq, x, y, yprime, 270 & info, rwork, iwork, ires, rwork, rwork) 271 do 200 n=j, neq, mband 272 k = (n-j)/mband + 1 273 y(n) = wm(isave+k) 274 yprime(n) = wm(ipsave+k) 275 del = squr*max(abs(y(n)),max(abs(h*yprime(n)),abs(wt(n)))) 276 del = sign(del,h*yprime(n)) 277 del = (y(n)+del) - y(n) 278 delinv = 1.0d0/del 279 i1 = max(1,(n-iwork(lmu))) 280 i2 = min(neq,(n+iwork(lml))) 281 ii = n*meb1 - iwork(lml) 282 do 190 i=i1, i2 283 190 wm(ii+i) = (e(i)-delta(i))*delinv 284 200 continue 285 go to 150 286c 287c DO LU DECOMPOSITION OF BANDED PD 288c REVERSE ENTRY 3 289 210 if (info(idb) .ne. 0) call ddasdb(3, neq, x, y, yprime, 290 & info, rwork, iwork, ires, rwork, rwork) 291 220 continue 292 if (ires .lt. 0) then 293 if (ires .ge. -2) go to 240 294 info(idb) = -ires 295 ires = 0 296 end if 297 if (lmata .gt. 4) then 298 ires = 3 299 if (lmata .le. 10) go to 260 300 go to 290 301 else 302 call dgbfa (wm, meband, neq, iwork(lml), iwork(lmu), 303 & iwork(lipvt), ier) 304 if (ier .eq. 0) ires = 0 305 end if 306 240 continue 307 return 308c 309c User defined matrix, Get Jacobian and factor in one step 310 250 continue 311 ires = 2 312c Enters here if already have matrix and want to factor 313 260 continue 314 if (info(idb) .ne. 0) call ddasdb(2, neq, x, y, yprime, 315 & info, rwork, iwork, ires, rwork, rwork) 316 call ddasf (x, y, yprime, e, wm, LDD, rwork(lcj), ires, 317 & rwork, iwork) 318 319 if (info(idb) .ne. 0) call ddasdb(3, neq, x, y, yprime, 320 & info, rwork, iwork, ires, rwork, rwork) 321 if (ires .lt. -2) then 322 info(idb) = -ires 323 ires = 0 324 end if 325 go to 240 326c User defined but with reverse communication 327 280 ires = 2 328c Enter here if matrix is computed. 329 290 continue 330 iwork(lires) = ires 331 iwork(revloc) = -1 332 go to 240 333c------ END OF SUBROUTINE DDASJ ------ 334 end 335