/[ascend]/trunk/sdassl/isrdassl.for
ViewVC logotype

Contents of /trunk/sdassl/isrdassl.for

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (21 years ago) by aw0a
File size: 220225 byte(s) 
Setting up web subdirectory in repository
1 DOUBLE PRECISION FUNCTION D1MACH (IDUM)
2 INTEGER IDUM
3 C
4 C***********************************************************************
5 C
6 C SDASSL routines
7 C by Linda Petzold, Andreas Kroener, Wolfgang Marquardt
8 C Created: 15/03/83
9 C Version: 1.1 Rev: 1989/12/11
10 C Date last modified: 1994/09/02
11 C
12 C This file is part of the SDASSL differential/algebraic system solver.
13 C
14 C Copyright (C) 1983, 1989, 1994 Linda Petzold, Andreas Kroener,
15 C Wolfgang Marquardt
16 C
17 C The SDASSL differential/algebraic system solver is free software;
18 C you can redistribute it and/or modify it under the terms of the GNU
19 C General Public License as published by the Free Software Foundation;
20 C either version 2 of the License, or (at your option) any later version.
21 C
22 C The SDASSL system solver is distributed in hope that it will be
23 C useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
24 C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25 C General Public License for more details.
26 C
27 C You should have received a copy of the GNU General Public License along
28 C with the program; if not, write to the Free Software Foundation, Inc.,
29 C 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
30 C
31 C***********************************************************************
32 C
33 C-----------------------------------------------------------------------
34 C THIS ROUTINE COMPUTES THE UNIT ROUNDOFF OF THE MACHINE IN DOUBLE
35 C PRECISION. THIS IS DEFINED AS THE SMALLEST POSITIVE MACHINE NUMBER
36 C U SUCH THAT 1.0D0 + U .NE. 1.0D0 (IN DOUBLE PRECISION).
37 C-----------------------------------------------------------------------
38 DOUBLE PRECISION U, COMP
39 U = 1.0D0
40 10 U = U*0.5D0
41 COMP = 1.0D0 + U
42 IF (COMP .NE. 1.0D0) GO TO 10
43 D1MACH = U*2.0D0
44 RETURN
45 C----------------------- END OF FUNCTION D1MACH ------------------------
46 END
47 SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
48 C
49 C CONSTANT TIMES A VECTOR PLUS A VECTOR.
50 C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
51 C JACK DONGARRA, LINPACK, 3/11/78.
52 C
53 DOUBLE PRECISION DX(1),DY(1),DA
54 INTEGER I,INCX,INCY,IX,IY,M,MP1,N
55 C
56 IF(N.LE.0)RETURN
57 IF (DA .EQ. 0.0D0) RETURN
58 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
59 C
60 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
61 C NOT EQUAL TO 1
62 C
63 IX = 1
64 IY = 1
65 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
66 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
67 DO 10 I = 1,N
68 DY(IY) = DY(IY) + DA*DX(IX)
69 IX = IX + INCX
70 IY = IY + INCY
71 10 CONTINUE
72 RETURN
73 C
74 C CODE FOR BOTH INCREMENTS EQUAL TO 1
75 C
76 C
77 C CLEAN-UP LOOP
78 C
79 20 M = MOD(N,4)
80 IF( M .EQ. 0 ) GO TO 40
81 DO 30 I = 1,M
82 DY(I) = DY(I) + DA*DX(I)
83 30 CONTINUE
84 IF( N .LT. 4 ) RETURN
85 40 MP1 = M + 1
86 DO 50 I = MP1,N,4
87 DY(I) = DY(I) + DA*DX(I)
88 DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
89 DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
90 DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
91 50 CONTINUE
92 RETURN
93 END
94 SUBROUTINE DDAINI(X,Y,YPRIME,NEQ,
95 * RES,JAC,H,WT,IDID,RPAR,IPAR,
96 * PHI,DELTA,E,WM,IWM,
97 * HMIN,UROUND,NONNEG,PTN)
98 C
99 C***BEGIN PROLOGUE DDAINI
100 C***REFER TO DDASSL
101 C***ROUTINES CALLED DDANRM,DDAJAC,DDASLV
102 C***COMMON BLOCKS DDA001
103 C***DATE WRITTEN 830315 (YYMMDD)
104 C***REVISION DATE 830315 (YYMMDD)
105 C***END PROLOGUE DDAINI
106 C
107 C-------------------------------------------------------
108 C ddaini takes one step of size h or smaller
109 C with the backward euler method, to
110 C find yprime at the initial time x. a modified
111 C damped newton iteration is used to
112 C solve the corrector iteration.
113 C
114 C the initial guess yprime is used in the
115 C prediction, and in forming the iteration
116 C matrix, but is not involved in the
117 C error test. this may have trouble
118 C converging if the initial guess is no
119 C good, or if g(xy,yprime) depends
120 C nonlinearly on yprime.
121 C
122 C the parameters represent:
123 C x -- independent variable
124 C y -- solution vector at x
125 C yprime -- derivative of solution vector
126 C neq -- number of equations
127 C h -- stepsize. imder may use a stepsize
128 C smaller than h.
129 C wt -- vector of weights for error
130 C criterion
131 C idid -- completion code with the following meanings
132 C idid= 1 -- yprime was found successfully
133 C idid=-12 -- ddaini failed to find yprime
134 C rpar,ipar -- real and integer parameter arrays
135 C that are not altered by ddaini
136 C phi -- work space for ddaini
137 C delta,e -- work space for ddaini
138 C wm,iwm -- real and integer arrays storing
139 C matrix information
140 C
141 C-----------------------------------------------------------------
142 C
143 C C_1 R. KOE EINBAU VON PTN (EXT, PARAMETERKLAMMERN DDAINI,DDAJAC)
144 C
145 C
146 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
147 LOGICAL CONVGD
148 DIMENSION Y(1),YPRIME(1),WT(1)
149 DIMENSION PHI(NEQ,1),DELTA(1),E(1)
150 DIMENSION WM(1),IWM(1)
151 DIMENSION RPAR(1),IPAR(1)
152 EXTERNAL RES,JAC,PTN
153 COMMON/DDA001/NPD,NTEMP,
154 * LML,LMU,LMXORD,LMTYPE,
155 * LNST,LNRE,LNJE,LETF,LCTF,LIPVT
156
157 DATA MAXIT/10/,MJAC/5/
158 DATA DAMP/0.75D0/
159
160 C
161 C
162 C---------------------------------------------------
163 C block 1.
164 C initializations.
165 C---------------------------------------------------
166 C
167 IDID=1
168 NEF=0
169 NCF=0
170 NSF=0
171 YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR)
172 C
173 C save y and yprime in phi
174 DO 100 I=1,NEQ
175 PHI(I,1)=Y(I)
176 100 PHI(I,2)=YPRIME(I)
177
178 C
179 C
180 C----------------------------------------------------
181 C block 2.
182 C do one backward euler step.
183 C----------------------------------------------------
184 C
185 C set up for start of corrector iteration
186 200 CJ=1.0D0/H
187 XNEW=X+H
188 C
189 C predict solution and derivative
190
191 DO 250 I=1,NEQ
192 250 Y(I)=Y(I)+H*YPRIME(I)
193 C
194 JCALC=-1
195 M=0
196 CONVGD=.TRUE.
197 C
198 C
199 C corrector loop.
200 300 IWM(LNRE)=IWM(LNRE)+1
201 IRES=0
202
203 CALL RES(XNEW,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
204 IF (IRES.LT.0) GO TO 430
205 C
206 C
207 C evaluate the iteration matrix
208 IF (JCALC.NE.-1) GO TO 310
209 IWM(LNJE)=IWM(LNJE)+1
210 JCALC=0
211 CALL DDAJAC(NEQ,XNEW,Y,YPRIME,DELTA,CJ,H,
212 * IER,WT,E,WM,IWM,RES,IRES,
213 * UROUND,JAC,RPAR,IPAR,PTN)
214
215 S=1000000.D0
216 IF (IRES.LT.0) GO TO 430
217 IF (IER.NE.0) GO TO 430
218 NSF=0
219
220 C
221 C
222 C
223 C multiply residual by damping factor
224 310 CONTINUE
225 DO 320 I=1,NEQ
226 320 DELTA(I)=DELTA(I)*DAMP
227
228 C
229 C compute a new iterate (back substitution)
230 C store the correction in delta
231
232 CALL DDASLV(NEQ,DELTA,WM,IWM)
233
234 C
235 C update y and yprime
236
237 DO 330 I=1,NEQ
238 Y(I)=Y(I)-DELTA(I)
239 330 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
240
241 C
242 C test for convergence of the iteration.
243
244 DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
245 IF (DELNRM.LE.100.D0*UROUND*YNORM)
246 * GO TO 400
247
248 IF (M.GT.0) GO TO 340
249 OLDNRM=DELNRM
250 GO TO 350
251
252 340 RATE=(DELNRM/OLDNRM)**(1.0D0/DFLOAT(M))
253 IF (RATE.GT.0.90D0) GO TO 430
254 S=RATE/(1.0D0-RATE)
255
256 350 IF (S*DELNRM .LE. 0.33D0) GO TO 400
257 C
258 C
259 C the corrector has not yet converged. update
260 C m and and test whether the maximum
261 C number of iterations have been tried.
262 C every mjac iterations, get a new
263 C iteration matrix.
264
265 M=M+1
266 IF (M.GE.MAXIT) GO TO 430
267
268 IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1
269
270 GO TO 300
271
272 C
273 C
274 C the iteration has converged.
275 C check nonnegativity constraints
276 400 IF (NONNEG.EQ.0) GO TO 450
277 DO 410 I=1,NEQ
278 410 DELTA(I)=DMIN1(Y(I),0.0D0)
279
280 DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
281 IF (DELNRM.GT.0.33D0) GO TO 430
282
283 DO 420 I=1,NEQ
284 Y(I)=Y(I)-DELTA(I)
285 420 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
286 GO TO 450
287 C
288 C
289 C exits from corrector loop.
290 430 CONVGD=.FALSE.
291 450 IF (.NOT.CONVGD) GO TO 600
292 C
293 C
294 C
295 C-----------------------------------------------------
296 C block 3.
297 C the corrector iteration converged.
298 C do error test.
299 C-----------------------------------------------------
300 C
301
302 DO 510 I=1,NEQ
303 510 E(I)=Y(I)-PHI(I,1)
304
305 ERR=DDANRM(NEQ,E,WT,RPAR,IPAR)
306
307 IF (ERR.LE.1.0D0) RETURN
308
309 C
310 C
311 C
312 C--------------------------------------------------------
313 C block 4.
314 C the backward euler step failed. restore y
315 C and yprime to their original values.
316 C reduce stepsize and try again, if
317 C possible.
318 C---------------------------------------------------------
319 C
320
321 600 CONTINUE
322 DO 610 I=1,NEQ
323 Y(I)=PHI(I,1)
324 610 YPRIME(I)=PHI(I,2)
325
326 IF (CONVGD) GO TO 640
327 IF (IER.EQ.0) GO TO 620
328 NSF=NSF+1
329 H=H*0.25D0
330 IF (NSF.LT.3.AND.DABS(H).GE.HMIN) GO TO 690
331 IDID=-12
332 RETURN
333 620 IF (IRES.GT.-2) GO TO 630
334 IDID=-12
335 RETURN
336 630 NCF=NCF+1
337 H=H*0.25D0
338 IF (NCF.LT.10.AND.DABS(H).GE.HMIN) GO TO 690
339 IDID=-12
340 RETURN
341
342 640 NEF=NEF+1
343 R=0.90D0/(2.0D0*ERR+0.0001D0)
344 R=DMAX1(0.1D0,DMIN1(0.5D0,R))
345 H=H*R
346 IF (DABS(H).GE.HMIN.AND.NEF.LT.10) GO TO 690
347 IDID=-12
348 RETURN
349 690 GO TO 200
350
351 C-------------end of subroutine ddaini----------------------
352 END
353 SUBROUTINE DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
354 * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,IPAR,PTN)
355 C
356 C***BEGIN PROLOGUE DDAJAC
357 C***REFER TO DDASSL
358 C***ROUTINES CALLED DGEFA,DGBFA
359 C***COMMON BLOCKS DDA001
360 C***DATE WRITTEN 830315 (YYMMDD)
361 C***REVISION DATE 830315 (YYMMDD)
362 C***END PROLOGUE DDAJAC
363 C-----------------------------------------------------------------------
364 C this routine computes the iteration matrix
365 C pd=dg/dy+cj*dg/dyprime (where g(x,y,yprime)=0).
366 C here pd is computed by the user-supplied
367 C routine jac if iwm(mtype) is 1 or 4, and
368 C it is computed by numerical finite differencing
369 C if iwm(mtype)is 2 or 5
370 C the parameters have the following meanings.
371 C y = array containing predicted values
372 C yprime = array containing predicted derivatives
373 C delta = residual evaluated at (x,y,yprime)
374 C_2 (used only if iwm(mtype)=2 or 5)
375 C (used only if iwm(itype)=2 or 3)
376 C cj = scalar parameter defining iteration matrix
377 C h = current stepsize in integration
378 C ier = variable which is .ne. 0
379 C if iteration matrix is singular or could not be
380 C decomposed for another reason,
381 C and 0 otherwise.
382 C wt = vector of weights for computing norms
383 C e = work space (temporary) of length neq
384 C wm = real work space for matrices. on
385 C output it contains the lu decomposition
386 C of the iteration matrix.
387 C iwm = integer work space containing
388 C matrix information
389 C res = name of the external user-supplied routine
390 C to evaluate the residual function g(x,y,yprime)
391 C ires = flag which is equal to zero if no illegal values
392 C in res, and less than zero otherwise. (if ires
393 C is less than zero, the matrix was not completed)
394 C in this case (if ires .lt. 0), then ier = 0.
395 C uround = the unit roundoff error of the machine being used.
396 C jac = name of the external user-supplied routine
397 C to evaluate the iteration matrix (this routine
398 C is only used if iwm(mtype) is 1 or 4)
399 C-----------------------------------------------------------------------C
400 C
401 C C_2 R.KOENIGSDORFF 18.9.86 BER. D. NUM. ABLEITUNG FUER SPARSE MATRIX
402 C C_3 " " " " " " EINBAU VON PTN IN PARAMETERKLAMMER UND EXT
403 C C_5 " " " 26.9.86 ERWEITERUNG COMMON/DDA001/
404 C C_6 " " " 20.10.86 EINBAU VON DDALDJ UND MA30BD
405 C C_9 " " " 12.01.87 AENDERUNG DER SKELETT-MATRIX FUER PIVOT
406 C C_10 " " " 16.05.87 ZURUECK ZU MA30AD BEI FEHLER IN MA30BD
407 C
408 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
409 EXTERNAL RES,JAC,PTN
410 DIMENSION Y(1),YPRIME(1),DELTA(1),WT(1),E(1)
411 DIMENSION WM(1),IWM(1),RPAR(1),IPAR(1)
412 C_5
413 COMMON/DDA001/NPD,NTEMP,
414 * LML,LMU,LMXORD,LMTYPE,
415 * LNST,LNRE,LNJE,LETF,LCTF,LIPVT,
416 * KWKN,KJAC,KLUD,KIPTR,KLENB,KICNB,KLENR,KLENRL,
417 * KIPA,KIQA,KICN,JIRN,KLENC,JIFIRST,JLASTR,JNEXTR,JLASTC,JNEXTC,
418 * JIPC,IDISP(2),LENOFF(1),JACSIZ,LUDSIZ,PT1,NCG,KICNG,KICGP,LIRN
419 C_5 END
420 COMMON/MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN
421 C
422 IER = 0
423 C_10
424 IFMA30=0
425 C_10
426 NPDM1=NPD-1
427 MTYPE=IWM(LMTYPE)
428 GO TO (100,200,300),MTYPE
429 C
430 C
431 C dense user-supplied matrix
432 100 LENPD=NEQ*NEQ
433 DO 110 I=1,LENPD
434 110 WM(NPDM1+I)=0.0D0
435 CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
436 GO TO 230
437 C
438 C
439 C dense finite-difference-generated matrix
440 200 IRES=0
441 NROW=NPDM1
442 SQUR = DSQRT(UROUND)
443 DO 210 I=1,NEQ
444 DEL=SQUR*DMAX1(DABS(Y(I)),DABS(H*YPRIME(I)),
445 * DABS(WT(I)))
446 DEL=DSIGN(DEL,H*YPRIME(I))
447 DEL=(Y(I)+DEL)-Y(I)
448 YSAVE=Y(I)
449 YPSAVE=YPRIME(I)
450 Y(I)=Y(I)+DEL
451 YPRIME(I)=YPRIME(I)+CJ*DEL
452 CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
453 IF (IRES .LT. 0) RETURN
454 DELINV=1.0D0/DEL
455 DO 220 L=1,NEQ
456 220 WM(NROW+L)=(E(L)-DELTA(L))*DELINV
457 NROW=NROW+NEQ
458 Y(I)=YSAVE
459 YPRIME(I)=YPSAVE
460 210 CONTINUE
461 C
462 C
463 C do dense-matrix lu decomposition on pd
464 230 CALL DGEFA(WM(NPD),NEQ,NEQ,IWM(LIPVT),IER)
465 RETURN
466 C
467 C
468 C_2 dummy section for iwm(mtype)=3
469 C sparse finite-difference-generated matrix
470 300 IF (IWM(LNJE) .GT. 1) GO TO 310
471 CALL PTN(IWM(KLENB),NEQ,IWM(KICNB),JACSIZ,NIRN,RPAR,IPAR)
472 IF (JACSIZ.LT.NIRN) STOP 'DASSL: JACSIZ.LT.NIRN'
473 NEQ1=NEQ+1
474 NEQ2=2*NEQ
475 CALL DDASCO(IWM(KICNB),IWM(KLENB),NEQ,NEQ,IWM(KICNG),
476 * IWM(KICGP),NCG,IWM(KLENR),NEQ2,NEQ1,NIRN)
477 C
478 310 KLUDP=KLUD+NEQ
479 CALL DDADIF (NEQ,X,Y,YPRIME,DELTA,CJ,H,WT,E,WM(KJAC),WM(KLUD),
480 * WM(KWKN),WM(KLUDP),NCG,IWM(KIPTR),IWM(KLENB),IWM(KICNG),
481 * IWM(KICGP),IWM(KICNB),RES,IRES,UROUND,RPAR,IPAR)
482 IF (IRES .LT. 0) RETURN
483 IF (IWM(LNJE).GT.1) GO TO 340
484 350 CALL DDALDS(WM(KJAC),IWM(KLENB),IWM(KICNB),NIRN,WM(KLUD),
485 * IWM(KICN),LUDSIZ,IWM(KLENR),IDISP,IWM(KIPA),IWM(KIQA),
486 * NEQ, NEQ)
487 C_9 * NEQ, 0)
488 U=PT1
489 CALL MA30AD(NEQ,IWM(KICN),WM(KLUD),LUDSIZ,IWM(KLENR),
490 * IWM(KLENRL),IDISP,IWM(KIPA),IWM(KIQA),IWM(JIRN),LIRN,
491 * IWM(KLENC),IWM(JIFIRST),IWM(JLASTR),IWM(JNEXTR),IWM(JLASTC),
492 * IWM(JNEXTC),IWM(KIPTR),IWM(JIPC),U,IFLAG)
493 IF (IFLAG.NE.0) THEN
494 WRITE(*,*) ' DASSL, DDAJAC: MA30AD, IFLAG.NE.0, IFLAG =',
495 * IFLAG
496 IER = - 1
497 ENDIF
498 IF (IER .NE. 0) RETURN
499 C
500 C
501 C do sparse-matrix lu decomposition on pd
502 C_6
503 340 CALL DDALDJ(WM(KJAC),IWM(KICNB),NIRN,IWM(KLENB),WM(KLUD),
504 * IWM(KICN),LUDSIZ,IWM(KLENR),IWM(KIPA),IWM(KIQA),IDISP,
505 * WM(KWKN),IWM(KIPTR),NEQ)
506 CALL MA30BD(NEQ,IWM(KICN),WM(KLUD),LUDSIZ,IWM(KLENR),
507 * IWM(KLENRL),IDISP,IWM(KIPA),IWM(KIQA),WM(KWKN),
508 * IWM(KIPTR),IFLAG)
509 IF (IFLAG.NE.0) WRITE(*,*) ' DASSL, DDAJAC: MA30BD, IFLAG.NE.0'
510 IF (IFLAG.GT.0) WRITE(*,*) ' Pivot I is very small, I =',IFLAG
511 IF (IFLAG.LT.0) WRITE(*,*) ' Unexpected singularity at stage I
512 * of the decomposition, I =',IFLAG
513 C_6
514 C_10
515 IF (IFLAG.NE.0 .AND. IFMA30.EQ.0) THEN
516 IFMA30=1
517 GOTO 350
518 ELSE IF (IFLAG .NE. 0) THEN
519 IER = -2
520 ENDIF
521 C_10
522 RETURN
523 C------end of subroutine ddajac------
524 END
525 DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
526 C
527 C***BEGIN PROLOGUE DDANRM
528 C***REFER TO DDASSL
529 C***ROUTINES CALLED (NONE)
530 C***DATE WRITTEN 830315 (YYMMDD)
531 C***REVISION DATE 830315 (YYMMDD)
532 C***END PROLOGUE DDANRM
533 C-----------------------------------------------------------------------
534 C this function routine computes the weighted
535 C root-mean-square norm of the vector of length
536 C neq contained in the array v,with weights
537 C contained in the array wt of length neq.
538 C ddanrm=sqrt((1/neq)*sum(v(i)/wt(i))**2)
539 C-----------------------------------------------------------------------
540 C
541 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
542 DIMENSION V(NEQ),WT(NEQ)
543 DIMENSION RPAR(1),IPAR(1)
544 DDANRM = 0.0D0
545 VMAX = 0.0D0
546 DO 10 I = 1,NEQ
547 10 IF(DABS(V(I)/WT(I)) .GT. VMAX) VMAX = DABS(V(I)/WT(I))
548 IF(VMAX .LE. 0.0D0) GO TO 30
549 SUM = 0.0D0
550 DO 20 I = 1,NEQ
551 20 SUM = SUM + ((V(I)/WT(I))/VMAX)**2
552 DDANRM = VMAX*DSQRT(SUM/DFLOAT(NEQ))
553 30 CONTINUE
554 RETURN
555 C------end of function ddanrm------
556 END
557 SUBROUTINE DDASLV(NEQ,DELTA,WM,IWM)
558 C
559 C***BEGIN PROLOGUE DDASLV
560 C***REFER TO DDASSL
561 C***ROUTINES CALLED DGESL,DGBSL
562 C***COMMON BLOCKS DDA001
563 C***DATE WRITTEN 830315 (YYMMDD)
564 C***REVISION DATE 830315 (YYMMDD)
565 C***END PROLOGUE DDASLV
566 C-----------------------------------------------------------------------
567 C this routine manages the solution of the linear
568 C system arising in the newton iteration.
569 C matrices and real temporary storage and
570 C real information are stored in the array wm.
571 C integer matrix information is stored in
572 C the array iwm.
573 C for a dense matrix, the linpack routine
574 C dgesl is called.
575 C for a banded matrix,the linpack routine
576 C dgbsl is called
577 C-----------------------------------------------------------------------
578 C
579 C C_4 R.KOENIGSDORFF 24.9.86 MA30CD EINGEBAUT
580 C C_5 " " " " 26.9.86 ERWEITERUNG COMMON/DDA001/
581 C
582 C
583 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
584 DIMENSION DELTA(1),WM(1),IWM(1)
585 C_5
586 COMMON/DDA001/NPD,NTEMP,LML,LMU,
587 * LMXORD,LMTYPE,
588 * LNST,LNRE,LNJE,LETF,LCTF,LIPVT,
589 * KWKN,KJAC,KLUD,KIPTR,KLENB,KICNB,KLENR,KLENRL,
590 * KIPA,KIQA,KICN,JIRN,KLENC,JIFIRST,JLASTR,JNEXTR,JLASTC,JNEXTC,
591 * JIPC,IDISP(2),LENOFF(1),JACSIZ,LUDSIZ,PT1,NCG,KICNG,KICGP,LIRN
592 C_5 END
593 C
594 MTYPE=IWM(LMTYPE)
595 GO TO(100,100,300),MTYPE
596 C
597 C dense matrix
598 100 CALL DGESL(WM(NPD),NEQ,NEQ,IWM(LIPVT),DELTA,0)
599 RETURN
600 C
601 C_4 sparse matrix
602 300 CALL MA30CD (NEQ,IWM(KICN),WM(KLUD),LUDSIZ,IWM(KLENR),IWM(KLENRL),
603 * LENOFF,IDISP,IWM(KIPA),IWM(KIQA),DELTA,WM(KWKN),2)
604 RETURN
605 C------end of subroutine ddaslv------
606 END
607 SUBROUTINE DDALDJ(B,ICNB,JACLEN,LENB,A,ICN,LUDLEN,LENR
608 1,IP,IQ,IDISP,W,IW,N)
609 DOUBLE PRECISION B(JACLEN), A(LUDLEN), W(N), ZERO
610 INTEGER DBLK, IW(N), IDISP(2)
611 INTEGER ICNB(JACLEN), LENB(N), ICN(LUDLEN), LENR(N)
612 1, IP(N),IQ(N)
613 DATA ZERO /0.D0/
614 C
615 C RELOAD TRANSPOSED NEWTON MATRIX, PARTITIONED (NALGB,NEQ-NALGB)
616 C
617 C T
618 C ( J(1,1) J(1,2) )
619 C ( J(2,1) J(2,2) )
620 C
621 C FOR DECOMPOSITION BY MA30B USING OLD PIVOTAL SEQUENCE
622 C
623 C IDISP IS THE POSITION IN A/ICN OF THE FIRST ELEMENT
624 LTF=1
625 DBLK=IDISP(1)
626 IW(1)=1
627 DO 10 I=1,N
628 IF(I.LT.N)IW(I+1)=IW(I)+LENB(I)
629 10 W(I)=ZERO
630 C EACH PASS THROUGH THIS MAIN LOOP PUTS ROW I OF THE PERMUTED FORM
631 C (ROW IOLD IN THE ORIGINAL MATRIX) INTO THE RIGHT PLACE IN A
632 DO 60 I=1,N
633 C LOAD ROW IOLD OF B TRANSPOSED INTO VECTOR W.
634 IOLD=IABS(IP(I)+0)
635 J1=IW(IOLD)
636 J2=J1+LENB(IOLD)-1
637 IF(J1.GT.J2)GO TO 30
638 DO 20 JJ=J1,J2
639 J=ICNB(JJ)
640 W(J)=B(JJ)
641 20 CONTINUE
642 30 IF (LENR(I).EQ.0) GO TO 60
643 C UNLOAD ROW IOLD (ROW I IN PERMUTED FORM)
644 C FROM W INTO APPROPRIATE PART OF A.
645 J1=DBLK
646 J2=J1+LENR(I)-1
647 DO 50 JJ=J1,J2
648 K=ICN(JJ)
649 J=IQ(K)
650 A(JJ)=W(J)
651 50 W(J)=ZERO
652 DBLK=J2+1
653 60 CONTINUE
654 RETURN
655 END
656 SUBROUTINE DDALDS(B,LENB,ICNB,JACLEN,A,ICN,LUDLEN,LENR,IDX
657 1,IP,IQ,N,IMPLI)
658 DOUBLE PRECISION B(JACLEN), A(LUDLEN), ZERO, ONE
659 INTEGER LENB(N), ICNB(JACLEN), ICN(LUDLEN), LENR(N),
660 1 IP(N), IQ(N)
661 INTEGER IDX(2)
662 DATA ZERO, ONE /0.D0, 1.D0/
663 C
664 C LOADS SKELETON NEWTON MATRIX (PARTITIONED IF NALGB .GT. 0)
665 C
666 C T
667 C ( I 0 )
668 C ( 0 I )
669 C
670 C FOR DECOMPOSITION BY MA30A TO CHOOSE PIVOTAL SEQUENCE
671 C
672 K2=0
673 J2=0
674 DO 5 I=1,N
675 J1=J2+1
676 J2=J2+LENB(I)
677 LR=LENB(I)
678 IF(I.LE.IMPLI) GO TO 3
679 LR=LR+1
680 IF(J1.GT.J2)GO TO 3
681 DO 2 J=J1,J2
682 IF(ICNB(J).EQ.I)LR=LR-1
683 2 CONTINUE
684 3 LENR(I)=LR
685 IP(I)=I
686 IQ(I)=I
687 5 K2=K2+LR
688 K1=LUDLEN+1-K2
689 IDX(1)=1
690 IDX(2)=K1
691 J2=0
692 DO 200 I=1,N
693 JM=LENR(I)
694 J1=J2+1
695 J2=J2+LENB(I)
696 IF(I.GT.IMPLI)GO TO 8
697 IF(J1.GT.J2)GO TO 20
698 DO 6 J=J1,J2
699 IC=ICNB(J)
700 A(K1)=ZERO
701 IF(IC.LE.IMPLI)A(K1)=B(J)
702 ICN(K1)=IC
703 6 K1=K1+1
704 GO TO 20
705 8 IF(J1.GT.J2)GO TO 15
706 J0=0
707 DO 10 J=J1,J2
708 IC=ICNB(J)
709 IF(IC.LE.IMPLI)GO TO 63
710 IF(IC-I)63,62,61
711 61 IF(J0.GT.0)GO TO 63
712 J0=1
713 A(K1)=ONE
714 ICN(K1)=I
715 K1=K1+1
716 GO TO 63
717 62 A(K1)=ONE
718 J0=1
719 GO TO 64
720 63 A(K1)=ZERO
721 64 ICN(K1)=IC
722 K1=K1+1
723 10 CONTINUE
724 IF(J0.NE.0)GO TO 20
725 15 A(K1)=ONE
726 ICN(K1)=I
727 K1=K1+1
728 20 CONTINUE
729 200 CONTINUE
730 RETURN
731 END
732 SUBROUTINE DDASCO(IR,IP,M,N,IC,IPC,NC,IW,MPN,NP1,JACSIZ)
733 C
734 C GROUPS COLUMNS OF SPARSE JACOBIAN MATRIX FOR FINITE
735 C DIFFERENCE EVALUATION
736 C
737 C M IS NUMBER OF ROWS
738 C N IS NUMBER OF COLUMNS
739 C IR CONTAINS ROW NUMBERS OF NON-ZEROS.
740 C IP CONTAINS NUMBERS OF NONZEROS IN EACH COLUMN.
741 C IC IS SET TO COLUMN NUMBERS WITHIN GROUPS
742 C AND MUST HAVE N ENTRIES.
743 C IPC IS SET TO POINT TO FIRST ENTRY OF IC FOR EACH
744 C GROUP, AND MAY NEED N+1 ENTRIES.
745 C NC IS SET TO (NUMBER OF GROUPS+1).
746 Change
747 C INTEGER*2 IR(JACSIZ),IP(N),IC(N),IPC(NP1)
748 INTEGER IR(JACSIZ),IP(N),IC(N),IPC(NP1)
749 Change
750 C IW IS WORKSPACE, IT NEEDS (M+N) ENTRIES
751 INTEGER IW(MPN)
752 NC=1
753 ICC=1
754 NM=N+M
755 DO 1 J=1,NM
756 1 IW(J)=0
757 10 IPC(NC)=ICC
758 DO 2 J=1,M
759 2 IW(J)=0
760 KCOL=1
761 JST=1
762 50 JND=JST+IP(KCOL)-1
763 K=KCOL+M
764 IF(IW(K).NE.0)GO TO 20
765 IF(JND.LT.JST)GO TO 30
766 DO 3 J=JST,JND
767 K=IR(J)
768 IF(IW(K).NE.0)GO TO 40
769 3 CONTINUE
770 C ACCEPT COLUMN
771 DO 4 J=JST,JND
772 K=IR(J)
773 4 IW(K)=1
774 30 IC(ICC)=KCOL
775 ICC=ICC+1
776 K=KCOL+M
777 IW(K)=1
778 C REJECT COLUMN
779 40 CONTINUE
780 C COLUMN ALREADY USED
781 20 KCOL=KCOL+1
782 JST=JND+1
783 IF(KCOL.LE.N)GO TO 50
784 IF(ICC.EQ.IPC(NC))GO TO 60
785 NC=NC+1
786 GO TO 10
787 60 RETURN
788 END
789 SUBROUTINE DDADIF(NEQ,X,Y,YPRIME,DELTA,CJ,H,WT,E,B,DELY,YSAVE,
790 * YPSAVE,NCG,IPTR,LENB,ICNG,ICGP,ICNB,RES,IRES,UROUND,RPAR,IPAR)
791 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
792 EXTERNAL RES
793 DIMENSION Y(1),YPRIME(1),DELTA(1),WT(1),E(1),YSAVE(1),YPSAVE(1)
794 DIMENSION IPTR(1),LENB(1),ICNG(1),ICGP(1),ICNB(1)
795 DIMENSION B(1),DELY(1),RPAR(1),IPAR(1)
796 C
797 C R.KOENIGSDORFF 24.10.86
798 C
799 IRES=0
800 SQUR=DSQRT(UROUND)
801 IPTR(1)=1
802 N1=NEQ-1
803 DO 5 I=1,N1
804 IPTR(I+1)=IPTR(I)+LENB(I)
805 5 CONTINUE
806 DO 10 KK=2,NCG
807 JST=ICGP(KK-1)
808 JND=ICGP(KK)-1
809 DO 20 K=JST,JND
810 J=ICNG(K)
811 DELY(J)=SQUR*DMAX1(DABS(Y(J)),DABS(H*YPRIME(J)),DABS(WT(J)))
812 DELY(J)=DSIGN(DELY(J),H*YPRIME(J))
813 DELY(J)=(Y(J)+DELY(J))-Y(J)
814 YSAVE(J)=Y(J)
815 YPSAVE(J)=YPRIME(J)
816 Y(J)=Y(J)+DELY(J)
817 YPRIME(J)=YPRIME(J)+CJ*DELY(J)
818 20 CONTINUE
819 CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
820 IF (IRES .LT. 0) RETURN
821 DO 30 K=JST,JND
822 J=ICNG(K)
823 IST=IPTR(J)
824 IND=IST+LENB(J)-1
825 DO 40 II=IST,IND
826 I=ICNB(II)
827 DELINV=1.0D0/DELY(J)
828 B(II)=(E(I)-DELTA(I))*DELINV
829 40 CONTINUE
830 Y(J)=YSAVE(J)
831 YPRIME(J)=YPSAVE(J)
832 30 CONTINUE
833 10 CONTINUE
834 RETURN
835 END
836 SUBROUTINE DDASSL (RES,NEQ,T,Y,YPRIME,TOUT,
837 * INFO,RTOL,ATOL,IDID,
838 * RWORK,LRW,IWORK,LIW,RPAR,IPAR,
839 * JAC,PTN,ICHAN)
840 C
841 C***BEGIN PROLOGUE DDASSL
842 C
843 C :
844 C :
845 C :
846 C :
847 C
848 CC***ROUTINES CALLED DDASTP,DDAINI,DDANRM,DDAWTS,DDATRP,XERRWV,D1MACH
849 C***COMMON BLOCKS DDA001
850 C***END PROLOGUE DDASSL
851 C
852 C
853 C C_2 " " " 17.9.86 SPARSE OPTION:MTYPE,LENPD,LENRW,LENIW
854 C C_3 " " " " " " EINBAU VON PTN IN PARAMETERKLAMM. U.EXT
855 C C_5 " " " 26.9.86 ERWEITERUNG VON COMMON/DDA001/
856 C C_8 R,KOENIGSDORFF 12.01.87 EINBAU VON COMMON/JACV1/
857 C C_10 A.K. 14.07.87 JACOBIDIMENSION. MIT LENJVD+LENJVS
858 C C_11 A.K. 19.11.87 fuer ICHAN=1 Neustart ohne Patterbest.
859 C C_12 A.K. 26.05.89 Adressrechng. ueberprueft, KWKN=1 da
860 C Restlaenge von RWORK uebergeben wird.
861 C
862 C
863 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
864 LOGICAL DONE
865 EXTERNAL RES,JAC,PTN
866 DIMENSION Y(1),YPRIME(1)
867 DIMENSION INFO(15)
868 DIMENSION RWORK(1),IWORK(1)
869 DIMENSION RTOL(1),ATOL(1)
870 DIMENSION RPAR(1),IPAR(1)
871 C_5
872 COMMON/DDA001/NPD,NTEMP,
873 * LML,LMU,LMXORD,LMTYPE,
874 * LNST,LNRE,LNJE,LETF,LCTF,LIPVT,
875 * KWKN,KJAC,KLUD,KIPTR,KLENB,KICNB,KLENR,KLENRL,
876 * KIPA,KIQA,KICN,JIRN,KLENC,JIFIRST,JLASTR,JNEXTR,JLASTC,JNEXTC,
877 * JIPC,IDISPX(2),LENOFF(1),JACSIZ,LUDSIZ,PT1,NCG,KICNG,KICGP,LIRN
878 C_5 END
879 C_8
880 C
881 COMMON /ANVE / LENAVX,LENAVU,LENAVP,LENB
882 C
883 COMMON / JACV / LENFX,LENBX,LENRX,LENBU
884 C COMMON / JACV1 / LENJVG,LENJVD,LENJVS,LENJAD,LENJAS,LENJMD,LENJMS
885 C_8 END
886 DATA LTSTOP,LHMAX,LH,LTN,
887 * LCJ,LCJOLD,LHOLD,LS,LROUND,
888 * LALPHA,LBETA,LGAMMA,
889 * LPSI,LSIGMA,LDELTA
890 * /1,2,3,4,
891 * 5,6,7,8,9,
892 * 11,17,23,
893 * 29,35,41/
894 IF(INFO(1).NE.0)GO TO 100
895 C
896 C-----------------------------------------------------------------------
897 C this block is executed for the initial call only.
898 C it contains checking of inputs and initializations.
899 C-----------------------------------------------------------------------
900 C
901 C first check info array to make sure all elements of info
902 C are either zero or one.
903 DO 10 I=2,11
904 IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701
905 10 CONTINUE
906 C
907 IF(NEQ.LE.0)GO TO 702
908 C
909 C set pointers into iwork
910 LML=1
911 LMU=2
912 LMXORD=3
913 LMTYPE=4
914 LJCALC=5
915 LPHASE=6
916 LK=7
917 LKOLD=8
918 LNS=9
919 LNSTL=10
920 LNST=11
921 LNRE=12
922 LNJE=13
923 LETF=14
924 LCTF=15
925 LIPVT=21
926 LIWM=1
927 C
928 C check and compute maximum order
929 MXORD=5
930 IF(INFO(9).EQ.0)GO TO 20
931 MXORD=IWORK(LMXORD)
932 IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703
933 20 IWORK(LMXORD)=MXORD
934 C
935 C compute mtype,lenpd,lenrw.check ml and mu.
936 IF(INFO(6).NE.0)GO TO 40
937 LENPD=NEQ**2
938 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
939 LENIW=20+NEQ
940 IF(INFO(5).NE.0)GO TO 30
941 IWORK(LMTYPE)=2
942 GO TO 60
943 30 IWORK(LMTYPE)=1
944 GO TO 60
945 C_2 SPARSE CASE
946 40 IWORK(LMTYPE)=3
947 NII2=1
948 C_10>
949 NJAC=LENFX + LENBX + LENB
950 C_ak160689 NLUD=2*(LENFX + LENBX + LENB)
951 NLUD=3*(LENFX + LENBX + LENB)
952 C_10<
953 MLUD=NLUD
954 LENPD=NEQ+NJAC+NLUD
955 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
956 LENIW=20+NEQ*(2+13/NII2)+(NJAC+NLUD+MLUD)/NII2
957 C
958 C check lengths of rwork and iwork
959 C_2 60 LENIW=20+NEQ
960 60 IF(LRW.LT.LENRW)GO TO 704
961 IF(LIW.LT.LENIW)GO TO 705
962 C
963 C check to see that tout is different from t
964 IF(TOUT .EQ. T)GO TO 719
965 C
966 C check hmax
967 IF(INFO(7).EQ.0)GO TO 70
968 HMAX=RWORK(LHMAX)
969 IF(HMAX.LE.0.0D0)GO TO 710
970 70 CONTINUE
971 C
972 C initialize counters
973 IWORK(LNST)=0
974 IWORK(LNRE)=0
975 C_11>
976 IF (ICHAN .EQ. 1) THEN
977 IWORK(LNJE)=1
978 ELSE
979 IWORK(LNJE)=0
980 ENDIF
981 C_11<
982 IWORK(LNSTL)=0
983 IDID=1
984 GO TO 200
985 C
986 C-----------------------------------------------------------------------
987 C this block is for continuation calls
988 C only. here we check info(1),and if the
989 C last step was interrupted we check whether
990 C appropriate action was taken.
991 C-----------------------------------------------------------------------
992 C
993 100 CONTINUE
994 IF(INFO(1).EQ.1)GO TO 110
995 IF(INFO(1).NE.-1)GO TO 701
996 C if we are here, the last step was interrupted
997 C by an error condition from ddastp,and
998 C appropriate action was not taken. this
999 C is a fatal error.
1000 CALL XERRWV(
1001 *49HDASSL-- THE LAST STEP TERMINATED WITH A NEGATIVE,
1002 *49,201,0,0,0,0,0,0.0D0,0.0D0)
1003 CALL XERRWV(
1004 *47HDASSL-- VALUE (=I1) OF IDID AND NO APPROPRIATE,
1005 *47,202,0,1,IDID,0,0,0.0D0,0.0D0)
1006 CALL XERRWV(
1007 *41HDASSL-- ACTION WAS TAKEN. RUN TERMINATED,
1008 *41,203,1,0,0,0,0,0.0D0,0.0D0)
1009 RETURN
1010 110 CONTINUE
1011 IWORK(LNSTL)=IWORK(LNST)
1012 C
1013 C-----------------------------------------------------------------------
1014 C this block is executed on all calls.
1015 C the error tolerance parameters are
1016 C checked, and the work array pointers
1017 C are set.
1018 C-----------------------------------------------------------------------
1019 C
1020 200 CONTINUE
1021 C check rtol,atol
1022 NZFLG=0
1023 RTOLI=RTOL(1)
1024 ATOLI=ATOL(1)
1025 DO 210 I=1,NEQ
1026 IF(INFO(2).EQ.1)RTOLI=RTOL(I)
1027 IF(INFO(2).EQ.1)ATOLI=ATOL(I)
1028 IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1
1029 IF(RTOLI.LT.0.0D0)GO TO 706
1030 IF(ATOLI.LT.0.0D0)GO TO 707
1031 210 CONTINUE
1032 IF(NZFLG.EQ.0)GO TO 708
1033 C
1034 C set up rwork storage.iwork storage is fixed
1035 C in data statement.
1036 LE=LDELTA+NEQ
1037 LWT=LE+NEQ
1038 LPHI=LWT+NEQ
1039 LPD=LPHI+(IWORK(LMXORD)+1)*NEQ
1040 LWM=LPD
1041 C_7
1042 IF (IWORK(LMTYPE).EQ.3) THEN
1043 LENOFF(1)=-1
1044 PT1=0.1D0
1045 JACSIZ=NJAC
1046 C_12 KWKN=LWM
1047 KWKN=1
1048 KJAC=KWKN+NEQ
1049 KLUD=KJAC+JACSIZ
1050 KIPTR=LIPVT
1051 KLENB=KIPTR+NEQ
1052 MYY=(NEQ-1)/NII2+1
1053 KICNB=KLENB+MYY
1054 KICNG=KICNB+(JACSIZ-1)/NII2+1
1055 KICGP=KICNG+MYY
1056 KLENR=KICGP+MYY+1
1057 KLENRL=MYY+KLENR
1058 KIPA=KLENRL+MYY
1059 KIQA=KIPA+MYY
1060 KICN=KIQA+MYY
1061 JIPC=LIW+1-NEQ
1062 KLENC=JIPC-MYY
1063 JIFIRST=KLENC-MYY
1064 JLASTR=JIFIRST-MYY
1065 JNEXTR=JLASTR-MYY
1066 JLASTC=JNEXTR-MYY
1067 JNEXTC=JLASTC-MYY
1068 C remaning length on IWORK: KSIZ
1069 KSIZ=JNEXTC-KICN
1070 KSIZN=(JACSIZ+1+2*NEQ)/NII2
1071 KSIZM=KSIZN
1072 IF(KSIZ.LT.2*KSIZN) KSIZM=KSIZ/2
1073 KSIZM=(KSIZM+KSIZ)/3
1074 C remaning length on RWORK: LIAN
1075 LIAN=(LRW-KLUD-LWM+1)/NII2
1076 C_12 LIAN=(LRW-KLUD+1)/NII2
1077 LUDSIZ=MIN0(KSIZ-KSIZM,LIAN)
1078 LIRN=MIN0(KSIZ-LUDSIZ,LIAN)
1079 JIRN=KICN+LUDSIZ
1080 LUDSIZ=LUDSIZ*NII2
1081 LIRN=LIRN*NII2
1082 ENDIF
1083 C_7 END
1084 NPD=1
1085 NTEMP=NPD+LENPD
1086 IF(INFO(1).EQ.1)GO TO 400
1087 C
1088 C-----------------------------------------------------------------------
1089 C this block is executed on the initial call
1090 C only. set the initial step size, and
1091 C the error weight vector, and phi.
1092 C compute initial yprime, if necessary.
1093 C-----------------------------------------------------------------------
1094 C
1095 300 CONTINUE
1096 TN=T
1097 IDID=1
1098 C
1099 C set error weight vector wt
1100 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
1101 DO 305 I = 1,NEQ
1102 IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713
1103 305 CONTINUE
1104 C
1105 C compute unit roundoff and hmin
1106 UROUND = D1MACH(4)
1107 RWORK(LROUND) = UROUND
1108 HMIN = 4.0D0*UROUND*DMAX1(DABS(T),DABS(TOUT))
1109 C
1110 C check initial interval to see that it is long enough
1111 TDIST = DABS(TOUT - T)
1112 IF(TDIST .LT. HMIN) GO TO 714
1113 C
1114 C check ho, if this was input
1115 IF (INFO(8) .EQ. 0) GO TO 310
1116 HO = RWORK(LH)
1117 IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711
1118 IF (HO .EQ. 0.0D0) GO TO 712
1119 GO TO 320
1120 310 CONTINUE
1121 C
1122 C compute initial stepsize, to be used by either
1123 C ddastp or ddaini, depending on info(11)
1124 HO = 0.001D0*TDIST
1125 YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR)
1126 IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM
1127 HO = DSIGN(HO,TOUT-T)
1128 C adjust ho if necessary to meet hmax bound
1129 320 IF (INFO(7) .EQ. 0) GO TO 330
1130 RH = DABS(HO)/HMAX
1131 IF (RH .GT. 1.0D0) HO = HO/RH
1132 C compute tstop, if applicable
1133 330 IF (INFO(4) .EQ. 0) GO TO 340
1134 TSTOP = RWORK(LTSTOP)
1135 IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715
1136 IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T
1137 IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709
1138 C
1139 C compute initial derivative, if applicable
1140 340 IF (INFO(11) .EQ. 0) GO TO 350
1141 C_020490ak
1142 C test : rwork(lphi) is the first address of (iwrok(maxord)+1)*neq elements
1143 C on rwork to hold the Nordsieck vector later on. Now 4*neq elements are
1144 C used to store old and intermediate values of Y and YPRIME during
1145 C determination of consistent intial conditions.
1146 if (iwork(lmxord) .lt. 3) then
1147 write(6,*) 'Array PHI is too small for DDAINI'
1148 write(6,*) 'maximum order should be greater than 3'
1149 stop '****DDASSL just before call to DDANINI *****'
1150 endif
1151 C_020490ak
1152 CALL DDAINI(T,Y,YPRIME,NEQ,
1153 * RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR,
1154 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1155 * RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),INFO(10),PTN)
1156 IF (IDID .LT. 0) GO TO 390
1157 C
1158 C load h with ho. store h in rwork(lh)
1159 350 H = HO
1160 RWORK(LH) = H
1161 C
1162 C load y and h*yprime into phi(*,1) and phi(*,2)
1163 360 ITEMP = LPHI + NEQ
1164 DO 370 I = 1,NEQ
1165 RWORK(LPHI + I - 1) = Y(I)
1166 370 RWORK(ITEMP + I - 1) = H*YPRIME(I)
1167 C
1168 390 GO TO 500
1169 C
1170 C-------------------------------------------------------
1171 C this block is for continuation calls only. its
1172 C purpose is to check stop conditions before
1173 C taking a step.
1174 C adjust h if necessary to meet hmax bound
1175 C-------------------------------------------------------
1176 C
1177 400 CONTINUE
1178 DONE = .FALSE.
1179 TN=RWORK(LTN)
1180 H=RWORK(LH)
1181 IF(INFO(7) .EQ. 0) GO TO 410
1182 RH = DABS(H)/HMAX
1183 IF(RH .GT. 1.0D0) H = H/RH
1184 410 CONTINUE
1185 IF(T .EQ. TOUT) GO TO 719
1186 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
1187 IF(INFO(4) .EQ. 1) GO TO 430
1188 IF(INFO(3) .EQ. 1) GO TO 420
1189 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
1190 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1191 * RWORK(LPHI),RWORK(LPSI))
1192 T=TOUT
1193 IDID = 3
1194 DONE = .TRUE.
1195 GO TO 490
1196 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490
1197 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425
1198 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1199 * RWORK(LPHI),RWORK(LPSI))
1200 T = TN
1201 IDID = 1
1202 DONE = .TRUE.
1203 GO TO 490
1204 425 CONTINUE
1205 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1206 * RWORK(LPHI),RWORK(LPSI))
1207 T = TOUT
1208 IDID = 3
1209 DONE = .TRUE.
1210 GO TO 490
1211 430 IF(INFO(3) .EQ. 1) GO TO 440
1212 TSTOP=RWORK(LTSTOP)
1213 IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
1214 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
1215 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
1216 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1217 * RWORK(LPHI),RWORK(LPSI))
1218 T=TOUT
1219 IDID = 3
1220 DONE = .TRUE.
1221 GO TO 490
1222 440 TSTOP = RWORK(LTSTOP)
1223 IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
1224 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
1225 IF((TN-T)*H .LE. 0.0D0) GO TO 450
1226 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445
1227 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1228 * RWORK(LPHI),RWORK(LPSI))
1229 T = TN
1230 IDID = 1
1231 DONE = .TRUE.
1232 GO TO 490
1233 445 CONTINUE
1234 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1235 * RWORK(LPHI),RWORK(LPSI))
1236 T = TOUT
1237 IDID = 3
1238 DONE = .TRUE.
1239 GO TO 490
1240 450 CONTINUE
1241 C check whether we are with in roundoff of tstop
1242 IF(DABS(TN-TSTOP).GT.100.0D0*UROUND*
1243 * (DABS(TN)+DABS(H)))GO TO 460
1244 IDID=2
1245 T=TSTOP
1246 DONE = .TRUE.
1247 GO TO 490
1248 460 TNEXT=TN+H*(1.0D0+4.0D0*UROUND)
1249 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
1250 H=(TSTOP-TN)*(1.0D0-4.0D0*UROUND)
1251 RWORK(LH)=H
1252 C
1253 490 IF (DONE) GO TO 590
1254 C
1255 C-------------------------------------------------------
1256 C the next block contains the call to the
1257 C one-step integrator ddastp.
1258 C this is a looping point for the integration
1259 C steps.
1260 C check for too many steps.
1261 C update wt.
1262 C check for too much accuracy requested.
1263 C compute minimum stepsize.
1264 C-------------------------------------------------------
1265 C
1266 500 CONTINUE
1267 C check for failure to compute initial yprime
1268 IF (IDID .EQ. -12) GO TO 527
1269 C
1270 C check for too many steps
1271 IF((IWORK(LNST)-IWORK(LNSTL)).LT.500)
1272 * GO TO 510
1273 IDID=-1
1274 GO TO 527
1275 C
1276 C update wt
1277 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),
1278 * RWORK(LWT),RPAR,IPAR)
1279 DO 520 I=1,NEQ
1280 IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520
1281 IDID=-3
1282 GO TO 527
1283 520 CONTINUE
1284 C
1285 C test for too much accuracy requested.
1286 R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*
1287 * 100.0D0*UROUND
1288 IF(R.LE.1.0D0)GO TO 525
1289 C multiply rtol and atol by r and return
1290 IF(INFO(2).EQ.1)GO TO 523
1291 RTOL(1)=R*RTOL(1)
1292 ATOL(1)=R*ATOL(1)
1293 IDID=-2
1294 GO TO 527
1295 523 DO 524 I=1,NEQ
1296 RTOL(I)=R*RTOL(I)
1297 524 ATOL(I)=R*ATOL(I)
1298 IDID=-2
1299 GO TO 527
1300 525 CONTINUE
1301 C
1302 C compute minimum stepsize
1303 HMIN=4.0D0*UROUND*DMAX1(DABS(TN),DABS(TOUT))
1304 C
1305 CALL DDASTP(TN,Y,YPRIME,NEQ,
1306 * RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR,
1307 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1308 * RWORK(LWM),IWORK(LIWM),
1309 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
1310 * RWORK(LPSI),RWORK(LSIGMA),
1311 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),
1312 * RWORK(LS),HMIN,RWORK(LROUND),
1313 * IWORK(LPHASE),IWORK(LJCALC),IWORK(LK),
1314 * IWORK(LKOLD),IWORK(LNS),INFO(10),PTN)
1315 527 IF(IDID.LT.0)GO TO 600
1316 C
1317 C------------------------------------------------------
1318 C this block handles the case of a successful
1319 C return from ddastp (idid=1) test for
1320 C stop conditions.
1321 C--------------------------------------------------------
1322 C
1323 IF(INFO(4).NE.0)GO TO 540
1324 IF(INFO(3).NE.0)GO TO 530
1325 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
1326 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1327 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1328 IDID=3
1329 T=TOUT
1330 GO TO 580
1331 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
1332 T=TN
1333 IDID=1
1334 GO TO 580
1335 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1336 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1337 IDID=3
1338 T=TOUT
1339 GO TO 580
1340 540 IF(INFO(3).NE.0)GO TO 550
1341 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
1342 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1343 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1344 T=TOUT
1345 IDID=3
1346 GO TO 580
1347 542 IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*
1348 * (DABS(TN)+DABS(H)))GO TO 545
1349 TNEXT=TN+H*(1.0D0+4.0D0*UROUND)
1350 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
1351 H=(TSTOP-TN)*(1.0D0-4.0D0*UROUND)
1352 GO TO 500
1353 545 IDID=2
1354 T=TSTOP
1355 GO TO 580
1356 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
1357 IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*(DABS(TN)+DABS(H)))GO TO 552
1358 T=TN
1359 IDID=1
1360 GO TO 580
1361 552 IDID=2
1362 T=TSTOP
1363 GO TO 580
1364 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1365 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1366 T=TOUT
1367 IDID=3
1368 580 CONTINUE
1369 C
1370 C--------------------------------------------------------
1371 C all successful returns from ddassl are made from
1372 C this block.
1373 C--------------------------------------------------------
1374 C
1375 590 CONTINUE
1376 RWORK(LTN)=TN
1377 RWORK(LH)=H
1378 RETURN
1379 C
1380 C-----------------------------------------------------------------------
1381 C this block handles all unsuccessful
1382 C returns other than for illegal input.
1383 C-----------------------------------------------------------------------
1384 C
1385 600 CONTINUE
1386 ITEMP=-IDID
1387 GO TO (610,620,630,690,690,640,650,660,670,675,
1388 * 680,685), ITEMP
1389 C
1390 C the maximum number of steps was taken before
1391 C reaching tout
1392 610 CALL XERRWV(
1393 *38HDASSL-- AT CURRENT T (=R1) 500 STEPS,
1394 *38,610,0,0,0,0,1,TN,0.0D0)
1395 CALL XERRWV(48HDASSL-- TAKEN ON THIS CALL BEFORE REACHING TOUT,
1396 *48,611,0,0,0,0,0,0.0D0,0.0D0)
1397 GO TO 690
1398 C
1399 C too much accuracy for machine precision
1400 620 CALL XERRWV(
1401 *47HDASSL-- AT T (=R1) TOO MUCH ACCURACY REQUESTED,
1402 *47,620,0,0,0,0,1,TN,0.0D0)
1403 CALL XERRWV(
1404 *48HDASSL-- FOR PRECISION OF MACHINE. RTOL AND ATOL,
1405 *48,621,0,0,0,0,0,0.0D0,0.0D0)
1406 CALL XERRWV(
1407 *45HDASSL-- WERE INCREASED TO APPROPRIATE VALUES,
1408 *45,622,0,0,0,0,0,0.0D0,0.0D0)
1409 C
1410 GO TO 690
1411 C wt(i) .le. 0.0d0 for some i (not at start of problem)
1412 630 CALL XERRWV(
1413 *38HDASSL-- AT T (=R1) SOME ELEMENT OF WT,
1414 *38,630,0,0,0,0,1,TN,0.0D0)
1415 CALL XERRWV(28HDASSL-- HAS BECOME .LE. 0.0,
1416 *28,631,0,0,0,0,0,0.0D0,0.0D0)
1417 GO TO 690
1418 C
1419 C error test failed repeatedly or with h=hmin
1420 640 CALL XERRWV(
1421 *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE,
1422 *44,640,0,0,0,0,2,TN,H)
1423 CALL XERRWV(
1424 *57HDASSL-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN,
1425 *57,641,0,0,0,0,0,0.0D0,0.0D0)
1426 GO TO 690
1427 C
1428 C corrector convergence failed repeatedly or with h=hmin
1429 650 CALL XERRWV(
1430 *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE,
1431 *44,650,0,0,0,0,2,TN,H)
1432 CALL XERRWV(
1433 *48HDASSL-- CORRECTOR FAILED TO CONVERGE REPEATEDLY,
1434 *48,651,0,0,0,0,0,0.0D0,0.0D0)
1435 CALL XERRWV(
1436 *28HDASSL-- OR WITH ABS(H)=HMIN,
1437 *28,652,0,0,0,0,0,0.0D0,0.0D0)
1438 GO TO 690
1439 C
1440 C the iteration matrix is singular
1441 660 CALL XERRWV(
1442 *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE,
1443 *44,660,0,0,0,0,2,TN,H)
1444 CALL XERRWV(
1445 *37HDASSL-- ITERATION MATRIX IS SINGULAR,
1446 *37,661,0,0,0,0,0,0.0D0,0.0D0)
1447 GO TO 690
1448 C
1449 C corrector failure preceeded by error test failures.
1450 670 CALL XERRWV(
1451 *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE,
1452 *44,670,0,0,0,0,2,TN,H)
1453 CALL XERRWV(
1454 *49HDASSL-- CORRECTOR COULD NOT CONVERGE. ALSO, THE,
1455 *49,671,0,0,0,0,0,0.0D0,0.0D0)
1456 CALL XERRWV(
1457 *38HDASSL-- ERROR TEST FAILED REPEATEDLY.,
1458 *38,672,0,0,0,0,0,0.0D0,0.0D0)
1459 GO TO 690
1460 C
1461 C corrector failure because ires = -1
1462 675 CALL XERRWV(
1463 *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE,
1464 *44,675,0,0,0,0,2,TN,H)
1465 CALL XERRWV(
1466 *45HDASSL-- CORRECTOR COULD NOT CONVERGE BECAUSE,
1467 *455,676,0,0,0,0,0,0.0D0,0.0D0)
1468 CALL XERRWV(
1469 *36HDASSL-- IRES WAS EQUAL TO MINUS ONE,
1470 *36,677,0,0,0,0,0,0.0D0,0.0D0)
1471 GO TO 690
1472 C
1473 C failure because ires = -2
1474 680 CALL XERRWV(
1475 *40HDASSL-- AT T (=R1) AND STEPSIZE H (=R2),
1476 *40,680,0,0,0,0,2,TN,H)
1477 CALL XERRWV(
1478 *36HDASSL-- IRES WAS EQUAL TO MINUS TWO,
1479 *36,681,0,0,0,0,0,0.0D0,0.0D0)
1480 GO TO 690
1481 C
1482 C failed to compute initial yprime
1483 685 CALL XERRWV(
1484 *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE,
1485 *44,685,0,0,0,0,2,TN,HO)
1486 CALL XERRWV(
1487 *45HDASSL-- INITIAL YPRIME COULD NOT BE COMPUTED,
1488 *45,686,0,0,0,0,0,0.0D0,0.0D0)
1489 GO TO 690
1490 690 CONTINUE
1491 INFO(1)=-1
1492 T=TN
1493 RWORK(LTN)=TN
1494 RWORK(LH)=H
1495 RETURN
1496 C-----------------------------------------------------------------------
1497 C this block handles all error returns due
1498 C to illegal input, as detected before calling
1499 C ddastp. first the error message routine is
1500 C called. if this happens twice in
1501 C succession, execution is terminated
1502 C
1503 C-----------------------------------------------------------------------
1504 701 CALL XERRWV(
1505 *55HDASSL-- SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE,
1506 *55,1,0,0,0,0,0,0.0D0,0.0D0)
1507 GO TO 750
1508 702 CALL XERRWV(25HDASSL-- NEQ (=I1) .LE. 0,
1509 *25,2,0,1,NEQ,0,0,0.0D0,0.0D0)
1510 GO TO 750
1511 703 CALL XERRWV(34HDASSL-- MAXORD (=I1) NOT IN RANGE,
1512 *34,3,0,1,MXORD,0,0,0.0D0,0.0D0)
1513 GO TO 750
1514 704 CALL XERRWV(
1515 *60HDASSL-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2),
1516 *60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)
1517 GO TO 750
1518 705 CALL XERRWV(
1519 *60HDASSL-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2),
1520 *60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)
1521 GO TO 750
1522 706 CALL XERRWV(
1523 *39HDASSL-- SOME ELEMENT OF RTOL IS .LT. 0,
1524 *39,6,0,0,0,0,0,0.0D0,0.0D0)
1525 GO TO 750
1526 707 CALL XERRWV(
1527 *39HDASSL-- SOME ELEMENT OF ATOL IS .LT. 0,
1528 *39,7,0,0,0,0,0,0.0D0,0.0D0)
1529 GO TO 750
1530 708 CALL XERRWV(
1531 *47HDASSL-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO,
1532 *47,8,0,0,0,0,0,0.0D0,0.0D0)
1533 GO TO 750
1534 709 CALL XERRWV(
1535 *54HDASSL-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2),
1536 *54,9,0,0,0,0,2,TSTOP,TOUT)
1537 GO TO 750
1538 710 CALL XERRWV(28HDASSL-- HMAX (=R1) .LT. 0.0,
1539 *28,10,0,0,0,0,1,HMAX,0.0D0)
1540 GO TO 750
1541 711 CALL XERRWV(34HDASSL-- TOUT (=R1) BEHIND T (=R2),
1542 *34,11,0,0,0,0,2,TOUT,T)
1543 GO TO 750
1544 712 CALL XERRWV(29HDASSL-- INFO(8)=1 AND H0=0.0,
1545 *29,12,0,0,0,0,0,0.0D0,0.0D0)
1546 GO TO 750
1547 713 CALL XERRWV(39HDASSL-- SOME ELEMENT OF WT IS .LE. 0.0,
1548 *39,13,0,0,0,0,0,0.0D0,0.0D0)
1549 GO TO 750
1550 714 CALL XERRWV(
1551 *61HDASSL-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION,
1552 *61,14,0,0,0,0,2,TOUT,T)
1553 GO TO 750
1554 715 CALL XERRWV(
1555 *49HDASSL-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2),
1556 *49,15,0,0,0,0,2,TSTOP,T)
1557 GO TO 750
1558 719 CALL XERRWV(
1559 *39HDASSL-- TOUT (=R1) IS EQUAL TO T (=R2),
1560 *39,19,0,0,0,0,2,TOUT,T)
1561 GO TO 750
1562 750 IF(INFO(1).EQ.-1) GO TO 760
1563 INFO(1)=-1
1564 IDID=-33
1565 RETURN
1566 760 CALL XERRWV(
1567 *46HDASSL-- REPEATED OCCURRENCES OF ILLEGAL INPUT,
1568 *46,801,0,0,0,0,0,0.0D0,0.0D0)
1569 770 CALL XERRWV(
1570 *47HDASSL-- RUN TERMINATED. APPARENT INFINITE LOOP,
1571 *47,802,1,0,0,0,0,0.0D0,0.0D0)
1572 RETURN
1573 C-----------end of subroutine ddassl-------------------------------------
1574 END
1575 SUBROUTINE DDASTP(X,Y,YPRIME,NEQ,
1576 * RES,JAC,H,WT,JSTART,IDID,RPAR,IPAR,
1577 * PHI,DELTA,E,WM,IWM,
1578 * ALPHA,BETA,GAMMA,PSI,SIGMA,
1579 * CJ,CJOLD,HOLD,S,HMIN,UROUND,
1580 * IPHASE,JCALC,K,KOLD,NS,NONNEG,PTN)
1581 C
1582 C***BEGIN PROLOGUE DDASTP
1583 C***REFER TO DDASSL
1584 C***ROUTINES CALLED DDANRM,DDAJAC,DDASLV,DDATRP
1585 C***COMMON BLOCKS DDA001
1586 C***DATE WRITTEN 830315 (YYMMDD)
1587 C***REVISION DATE 830315 (YYMMDD)
1588 C***END PROLOGUE DDASTP
1589 C
1590 C
1591 C-----------------------------------------------------------------------
1592 C dastep solves a system of differential/
1593 C algebraic equations of the form
1594 C g(x,y,yprime) = 0, for one step (normally
1595 C from x to x+h).
1596 C
1597 C the methods used are modified divided
1598 C difference,fixed leading coefficient
1599 C forms of backward differentiation
1600 C formulas. the code adjusts the stepsize
1601 C and order to control the local error per
1602 C step.
1603 C
1604 C
1605 C the parameters represent
1606 C x -- independent variable
1607 C y -- solution vector at x
1608 C yprime -- derivative of solution vector
1609 C after successful step
1610 C neq -- number of equations to be integrated
1611 C res -- external user-supplied subroutine
1612 C to evaluate the residual. the call is
1613 C call res(x,y,yprime,delta,ires,rpar,ipar)
1614 C x,y,yprime are input. delta is output.
1615 C on input, ires=0. res should alter ires only
1616 C if it encounters an illegal value of y or a
1617 C stop condition. set ires=-1 if an input value
1618 C of y is illegal, and dastep will try to solve
1619 C the problem without getting ires = -1. if
1620 C ires=-2, dastep returns control to the calling
1621 C program with idid = -11.
1622 C jac -- external user-supplied routine to evaluate
1623 C the iteration matrix (this is optional)
1624 C the call is of the form
1625 C call jac(x,y,yprime,pd,cj,rpar,ipar)
1626 C pd is the matrix of partial derivatives,
1627 C pd=dg/dy+cj*dg/dyprime
1628 C h -- appropriate step size for next step.
1629 C normally determined by the code
1630 C wt -- vector of weights for error criterion.
1631 C jstart -- integer variable set 0 for
1632 C first step, 1 otherwise.
1633 C idid -- completion code with the following meanings%
1634 C idid= 1 -- the step was completed successfully
1635 C idid=-6 -- the error test failed repeatedly
1636 C idid=-7 -- the corrector could not converge
1637 C idid=-8 -- the iteration matrix is singular
1638 C idid=-9 -- the corrector could not converge.
1639 C there were repeated error test
1640 C failures on this step.
1641 C idid=-10-- the corrector could not converge
1642 C because ires was equal to minus one
1643 C idid=-11-- ires equal to -2 was encountered,
1644 C and control is being returned to
1645 C the calling program
1646 C rpar,ipar -- real and integer parameter arrays that
1647 C are used for communication between the
1648 C calling program and external user routines
1649 C they are not altered by dastep
1650 C phi -- array of divided differences used by
1651 C dastep. the length is neq*(k+1),where
1652 C k is the maximum order
1653 C delta,e -- work vectors for dastep of length neq
1654 C wm,iwm -- real and integer arrays storing
1655 C matrix information such as the matrix
1656 C of partial derivatives,permutation
1657 C vector,and various other information.
1658 C
1659 C the other parameters are information
1660 C which is needed internally by dastep to
1661 C continue from step to step.
1662 C
1663 C-----------------------------------------------------------------------
1664 C_1 R.KOE EINBAU VON PTN IN PARAMETERKAMMERN DDASTP,DDAJAC U.EXT
1665 C
1666 C
1667 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1668 LOGICAL CONVGD
1669 DIMENSION Y(1),YPRIME(1),WT(1)
1670 DIMENSION PHI(NEQ,1),DELTA(1),E(1)
1671 DIMENSION WM(1),IWM(1)
1672 DIMENSION PSI(1),ALPHA(1),BETA(1),GAMMA(1),SIGMA(1)
1673 DIMENSION RPAR(1),IPAR(1)
1674 EXTERNAL RES,JAC,PTN
1675 COMMON/DDA001/NPD,NTEMP,
1676 * LML,LMU,LMXORD,LMTYPE,
1677 * LNST,LNRE,LNJE,LETF,LCTF,LIPVT
1678 DATA MAXIT/4/
1679 DATA XRATE/0.25D0/
1680 C
1681 C
1682 C
1683 C
1684 C
1685 C-----------------------------------------------------------------------
1686 C block 1.
1687 C initialize. on the first call,set
1688 C the order to 1 and initialize
1689 C other variables.
1690 C-----------------------------------------------------------------------
1691 C
1692 C initializations for all calls
1693 IDID=1
1694 XOLD=X
1695 NCF=0
1696 NSF=0
1697 NEF=0
1698 IF(JSTART .NE. 0) GO TO 120
1699 C
1700 C if this is the first step,perform
1701 C other initializations
1702 C AUCH BEI CONTINUATION CALLS!!!!!
1703 IWM(LETF) = 0
1704 IWM(LCTF) = 0
1705 K=1
1706 KOLD=0
1707 HOLD=0.0D0
1708 JSTART=1
1709 PSI(1)=H
1710 CJOLD = 1.0D0/H
1711 CJ = CJOLD
1712 S = 100.D0
1713 JCALC = -1
1714 DELNRM=1.0D0
1715 IPHASE = 0
1716 NS=0
1717 120 CONTINUE
1718 C
1719 C
1720 C
1721 C
1722 C
1723 C
1724 C-----------------------------------------------------------------------
1725 C block 2
1726 C compute coefficients of formulas for
1727 C this step.
1728 C-----------------------------------------------------------------------
1729 200 CONTINUE
1730 KP1=K+1
1731 KP2=K+2
1732 KM1=K-1
1733 XOLD=X
1734 IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
1735 NS=MIN0(NS+1,KOLD+2)
1736 NSP1=NS+1
1737 IF(KP1 .LT. NS)GO TO 230
1738 C
1739 BETA(1)=1.0D0
1740 ALPHA(1)=1.0D0
1741 TEMP1=H
1742 GAMMA(1)=0.0D0
1743 SIGMA(1)=1.0D0
1744 DO 210 I=2,KP1
1745 TEMP2=PSI(I-1)
1746 PSI(I-1)=TEMP1
1747 BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
1748 TEMP1=TEMP2+H
1749 ALPHA(I)=H/TEMP1
1750 SIGMA(I)=DFLOAT(I-1)*SIGMA(I-1)*ALPHA(I)
1751 GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
1752 210 CONTINUE
1753 PSI(KP1)=TEMP1
1754 230 CONTINUE
1755 C
1756 C compute alphas, alpha0
1757 ALPHAS = 0.0D0
1758 ALPHA0 = 0.0D0
1759 DO 240 I = 1,K
1760 ALPHAS = ALPHAS - 1.0D0/DFLOAT(I)
1761 ALPHA0 = ALPHA0 - ALPHA(I)
1762 240 CONTINUE
1763 C
1764 C compute leading coefficient cj
1765 CJLAST = CJ
1766 CJ = -ALPHAS/H
1767 C
1768 C compute variable stepsize error coefficient ck
1769 CK = DABS(ALPHA(KP1) + ALPHAS - ALPHA0)
1770 CK = DMAX1(CK,ALPHA(KP1))
1771 C
1772 C decide whether new jacobian is needed
1773 TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
1774 TEMP2 = 1.0D0/TEMP1
1775 IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
1776 IF (CJ .NE. CJLAST) S = 100.D0
1777 C
1778 C change phi to phi star
1779 IF(KP1 .LT. NSP1) GO TO 280
1780 DO 270 J=NSP1,KP1
1781 DO 260 I=1,NEQ
1782 260 PHI(I,J)=BETA(J)*PHI(I,J)
1783 270 CONTINUE
1784 280 CONTINUE
1785 C
1786 C update time
1787 X=X+H
1788 C
1789 C
1790 C
1791 C
1792 C
1793 C-----------------------------------------------------------------------
1794 C block 3
1795 C predict the solution and derivative,
1796 C and solve the corrector equation
1797 C-----------------------------------------------------------------------
1798 C
1799 C first,predict the solution and derivative
1800 300 CONTINUE
1801 DO 310 I=1,NEQ
1802 Y(I)=PHI(I,1)
1803 310 YPRIME(I)=0.0D0
1804 DO 330 J=2,KP1
1805 DO 320 I=1,NEQ
1806 Y(I)=Y(I)+PHI(I,J)
1807 320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
1808 330 CONTINUE
1809 PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR)
1810 C
1811 C
1812 C
1813 C solve the corrector equation using a
1814 C modified newton scheme.
1815 CONVGD= .TRUE.
1816 M=0
1817 IWM(LNRE)=IWM(LNRE)+1
1818 IRES = 0
1819 CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
1820 IF (IRES .LT. 0) GO TO 380
1821 C
1822 C
1823 C if indicated,reevaluate the
1824 C iteration matrix pd = dg/dy + cj*dg/dyprime
1825 C (where g(x,y,yprime)=0). set
1826 C jcalc to 0 as an indicator that
1827 C this has been done.
1828 IF(JCALC .NE. -1)GO TO 340
1829 IWM(LNJE)=IWM(LNJE)+1
1830 JCALC=0
1831 CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
1832 * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,IPAR,PTN)
1833 CJOLD=CJ
1834 S = 100.D0
1835 IF (IRES .LT. 0) GO TO 380
1836 IF(IER .NE. 0)GO TO 380
1837 NSF=0
1838 C
1839 C
1840 C initialize the error accumulation vector e.
1841 340 CONTINUE
1842 DO 345 I=1,NEQ
1843 345 E(I)=0.0D0
1844 C
1845 S = 100.E0
1846 C
1847 C
1848 C corrector loop.
1849 350 CONTINUE
1850 C
1851 C multiply residual by temp1 to accelerate convergence
1852 TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
1853 DO 355 I = 1,NEQ
1854 355 DELTA(I) = DELTA(I) * TEMP1
1855 C
1856 C compute a new iterate (back-substitution).
1857 C store the correction in delta.
1858 CALL DDASLV(NEQ,DELTA,WM,IWM)
1859 C
1860 C update y,e,and yprime
1861 DO 360 I=1,NEQ
1862 Y(I)=Y(I)-DELTA(I)
1863 E(I)=E(I)-DELTA(I)
1864 360 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
1865 C
1866 C test for convergence of the iteration
1867 DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
1868 IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375
1869 IF (M .GT. 0) GO TO 365
1870 OLDNRM = DELNRM
1871 GO TO 367
1872 365 RATE = (DELNRM/OLDNRM)**(1.0D0/DFLOAT(M))
1873 IF (RATE .GT. 0.90D0) GO TO 370
1874 S = RATE/(1.0D0 - RATE)
1875 367 IF (S*DELNRM .LE. 0.33D0) GO TO 375
1876 C
1877 C the corrector has not yet converged.
1878 C update m and test whether the
1879 C maximum number of iterations have
1880 C been tried.
1881 M=M+1
1882 IF(M.GE.MAXIT)GO TO 370
1883 C
1884 C evaluate the residual
1885 C and go back to do another iteration
1886 IWM(LNRE)=IWM(LNRE)+1
1887 IRES = 0
1888 CALL RES(X,Y,YPRIME,DELTA,IRES,
1889 * RPAR,IPAR)
1890 IF (IRES .LT. 0) GO TO 380
1891 GO TO 350
1892 C
1893 C
1894 C the corrector failed to converge in maxit
1895 C iterations. if the iteration matrix
1896 C is not current,re-do the step with
1897 C a new iteration matrix.
1898 370 CONTINUE
1899 IF(JCALC.EQ.0)GO TO 380
1900 JCALC=-1
1901 GO TO 300
1902 C
1903 C
1904 C the iteration has converged. if nonnegativity of solution is
1905 C required, set the solution nonnegative, if the perturbation
1906 C to do it is small enough. if the change is too large, then
1907 C consider the corrector iteration to have failed.
1908 375 IF(NONNEG .EQ. 0) GO TO 390
1909 DO 377 I = 1,NEQ
1910 377 DELTA(I) = DMIN1(Y(I),0.0D0)
1911 DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
1912 IF(DELNRM .GT. 0.33D0) GO TO 380
1913 DO 378 I = 1,NEQ
1914 378 E(I) = E(I) - DELTA(I)
1915 GO TO 390
1916 C
1917 C
1918 C exits from block 3
1919 C no convergence with current iteration
1920 C matrix,or singular iteration matrix
1921 380 CONVGD= .FALSE.
1922 390 JCALC = 1
1923 IF(.NOT.CONVGD)GO TO 600
1924 C
1925 C
1926 C
1927 C
1928 C
1929 C-----------------------------------------------------------------------
1930 C block 4
1931 C estimate the errors at orders k,k-1,k-2
1932 C as if constant stepsize was used. estimate
1933 C the local error at order k and test
1934 C whether the current step is successful.
1935 C-----------------------------------------------------------------------
1936 C
1937 C estimate errors at orders k,k-1,k-2
1938 ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR)
1939 ERK = SIGMA(K+1)*ENORM
1940 TERK = FLOAT(K+1)*ERK
1941 EST = ERK
1942 KNEW=K
1943 IF(K .EQ. 1)GO TO 430
1944 DO 405 I = 1,NEQ
1945 405 DELTA(I) = PHI(I,KP1) + E(I)
1946 ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
1947 TERKM1 = FLOAT(K)*ERKM1
1948 IF(K .GT. 2)GO TO 410
1949 IF(TERKM1 .LE. 0.5*TERK)GO TO 420
1950 GO TO 430
1951 410 CONTINUE
1952 DO 415 I = 1,NEQ
1953 415 DELTA(I) = PHI(I,K) + DELTA(I)
1954 ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
1955 TERKM2 = FLOAT(K-1)*ERKM2
1956 IF(DMAX1(TERKM1,TERKM2).GT.TERK)GO TO 430
1957 C lower the order
1958 420 CONTINUE
1959 KNEW=K-1
1960 EST = ERKM1
1961 C
1962 C
1963 C calculate the local error for the current step
1964 C to see if the step was successful
1965 430 CONTINUE
1966 ERR = CK * ENORM
1967 IF(ERR .GT. 1.0D0)GO TO 600
1968 C
1969 C
1970 C
1971 C
1972 C
1973 C-----------------------------------------------------------------------
1974 C block 5
1975 C the step is successful. determine
1976 C the best order and stepsize for
1977 C the next step. update the differences
1978 C for the next step.
1979 C-----------------------------------------------------------------------
1980 IDID=1
1981 IWM(LNST)=IWM(LNST)+1
1982 KDIFF=K-KOLD
1983 KOLD=K
1984 HOLD=H
1985 C
1986 C
1987 C estimate the error at order k+1 unless%
1988 C already decided to lower order, or
1989 C already using maximum order, or
1990 C stepsize not constant, or
1991 C order raised in previous step
1992 IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
1993 IF(IPHASE .EQ. 0)GO TO 545
1994 IF(KNEW.EQ.KM1)GO TO 540
1995 IF(K.EQ.IWM(LMXORD)) GO TO 550
1996 IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
1997 DO 510 I=1,NEQ
1998 510 DELTA(I)=E(I)-PHI(I,KP2)
1999 ERKP1 = (1.0D0/DFLOAT(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
2000 TERKP1 = FLOAT(K+2)*ERKP1
2001 IF(K.GT.1)GO TO 520
2002 IF(TERKP1.GE.0.5D0*TERK)GO TO 550
2003 GO TO 530
2004 520 IF(TERKM1.LE.DMIN1(TERK,TERKP1))GO TO 540
2005 IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
2006 C
2007 C raise order
2008 530 K=KP1
2009 EST = ERKP1
2010 GO TO 550
2011 C
2012 C lower order
2013 540 K=KM1
2014 EST = ERKM1
2015 GO TO 550
2016 C
2017 C if iphase = 0, increase order by one and multiply stepsize by
2018 C factor two
2019 545 K = KP1
2020 HNEW = H*2.0D0
2021 H = HNEW
2022 GO TO 575
2023 C
2024 C
2025 C determine the appropriate stepsize for
2026 C the next step.
2027 550 HNEW=H
2028 TEMP2=K+1
2029 R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
2030 IF(R .LT. 2.0D0) GO TO 555
2031 HNEW = 2.0D0*H
2032 GO TO 560
2033 555 IF(R .GT. 1.0D0) GO TO 560
2034 R = DMAX1(0.5D0,DMIN1(0.9D0,R))
2035 HNEW = H*R
2036 560 H=HNEW
2037 C
2038 C
2039 C update differences for next step
2040 575 CONTINUE
2041 IF(KOLD.EQ.IWM(LMXORD))GO TO 585
2042 DO 580 I=1,NEQ
2043 580 PHI(I,KP2)=E(I)
2044 585 CONTINUE
2045 DO 590 I=1,NEQ
2046 590 PHI(I,KP1)=PHI(I,KP1)+E(I)
2047 DO 595 J1=2,KP1
2048 J=KP1-J1+1
2049 DO 595 I=1,NEQ
2050 595 PHI(I,J)=PHI(I,J)+PHI(I,J+1)
2051 RETURN
2052 C
2053 C
2054 C
2055 C
2056 C
2057 C-----------------------------------------------------------------------
2058 C block 6
2059 C the step is unsuccessful. restore x,psi,phi
2060 C determine appropriate stepsize for
2061 C continuing the integration, or exit with
2062 C an error flag if there have been many
2063 C failures.
2064 C-----------------------------------------------------------------------
2065 600 IPHASE = 1
2066 C
2067 C restore x,phi,psi
2068 X=XOLD
2069 IF(KP1.LT.NSP1)GO TO 630
2070 DO 620 J=NSP1,KP1
2071 TEMP1=1.0D0/BETA(J)
2072 DO 610 I=1,NEQ
2073 610 PHI(I,J)=TEMP1*PHI(I,J)
2074 620 CONTINUE
2075 630 CONTINUE
2076 DO 640 I=2,KP1
2077 640 PSI(I-1)=PSI(I)-H
2078 C
2079 C
2080 C test whether failure is due to corrector iteration
2081 C or error test
2082 IF(CONVGD)GO TO 660
2083 IWM(LCTF)=IWM(LCTF)+1
2084 C
2085 C
2086 C the newton iteration failed to converge with
2087 C a current iteration matrix. determine the cause
2088 C of the failure and take appropriate action.
2089 IF(IER.EQ.0)GO TO 650
2090 C
2091 C the iteration matrix is singular. reduce
2092 C the stepsize by a factor of 4. if
2093 C this happens three times in a row on
2094 C the same step, return with an error flag
2095 NSF=NSF+1
2096 R = 0.25D0
2097 H=H*R
2098 IF (NSF .LT. 3 .AND. DABS(H) .GE. HMIN) GO TO 690
2099 IDID=-8
2100 GO TO 675
2101 C
2102 C
2103 C the newton iteration failed to converge for a reason
2104 C other than a singular iteration matrix. if ires = -2, then
2105 C return. otherwise, reduce the stepsize and try again, unless
2106 C too many failures have occured.
2107 650 CONTINUE
2108 IF (IRES .GT. -2) GO TO 655
2109 IDID = -11
2110 GO TO 675
2111 655 NCF = NCF + 1
2112 R = 0.25D0
2113 H = H*R
2114 IF (NCF .LT. 10 .AND. DABS(H) .GE. HMIN) GO TO 690
2115 IDID = -7
2116 IF (IRES .LT. 0) IDID = -10
2117 IF (NEF .GE. 3) IDID = -9
2118 GO TO 675
2119 C
2120 C
2121 C the newton scheme converged,and the cause
2122 C of the failure was the error estimate
2123 C exceeding the tolerance.
2124 660 NEF=NEF+1
2125 IWM(LETF)=IWM(LETF)+1
2126 IF (NEF .GT. 1) GO TO 665
2127 C
2128 C on first error test failure, keep current order or lower
2129 C order by one. compute new stepsize based on differences
2130 C of the solution.
2131 K = KNEW
2132 TEMP2 = K + 1
2133 R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
2134 R = DMAX1(0.25D0,DMIN1(0.9D0,R))
2135 H = H*R
2136 IF (DABS(H) .GE. HMIN) GO TO 690
2137 IDID = -6
2138 GO TO 675
2139 C
2140 C on second error test failure, use the current order or
2141 C decrease order by one. reduce the stepsize by a factor of
2142 C one quarter.
2143 665 IF (NEF .GT. 2) GO TO 670
2144 K = KNEW
2145 H = 0.25D0*H
2146 IF (DABS(H) .GE. HMIN) GO TO 690
2147 IDID = -6
2148 GO TO 675
2149 C
2150 C on third and subsequent error test failures, set the order to
2151 C one and reduce the stepsize by a factor of one quarter
2152 670 K = 1
2153 H = 0.25D0*H
2154 IF (DABS(H) .GE. HMIN) GO TO 690
2155 IDID = -6
2156 GO TO 675
2157 C
2158 C
2159 C
2160 C
2161 C for all crashes, restore y to its last value,
2162 C interpolate to find yprime at last x, and return
2163 675 CONTINUE
2164 CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
2165 RETURN
2166 C
2167 C
2168 C go back and try this step again
2169 690 GO TO 200
2170 C
2171 C------end of subroutine dastep------
2172 END
2173 SUBROUTINE DDATRP(X,XOUT,YOUT,YPOUT,NEQ,KOLD,PHI,PSI)
2174 C
2175 C***BEGIN PROLOGUE DDATRP
2176 C***REFER TO DDASSL
2177 C***ROUTINES CALLED (NONE)
2178 C***DATE WRITTEN 830315 (YYMMDD)
2179 C***REVISION DATE 830315 (YYMMDD)
2180 C***END PROLOGUE DDATRP
2181 C
2182 C-----------------------------------------------------------------------
2183 C the methods in subroutine dastep use polynomials
2184 C to approximate the solution. ddatrp approximates the
2185 C solution and its derivative at time xout by evaluating
2186 C one of these polynomials,and its derivative,there.
2187 C information defining this polynomial is passed from
2188 C dastep, so ddatrp cannot be used alone.
2189 C
2190 C the parameters are%
2191 C x the current time in the integration.
2192 C xout the time at which the solution is desired
2193 C yout the interpolated approximation to y at xout
2194 C (this is output)
2195 C ypout the interpolated approximation to yprime at xout
2196 C (this is output)
2197 C neq number of equations
2198 C kold order used on last successful step
2199 C phi array of scaled divided differences of y
2200 C psi array of past stepsize history
2201 C-----------------------------------------------------------------------
2202 C
2203 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2204 DIMENSION YOUT(1),YPOUT(1)
2205 DIMENSION PHI(NEQ,1),PSI(1)
2206 KOLDP1=KOLD+1
2207 TEMP1=XOUT-X
2208 DO 10 I=1,NEQ
2209 YOUT(I)=PHI(I,1)
2210 10 YPOUT(I)=0.0D0
2211 C=1.0D0
2212 D=0.0D0
2213 GAMMA=TEMP1/PSI(1)
2214 DO 30 J=2,KOLDP1
2215 D=D*GAMMA+C/PSI(J-1)
2216 C=C*GAMMA
2217 GAMMA=(TEMP1+PSI(J-1))/PSI(J)
2218 DO 20 I=1,NEQ
2219 YOUT(I)=YOUT(I)+C*PHI(I,J)
2220 20 YPOUT(I)=YPOUT(I)+D*PHI(I,J)
2221 30 CONTINUE
2222 RETURN
2223 C
2224 C------end of subroutine ddatrp------
2225 END
2226 SUBROUTINE DDAWTS(NEQ,IWT,RTOL,ATOL,Y,WT,RPAR,IPAR)
2227 C
2228 C***BEGIN PROLOGUE DDAWTS
2229 C***REFER TO DDASSL
2230 C***ROUTINES CALLED (NONE)
2231 C***DATE WRITTEN 830315 (YYMMDD)
2232 C***REVISION DATE 830315 (YYMMDD)
2233 C***END PROLOGUE DDAWTS
2234 C-----------------------------------------------------------------------
2235 C this subroutine sets the error weight vector
2236 C wt according to wt(i)=rtol(i)*abs(y(i))+atol(i),
2237 C i=1,-,n.
2238 C rtol and atol are scalars if iwt = 0,
2239 C and vectors if iwt = 1.
2240 C-----------------------------------------------------------------------
2241 C
2242 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2243 DIMENSION RTOL(1),ATOL(1),Y(1),WT(1)
2244 DIMENSION RPAR(1),IPAR(1)
2245 RTOLI=RTOL(1)
2246 ATOLI=ATOL(1)
2247 DO 20 I=1,NEQ
2248 IF (IWT .EQ.0) GO TO 10
2249 RTOLI=RTOL(I)
2250 ATOLI=ATOL(I)
2251 10 WT(I)=RTOLI*DABS(Y(I))+ATOLI
2252 20 CONTINUE
2253 RETURN
2254 C-----------end of subroutine ddawts-------------------------------------
2255 END
2256 DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
2257 C
2258 C FORMS THE DOT PRODUCT OF TWO VECTORS.
2259 C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2260 C JACK DONGARRA, LINPACK, 3/11/78.
2261 C
2262 DOUBLE PRECISION DX(1),DY(1),DTEMP
2263 INTEGER I,INCX,INCY,IX,IY,M,MP1,N
2264 C
2265 DDOT = 0.0D0
2266 DTEMP = 0.0D0
2267 IF(N.LE.0)RETURN
2268 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2269 C
2270 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2271 C NOT EQUAL TO 1
2272 C
2273 IX = 1
2274 IY = 1
2275 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2276 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2277 DO 10 I = 1,N
2278 DTEMP = DTEMP + DX(IX)*DY(IY)
2279 IX = IX + INCX
2280 IY = IY + INCY
2281 10 CONTINUE
2282 DDOT = DTEMP
2283 RETURN
2284 C
2285 C CODE FOR BOTH INCREMENTS EQUAL TO 1
2286 C
2287 C
2288 C CLEAN-UP LOOP
2289 C
2290 20 M = MOD(N,5)
2291 IF( M .EQ. 0 ) GO TO 40
2292 DO 30 I = 1,M
2293 DTEMP = DTEMP + DX(I)*DY(I)
2294 30 CONTINUE
2295 IF( N .LT. 5 ) GO TO 60
2296 40 MP1 = M + 1
2297 DO 50 I = MP1,N,5
2298 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
2299 * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
2300 50 CONTINUE
2301 60 DDOT = DTEMP
2302 RETURN
2303 END
2304 SUBROUTINE DGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB)
2305 INTEGER LDA,N,ML,MU,IPVT(1),JOB
2306 DOUBLE PRECISION ABD(LDA,1),B(1)
2307 C
2308 C DGBSL SOLVES THE DOUBLE PRECISION BAND SYSTEM
2309 C A * X = B OR TRANS(A) * X = B
2310 C USING THE FACTORS COMPUTED BY DGBCO OR DGBFA.
2311 C
2312 C ON ENTRY
2313 C
2314 C ABD DOUBLE PRECISION(LDA, N)
2315 C THE OUTPUT FROM DGBCO OR DGBFA.
2316 C
2317 C LDA INTEGER
2318 C THE LEADING DIMENSION OF THE ARRAY ABD .
2319 C
2320 C N INTEGER
2321 C THE ORDER OF THE ORIGINAL MATRIX.
2322 C
2323 C ML INTEGER
2324 C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
2325 C
2326 C MU INTEGER
2327 C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
2328 C
2329 C IPVT INTEGER(N)
2330 C THE PIVOT VECTOR FROM DGBCO OR DGBFA.
2331 C
2332 C B DOUBLE PRECISION(N)
2333 C THE RIGHT HAND SIDE VECTOR.
2334 C
2335 C JOB INTEGER
2336 C = 0 TO SOLVE A*X = B ,
2337 C = NONZERO TO SOLVE TRANS(A)*X = B , WHERE
2338 C TRANS(A) IS THE TRANSPOSE.
2339 C
2340 C ON RETURN
2341 C
2342 C B THE SOLUTION VECTOR X .
2343 C
2344 C ERROR CONDITION
2345 C
2346 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
2347 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY
2348 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
2349 C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE
2350 C CALLED CORRECTLY AND IF DGBCO HAS SET RCOND .GT. 0.0
2351 C OR DGBFA HAS SET INFO .EQ. 0 .
2352 C
2353 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
2354 C WITH P COLUMNS
2355 C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
2356 C IF (RCOND IS TOO SMALL) GO TO ...
2357 C DO 10 J = 1, P
2358 C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
2359 C 10 CONTINUE
2360 C
2361 C LINPACK. THIS VERSION DATED 08/14/78 .
2362 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
2363 C
2364 C SUBROUTINES AND FUNCTIONS
2365 C
2366 C BLAS DAXPY,DDOT
2367 C FORTRAN MIN0
2368 C
2369 C INTERNAL VARIABLES
2370 C
2371 DOUBLE PRECISION DDOT,T
2372 INTEGER K,KB,L,LA,LB,LM,M,NM1
2373 C
2374 M = MU + ML + 1
2375 NM1 = N - 1
2376 IF (JOB .NE. 0) GO TO 50
2377 C
2378 C JOB = 0 , SOLVE A * X = B
2379 C FIRST SOLVE L*Y = B
2380 C
2381 IF (ML .EQ. 0) GO TO 30
2382 IF (NM1 .LT. 1) GO TO 30
2383 DO 20 K = 1, NM1
2384 LM = MIN0(ML,N-K)
2385 L = IPVT(K)
2386 T = B(L)
2387 IF (L .EQ. K) GO TO 10
2388 B(L) = B(K)
2389 B(K) = T
2390 10 CONTINUE
2391 CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
2392 20 CONTINUE
2393 30 CONTINUE
2394 C
2395 C NOW SOLVE U*X = Y
2396 C
2397 DO 40 KB = 1, N
2398 K = N + 1 - KB
2399 B(K) = B(K)/ABD(M,K)
2400 LM = MIN0(K,M) - 1
2401 LA = M - LM
2402 LB = K - LM
2403 T = -B(K)
2404 CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
2405 40 CONTINUE
2406 GO TO 100
2407 50 CONTINUE
2408 C
2409 C JOB = NONZERO, SOLVE TRANS(A) * X = B
2410 C FIRST SOLVE TRANS(U)*Y = B
2411 C
2412 DO 60 K = 1, N
2413 LM = MIN0(K,M) - 1
2414 LA = M - LM
2415 LB = K - LM
2416 T = DDOT(LM,ABD(LA,K),1,B(LB),1)
2417 B(K) = (B(K) - T)/ABD(M,K)
2418 60 CONTINUE
2419 C
2420 C NOW SOLVE TRANS(L)*X = Y
2421 C
2422 IF (ML .EQ. 0) GO TO 90
2423 IF (NM1 .LT. 1) GO TO 90
2424 DO 80 KB = 1, NM1
2425 K = N - KB
2426 LM = MIN0(ML,N-K)
2427 B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
2428 L = IPVT(K)
2429 IF (L .EQ. K) GO TO 70
2430 T = B(L)
2431 B(L) = B(K)
2432 B(K) = T
2433 70 CONTINUE
2434 80 CONTINUE
2435 90 CONTINUE
2436 100 CONTINUE
2437 RETURN
2438 END
2439 SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
2440 INTEGER LDA,N,IPVT(1),INFO
2441 DOUBLE PRECISION A(LDA,1)
2442 C
2443 C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
2444 C
2445 C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
2446 C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
2447 C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
2448 C
2449 C ON ENTRY
2450 C
2451 C A DOUBLE PRECISION(LDA, N)
2452 C THE MATRIX TO BE FACTORED.
2453 C
2454 C LDA INTEGER
2455 C THE LEADING DIMENSION OF THE ARRAY A .
2456 C
2457 C N INTEGER
2458 C THE ORDER OF THE MATRIX A .
2459 C
2460 C ON RETURN
2461 C
2462 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
2463 C WHICH WERE USED TO OBTAIN IT.
2464 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
2465 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
2466 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
2467 C
2468 C IPVT INTEGER(N)
2469 C AN INTEGER VECTOR OF PIVOT INDICES.
2470 C
2471 C INFO INTEGER
2472 C = 0 NORMAL VALUE.
2473 C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
2474 C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
2475 C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
2476 C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE
2477 C INDICATION OF SINGULARITY.
2478 C
2479 C LINPACK. THIS VERSION DATED 08/14/78 .
2480 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
2481 C
2482 C SUBROUTINES AND FUNCTIONS
2483 C
2484 C BLAS DAXPY,DSCAL,IDAMAX
2485 C
2486 C INTERNAL VARIABLES
2487 C
2488 DOUBLE PRECISION T
2489 INTEGER IDAMAX,J,K,KP1,L,NM1
2490 C
2491 C
2492 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
2493 C
2494 INFO = 0
2495 NM1 = N - 1
2496 IF (NM1 .LT. 1) GO TO 70
2497 DO 60 K = 1, NM1
2498 KP1 = K + 1
2499 C
2500 C FIND L = PIVOT INDEX
2501 C
2502 L = IDAMAX(N-K+1,A(K,K),1) + K - 1
2503 IPVT(K) = L
2504 C
2505 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
2506 C
2507 IF (A(L,K) .EQ. 0.0D0) GO TO 40
2508 C
2509 C INTERCHANGE IF NECESSARY
2510 C
2511 IF (L .EQ. K) GO TO 10
2512 T = A(L,K)
2513 A(L,K) = A(K,K)
2514 A(K,K) = T
2515 10 CONTINUE
2516 C
2517 C COMPUTE MULTIPLIERS
2518 C
2519 T = -1.0D0/A(K,K)
2520 CALL DSCAL(N-K,T,A(K+1,K),1)
2521 C
2522 C ROW ELIMINATION WITH COLUMN INDEXING
2523 C
2524 DO 30 J = KP1, N
2525 T = A(L,J)
2526 IF (L .EQ. K) GO TO 20
2527 A(L,J) = A(K,J)
2528 A(K,J) = T
2529 20 CONTINUE
2530 CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
2531 30 CONTINUE
2532 GO TO 50
2533 40 CONTINUE
2534 INFO = K
2535 50 CONTINUE
2536 60 CONTINUE
2537 70 CONTINUE
2538 IPVT(N) = N
2539 IF (A(N,N) .EQ. 0.0D0) INFO = N
2540 RETURN
2541 END
2542 SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB)
2543 INTEGER LDA,N,IPVT(1),JOB
2544 DOUBLE PRECISION A(LDA,1),B(1)
2545 C
2546 C DGESL SOLVES THE DOUBLE PRECISION SYSTEM
2547 C A * X = B OR TRANS(A) * X = B
2548 C USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
2549 C
2550 C ON ENTRY
2551 C
2552 C A DOUBLE PRECISION(LDA, N)
2553 C THE OUTPUT FROM DGECO OR DGEFA.
2554 C
2555 C LDA INTEGER
2556 C THE LEADING DIMENSION OF THE ARRAY A .
2557 C
2558 C N INTEGER
2559 C THE ORDER OF THE MATRIX A .
2560 C
2561 C IPVT INTEGER(N)
2562 C THE PIVOT VECTOR FROM DGECO OR DGEFA.
2563 C
2564 C B DOUBLE PRECISION(N)
2565 C THE RIGHT HAND SIDE VECTOR.
2566 C
2567 C JOB INTEGER
2568 C = 0 TO SOLVE A*X = B ,
2569 C = NONZERO TO SOLVE TRANS(A)*X = B WHERE
2570 C TRANS(A) IS THE TRANSPOSE.
2571 C
2572 C ON RETURN
2573 C
2574 C B THE SOLUTION VECTOR X .
2575 C
2576 C ERROR CONDITION
2577 C
2578 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
2579 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY
2580 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
2581 C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE
2582 C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0
2583 C OR DGEFA HAS SET INFO .EQ. 0 .
2584 C
2585 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
2586 C WITH P COLUMNS
2587 C CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
2588 C IF (RCOND IS TOO SMALL) GO TO ...
2589 C DO 10 J = 1, P
2590 C CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
2591 C 10 CONTINUE
2592 C
2593 C LINPACK. THIS VERSION DATED 08/14/78 .
2594 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
2595 C
2596 C SUBROUTINES AND FUNCTIONS
2597 C
2598 C BLAS DAXPY,DDOT
2599 C
2600 C INTERNAL VARIABLES
2601 C
2602 DOUBLE PRECISION DDOT,T
2603 INTEGER K,KB,L,NM1
2604 C
2605 NM1 = N - 1
2606 IF (JOB .NE. 0) GO TO 50
2607 C
2608 C JOB = 0 , SOLVE A * X = B
2609 C FIRST SOLVE L*Y = B
2610 C
2611 IF (NM1 .LT. 1) GO TO 30
2612 DO 20 K = 1, NM1
2613 L = IPVT(K)
2614 T = B(L)
2615 IF (L .EQ. K) GO TO 10
2616 B(L) = B(K)
2617 B(K) = T
2618 10 CONTINUE
2619 CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
2620 20 CONTINUE
2621 30 CONTINUE
2622 C
2623 C NOW SOLVE U*X = Y
2624 C
2625 DO 40 KB = 1, N
2626 K = N + 1 - KB
2627 B(K) = B(K)/A(K,K)
2628 T = -B(K)
2629 CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
2630 40 CONTINUE
2631 GO TO 100
2632 50 CONTINUE
2633 C
2634 C JOB = NONZERO, SOLVE TRANS(A) * X = B
2635 C FIRST SOLVE TRANS(U)*Y = B
2636 C
2637 DO 60 K = 1, N
2638 T = DDOT(K-1,A(1,K),1,B(1),1)
2639 B(K) = (B(K) - T)/A(K,K)
2640 60 CONTINUE
2641 C
2642 C NOW SOLVE TRANS(L)*X = Y
2643 C
2644 IF (NM1 .LT. 1) GO TO 90
2645 DO 80 KB = 1, NM1
2646 K = N - KB
2647 B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
2648 L = IPVT(K)
2649 IF (L .EQ. K) GO TO 70
2650 T = B(L)
2651 B(L) = B(K)
2652 B(K) = T
2653 70 CONTINUE
2654 80 CONTINUE
2655 90 CONTINUE
2656 100 CONTINUE
2657 RETURN
2658 END
2659 SUBROUTINE DSCAL(N,DA,DX,INCX)
2660 C
2661 C SCALES A VECTOR BY A CONSTANT.
2662 C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
2663 C JACK DONGARRA, LINPACK, 3/11/78.
2664 C
2665 DOUBLE PRECISION DA,DX(1)
2666 INTEGER I,INCX,M,MP1,N,NINCX
2667 C
2668 IF(N.LE.0)RETURN
2669 IF(INCX.EQ.1)GO TO 20
2670 C
2671 C CODE FOR INCREMENT NOT EQUAL TO 1
2672 C
2673 NINCX = N*INCX
2674 DO 10 I = 1,NINCX,INCX
2675 DX(I) = DA*DX(I)
2676 10 CONTINUE
2677 RETURN
2678 C
2679 C CODE FOR INCREMENT EQUAL TO 1
2680 C
2681 C
2682 C CLEAN-UP LOOP
2683 C
2684 20 M = MOD(N,5)
2685 IF( M .EQ. 0 ) GO TO 40
2686 DO 30 I = 1,M
2687 DX(I) = DA*DX(I)
2688 30 CONTINUE
2689 IF( N .LT. 5 ) RETURN
2690 40 MP1 = M + 1
2691 DO 50 I = MP1,N,5
2692 DX(I) = DA*DX(I)
2693 DX(I + 1) = DA*DX(I + 1)
2694 DX(I + 2) = DA*DX(I + 2)
2695 DX(I + 3) = DA*DX(I + 3)
2696 DX(I + 4) = DA*DX(I + 4)
2697 50 CONTINUE
2698 RETURN
2699 END
2700 INTEGER FUNCTION IDAMAX(N,DX,INCX)
2701 C
2702 C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
2703 C JACK DONGARRA, LINPACK, 3/11/78.
2704 C
2705 DOUBLE PRECISION DX(1),DMAX
2706 INTEGER I,INCX,IX,N
2707 C
2708 IDAMAX = 0
2709 IF( N .LT. 1 ) RETURN
2710 IDAMAX = 1
2711 IF(N.EQ.1)RETURN
2712 IF(INCX.EQ.1)GO TO 20
2713 C
2714 C CODE FOR INCREMENT NOT EQUAL TO 1
2715 C
2716 IX = 1
2717 DMAX = DABS(DX(1))
2718 IX = IX + INCX
2719 DO 10 I = 2,N
2720 IF(DABS(DX(IX)).LE.DMAX) GO TO 5
2721 IDAMAX = I
2722 DMAX = DABS(DX(IX))
2723 5 IX = IX + INCX
2724 10 CONTINUE
2725 RETURN
2726 C
2727 C CODE FOR INCREMENT EQUAL TO 1
2728 C
2729 20 DMAX = DABS(DX(1))
2730 DO 30 I = 2,N
2731 IF(DABS(DX(I)).LE.DMAX) GO TO 30
2732 IDAMAX = I
2733 DMAX = DABS(DX(I))
2734 30 CONTINUE
2735 RETURN
2736 END
2737 SUBROUTINE MA30AD(NN, ICN, A, LICN, LENR, LENRL, IDISP, IP, IQ,
2738 * IRN, LIRN, LENC, IFIRST, LASTR, NEXTR, LASTC, NEXTC, IPTR, IPC, 0000000
2739 * U, IFLAG) 0000000
2740 c_270390
2741 EXTERNAL MA30$DATA
2742 c_270390
2743 C IF THE USER REQUIRES A MORE CONVENIENT DATA INTERFACE THEN THE MA28 0000000
2744 C PACKAGE SHOULD BE USED. THE MA28 SUBROUTINES CALL THE MA30 0000000
2745 C SUBROUTINES AFTER CHECKING THE USER'S INPUT DATA AND OPTIONALLY 0000000
2746 C USING MC23A/AD TO PERMUTE THE MATRIX TO BLOCK TRIANGULAR FORM. 0000000
2747 C THIS PACKAGE OF SUBROUTINES (MA30A/AD, MA30B/BD, MA30C/CD AND 0000000
2748 C MA30D/DD) PERFORMS OPERATIONS PERTINENT TO THE SOLUTION OF A 0000000
2749 C GENERAL SPARSE N BY N SYSTEM OF LINEAR EQUATIONS (I.E. SOLVE 0000000
2750 C AX=B). STRUCTUALLY SINGULAR MATRICES ARE PERMITTED INCLUDING 0000000
2751 C THOSE WITH ROW OR COLUMNS CONSISTING ENTIRELY OF ZEROS (I.E. 0000000
2752 C INCLUDING RECTANGULAR MATRICES). IT IS ASSUMED THAT THE 0000000
2753 C NON-ZEROS OF THE MATRIX A DO NOT DIFFER WIDELY IN SIZE. IF 0000000
2754 C NECESSARY A PRIOR CALL OF THE SCALING SUBROUTINE MC19A/AD MAY BE 0000000
2755 C MADE. 0000000
2756 C A DISCUSSION OF THE DESIGN OF THESE SUBROUTINES IS GIVEN BY DUFF AND 0000000
2757 C REID (ACM TRANS MATH SOFTWARE 5 PP 18-35,1979 (CSS 48)) WHILE 0000000
2758 C FULLER DETAILS OF THE IMPLEMENTATION ARE GIVEN IN DUFF (HARWELL 0000000
2759 C REPORT AERE-R 8730,1977). THE ADDITIONAL PIVOTING OPTION IN 0000000
2760 C MA30A/AD AND THE USE OF DROP TOLERANCES (SEE COMMON BLOCK 0000000
2761 C MA30I/ID) WERE ADDED TO THE PACKAGE AFTER JOINT WORK WITH REID, 0000000
2762 C SCHAUMBURG, WASNIEWSKI AND ZLATEV (DUFF, REID, SCHAUMBURG, 0000000
2763 C WASNIEWSKI AND ZLATEV, HARWELL REPORT CSS 135, 1983). 0000000
2764 C 0000000
2765 C MA30A/AD PERFORMS THE LU DECOMPOSITION OF THE DIAGONAL BLOCKS OF THE 0000000
2766 C PERMUTATION PAQ OF A SPARSE MATRIX A, WHERE INPUT PERMUTATIONS 0000000
2767 C P1 AND Q1 ARE USED TO DEFINE THE DIAGONAL BLOCKS. THERE MAY BE 0000000
2768 C NON-ZEROS IN THE OFF-DIAGONAL BLOCKS BUT THEY ARE UNAFFECTED BY 0000000
2769 C MA30A/AD. P AND P1 DIFFER ONLY WITHIN BLOCKS AS DO Q AND Q1. THE 0000000
2770 C PERMUTATIONS P1 AND Q1 MAY BE FOUND BY CALLING MC23A/AD OR THE 0000000
2771 C MATRIX MAY BE TREATED AS A SINGLE BLOCK BY USING P1=Q1=I. THE 0000000
2772 C MATRIX NON-ZEROS SHOULD BE HELD COMPACTLY BY ROWS, ALTHOUGH IT 0000000
2773 C SHOULD BE NOTED THAT THE USER CAN SUPPLY THE MATRIX BY COLUMNS 0000000
2774 C TO GET THE LU DECOMPOSITION OF A TRANSPOSE. 0000000
2775 C 0000000
2776 C THE PARAMETERS ARE... 0000000
2777 C THIS DESCRIPTION SHOULD ALSO BE CONSULTED FOR FURTHER INFORMATION ON 0000000
2778 C MOST OF THE PARAMETERS OF MA30B/BD AND MA30C/CD. 0000000
2779 C 0000000
2780 C N IS AN INTEGER VARIABLE WHICH MUST BE SET BY THE USER TO THE
2781 C ORDER OF THE MATRIX. BECAUSE OF THE USE OF INTEGER*2 ARRAYS IN
2782 C THE IBM VERSION, THE VALUE OF N SHOULD BE LESS THAN 32768. IT
2783 C IS NOT ALTERED BY MA30A/AD.
2784 C ICN IS AN INTEGER*2 ARRAY OF LENGTH LICN. POSITIONS IDISP(2) TO
2785 C LICN MUST BE SET BY THE USER TO CONTAIN THE COLUMN INDICES OF 0000000
2786 C THE NON-ZEROS IN THE DIAGONAL BLOCKS OF P1*A*Q1. THOSE BELONGING 0000000
2787 C TO A SINGLE ROW MUST BE CONTIGUOUS BUT THE ORDERING OF COLUMN 0000000
2788 C INDICES WITH EACH ROW IS UNIMPORTANT. THE NON-ZEROS OF ROW I 0000000
2789 C PRECEDE THOSE OF ROW I+1,I=1,...,N-1 AND NO WASTED SPACE IS 0000000
2790 C ALLOWED BETWEEN THE ROWS. ON OUTPUT THE COLUMN INDICES OF THE 0000000
2791 C LU DECOMPOSITION OF PAQ ARE HELD IN POSITIONS IDISP(1) TO 0000000
2792 C IDISP(2), THE ROWS ARE IN PIVOTAL ORDER, AND THE COLUMN INDICES 0000000
2793 C OF THE L PART OF EACH ROW ARE IN PIVOTAL ORDER AND PRECEDE THOSE 0000000
2794 C OF U. AGAIN THERE IS NO WASTED SPACE EITHER WITHIN A ROW OR 0000000
2795 C BETWEEN THE ROWS. ICN(1) TO ICN(IDISP(1)-1), ARE NEITHER 0000000
2796 C REQUIRED NOR ALTERED. IF MC23A/AD BEEN CALLED, THESE WILL HOLD 0000000
2797 C INFORMATION ABOUT THE OFF-DIAGONAL BLOCKS. 0000000
2798 C A IS A REAL/DOUBLE PRECISION ARRAY OF LENGTH LICN WHOSE ENTRIES 0000000
2799 C IDISP(2) TO LICN MUST BE SET BY THE USER TO THE VALUES OF THE 0000000
2800 C NON-ZERO ENTRIES OF THE MATRIX IN THE ORDER INDICATED BY ICN. 0000000
2801 C ON OUTPUT A WILL HOLD THE LU FACTORS OF THE MATRIX WHERE AGAIN 0000000
2802 C THE POSITION IN THE MATRIX IS DETERMINED BY THE CORRESPONDING 0000000
2803 C VALUES IN ICN. A(1) TO A(IDISP(1)-1) ARE NEITHER REQUIRED NOR 0000000
2804 C ALTERED. 0000000
2805 C LICN IS AN INTEGER VARIABLE WHICH MUST BE SET BY THE USER TO THE 0000000
2806 C LENGTH OF ARRAYS ICN AND A. IT MUST BE BIG ENOUGH FOR A AND ICN 0000000
2807 C TO HOLD ALL THE NON-ZEROS OF L AND U AND LEAVE SOME "ELBOW 0000000
2808 C ROOM". IT IS POSSIBLE TO CALCULATE A MINIMUM VALUE FOR LICN BY 0000000
2809 C A PRELIMINARY RUN OF MA30A/AD. THE ADEQUACY OF THE ELBOW ROOM 0000000
2810 C CAN BE JUDGED BY THE SIZE OF THE COMMON BLOCK VARIABLE ICNCP. IT 0000000
2811 C IS NOT ALTERED BY MA30A/AD. 0000000
2812 C LENR IS AN INTEGER*2 ARRAY OF LENGTH N. ON INPUT, LENR(I) SHOULD
2813 C EQUAL THE NUMBER OF NON-ZEROS IN ROW I, I=1,...,N OF THE 0000000
2814 C DIAGONAL BLOCKS OF P1*A*Q1. ON OUTPUT, LENR(I) WILL EQUAL THE 0000000
2815 C TOTAL NUMBER OF NON-ZEROS IN ROW I OF L AND ROW I OF U. 0000000
2816 C LENRL IS AN INTEGER*2 ARRAY OF LENGTH N. ON OUTPUT FROM MA30A/AD,
2817 C LENRL(I) WILL HOLD THE NUMBER OF NON-ZEROS IN ROW I OF L. 0000000
2818 C IDISP IS AN INTEGER ARRAY OF LENGTH 2. THE USER SHOULD SET IDISP(1) 0000000
2819 C TO BE THE FIRST AVAILABLE POSITION IN A/ICN FOR THE LU 0000000
2820 C DECOMPOSITION WHILE IDISP(2) IS SET TO THE POSITION IN A/ICN OF 0000000
2821 C THE FIRST NON-ZERO IN THE DIAGONAL BLOCKS OF P1*A*Q1. ON OUTPUT, 0000000
2822 C IDISP(1) WILL BE UNALTERED WHILE IDISP(2) WILL BE SET TO THE 0000000
2823 C POSITION IN A/ICN OF THE LAST NON-ZERO OF THE LU DECOMPOSITION. 0000000
2824 C IP IS AN INTEGER*2 ARRAY OF LENGTH N WHICH HOLDS A PERMUTATION OF
2825 C THE INTEGERS 1 TO N. ON INPUT TO MA30A/AD, THE ABSOLUTE VALUE OF 0000000
2826 C IP(I) MUST BE SET TO THE ROW OF A WHICH IS ROW I OF P1*A*Q1. A 0000000
2827 C NEGATIVE VALUE FOR IP(I) INDICATES THAT ROW I IS AT THE END OF A 0000000
2828 C DIAGONAL BLOCK. ON OUTPUT FROM MA30A/AD, IP(I) INDICATES THE ROW 0000000
2829 C OF A WHICH IS THE I TH ROW IN PAQ. IP(I) WILL STILL BE NEGATIVE 0000000
2830 C FOR THE LAST ROW OF EACH BLOCK (EXCEPT THE LAST). 0000000
2831 C IQ IS AN INTEGER*2 ARRAY OF LENGTH N WHICH AGAIN HOLDS A
2832 C PERMUTATION OF THE INTEGERS 1 TO N. ON INPUT TO MA30A/AD, IQ(J) 0000000
2833 C MUST BE SET TO THE COLUMN OF A WHICH IS COLUMN J OF P1*A*Q1. ON 0000000
2834 C OUTPUT FROM MA30A/AD, THE ABSOLUTE VALUE OF IQ(J) INDICATES THE 0000000
2835 C COLUMN OF A WHICH IS THE J TH IN PAQ. FOR ROWS, I SAY, IN WHICH 0000000
2836 C STRUCTURAL OR NUMERICAL SINGULARITY IS DETECTED IQ(I) IS 0000000
2837 C NEGATED. 0000000
2838 C IRN IS AN INTEGER*2 ARRAY OF LENGTH LIRN USED AS WORKSPACE BY
2839 C MA30A/AD. 0000000
2840 C LIRN IS AN INTEGER VARIABLE. IT SHOULD BE GREATER THAN THE 0000000
2841 C LARGEST NUMBER OF NON-ZEROS IN A DIAGONAL BLOCK OF P1*A*Q1 BUT 0000000
2842 C NEED NOT BE AS LARGE AS LICN. IT IS THE LENGTH OF ARRAY IRN AND 0000000
2843 C SHOULD BE LARGE ENOUGH TO HOLD THE ACTIVE PART OF ANY BLOCK, 0000000
2844 C PLUS SOME "ELBOW ROOM", THE A POSTERIORI ADEQUACY OF WHICH CAN 0000000
2845 C BE ESTIMATED BY EXAMINING THE SIZE OF COMMON BLOCK VARIABLE 0000000
2846 C IRNCP. 0000000
2847 C LENC,IFIRST,LASTR,NEXTR,LASTC,NEXTC ARE ALL INTEGER*2 ARRAYS OF
2848 C LENGTH N WHICH ARE USED AS WORKSPACE BY MA30A/AD. IF NSRCH IS 0000000
2849 C SET TO A VALUE LESS THAN OR EQUAL TO N, THEN ARRAYS LASTC AND 0000000
2850 C NEXTC ARE NOT REFERENCED BY MA30A/AD AND SO CAN BE DUMMIED IN 0000000
2851 C THE CALL TO MA30A/AD. 0000000
2852 C IPTR,IPC ARE INTEGER ARRAYS OF LENGTH N WHICH ARE USED AS WORKSPACE 0000000
2853 C BY MA30A/AD. 0000000
2854 C U IS A REAL/DOUBLE PRECISION VARIABLE WHICH SHOULD BE SET BY THE 0000000
2855 C USER TO A VALUE BETWEEN 0. AND 1.0. IF LESS THAN ZERO IT IS 0000000
2856 C RESET TO ZERO AND IF ITS VALUE IS 1.0 OR GREATER IT IS RESET TO 0000000
2857 C 0.9999 (0.999999999 IN D VERSION). IT DETERMINES THE BALANCE 0000000
2858 C BETWEEN PIVOTING FOR SPARSITY AND FOR STABILITY, VALUES NEAR 0000000
2859 C ZERO EMPHASIZING SPARSITY AND VALUES NEAR ONE EMPHASIZING 0000000
2860 C STABILITY. WE RECOMMEND U=0.1 AS A POSIBLE FIRST TRIAL VALUE. 0000000
2861 C THE STABILITY CAN BE JUDGED BY A LATER CALL TO MC24A/AD OR BY 0000000
2862 C SETTING LBIG TO .TRUE. 0000000
2863 C IFLAG IS AN INTEGER VARIABLE. IT WILL HAVE A NON-NEGATIVE VALUE IF 0000000
2864 C MA30A/AD IS SUCCESSFUL. NEGATIVE VALUES INDICATE ERROR 0000000
2865 C CONDITIONS WHILE POSITIVE VALUES INDICATE THAT THE MATRIX HAS 0000000
2866 C BEEN SUCCESSFULLY DECOMPOSED BUT IS SINGULAR. FOR EACH NON-ZERO 0000000
2867 C VALUE, AN APPROPRIATE MESSAGE IS OUTPUT ON UNIT LP. POSSIBLE 0000000
2868 C NON-ZERO VALUES FOR IFLAG ARE ... 0000000
2869 C 0000000
2870 C -1 THE MATRIX IS STRUCTUALLY SINGULAR WITH RANK GIVEN BY IRANK IN 0000000
2871 C COMMON BLOCK MA30F/FD. 0000000
2872 C +1 IF, HOWEVER, THE USER WANTS THE LU DECOMPOSITION OF A 0000000
2873 C STRUCTURALLY SINGULAR MATRIX AND SETS THE COMMON BLOCK VARIABLE 0000000
2874 C ABORT1 TO .FALSE., THEN, IN THE EVENT OF SINGULARITY AND A 0000000
2875 C SUCCESSFUL DECOMPOSITION, IFLAG IS RETURNED WITH THE VALUE +1 0000000
2876 C AND NO MESSAGE IS OUTPUT. 0000000
2877 C -2 THE MATRIX IS NUMERICALLY SINGULAR (IT MAY ALSO BE STRUCTUALLY 0000000
2878 C SINGULAR) WITH ESTIMATED RANK GIVEN BY IRANK IN COMMON BLOCK 0000000
2879 C MA30F/FD. 0000000
2880 C +2 THE USER CAN CHOOSE TO CONTINUE THE DECOMPOSITION EVEN WHEN A 0000000
2881 C ZERO PIVOT IS ENCOUNTERED BY SETTING COMMON BLOCK VARIABLE 0000000
2882 C ABORT2 TO .FALSE. IF A SINGULARITY IS ENCOUNTERED, IFLAG WILL 0000000
2883 C THEN RETURN WITH A VALUE OF +2, AND NO MESSAGE IS OUTPUT IF THE 0000000
2884 C DECOMPOSITION HAS BEEN COMPLETED SUCCESSFULLY. 0000000
2885 C -3 LIRN HAS NOT BEEN LARGE ENOUGH TO CONTINUE WITH THE 0000000
2886 C DECOMPOSITION. IF THE STAGE WAS ZERO THEN COMMON BLOCK VARIABLE 0000000
2887 C MINIRN GIVES THE LENGTH SUFFICIENT TO START THE DECOMPOSITION ON 0000000
2888 C THIS BLOCK. FOR A SUCCESSFUL DECOMPOSITION ON THIS BLOCK THE 0000000
2889 C USER SHOULD MAKE LIRN SLIGHTLY (SAY ABOUT N/2) GREATER THAN THIS 0000000
2890 C VALUE. 0000000
2891 C -4 LICN NOT LARGE ENOUGH TO CONTINUE WITH THE DECOMPOSITION. 0000000
2892 C -5 THE DECOMPOSITION HAS BEEN COMPLETED BUT SOME OF THE LU FACTORS 0000000
2893 C HAVE BEEN DISCARDED TO CREATE ENOUGH ROOM IN A/ICN TO CONTINUE 0000000
2894 C THE DECOMPOSITION. THE VARIABLE MINICN IN COMMON BLOCK MA30F/FD 0000000
2895 C THEN GIVES THE SIZE THAT LICN SHOULD BE TO ENABLE THE 0000000
2896 C FACTORIZATION TO BE SUCCESSFUL. IF THE USER SETS COMMON BLOCK 0000000
2897 C VARIABLE ABORT3 TO .TRUE., THEN THE SUBROUTINE WILL EXIT 0000000
2898 C IMMEDIATELY INSTEAD OF DESTROYING ANY FACTORS AND CONTINUING. 0000000
2899 C -6 BOTH LICN AND LIRN ARE TOO SMALL. TERMINATION HAS BEEN CAUSED BY 0000000
2900 C LACK OF SPACE IN IRN (SEE ERROR IFLAG= -3), BUT ALREADY SOME OF 0000000
2901 C THE LU FACTORS IN A/ICN HAVE BEEN LOST (SEE ERROR IFLAG= -5). 0000000
2902 C MINICN GIVES THE MINIMUM AMOUNT OF SPACE REQUIRED IN A/ICN FOR 0000000
2903 C DECOMPOSITION UP TO THIS POINT. 0000000
2904 C 0000000
2905 DOUBLE PRECISION A(LICN), U, AU, UMAX, AMAX, ZERO, PIVRAT, PIVR,
2906 * TOL, BIG, ANEW, AANEW, SCALE
2907 INTEGER IPTR(NN), PIVOT, PIVEND, DISPC, OLDPIV, OLDEND, PIVROW, 0000000
2908 * ROWI, IPC(NN), IDISP(2) 0000000
2909 Change
2910 C INTEGER*2 ICN(LICN), LENR(NN), LENRL(NN), IP(NN), IQ(NN),
2911 C * LENC(NN), IRN(LIRN), IFIRST(NN), LASTR(NN), NEXTR(NN), 0000000
2912 C * LASTC(NN), NEXTC(NN) 0000000
2913 INTEGER ICN(LICN), LENR(NN), LENRL(NN), IP(NN), IQ(NN),
2914 * LENC(NN), IRN(LIRN), IFIRST(NN), LASTR(NN), NEXTR(NN), 0000000
2915 * LASTC(NN), NEXTC(NN) 0000000
2916 Change
2917 LOGICAL ABORT1, ABORT2, ABORT3, LBIG 0000000
2918 C FOR COMMENTS OF COMMON BLOCK VARIABLES SEE BLOCK DATA SUBPROGRAM.
2919 COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
2920 COMMON /MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN
2921 COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
2922 C 0000000
2923 DATA UMAX/.999999999D0/
2924 DATA ZERO /0.0D0/
2925 MSRCH = NSRCH 0000000
2926 NDROP = 0 0000000
2927 MINIRN = 0 0000000
2928 MINICN = IDISP(1) - 1 0000000
2929 MOREI = 0 0000000
2930 IRANK = NN 0000000
2931 IRNCP = 0 0000000
2932 ICNCP = 0 0000000
2933 IFLAG = 0 0000000
2934 C RESET U IF NECESSARY. 0000000
2935 U = DMIN1(U,UMAX)
2936 C IBEG IS THE POSITION OF THE NEXT PIVOT ROW AFTER ELIMINATION STEP 0000000
2937 C USING IT. 0000000
2938 U = DMAX1(U,ZERO)
2939 IBEG = IDISP(1) 0000000
2940 C IACTIV IS THE POSITION OF THE FIRST ENTRY IN THE ACTIVE PART OF A/ICN.0000000
2941 IACTIV = IDISP(2) 0000000
2942 C NZROW IS CURRENT NUMBER OF NON-ZEROS IN ACTIVE AND UNPROCESSED PART 0000000
2943 C OF ROW FILE ICN. 0000000
2944 NZROW = LICN - IACTIV + 1 0000000
2945 MINICN = NZROW + MINICN 0000000
2946 C 0000000
2947 C COUNT THE NUMBER OF DIAGONAL BLOCKS AND SET UP POINTERS TO THE 0000000
2948 C BEGINNINGS OF THE ROWS. 0000000
2949 C NUM IS THE NUMBER OF DIAGONAL BLOCKS. 0000000
2950 NUM = 1 0000000
2951 IPTR(1) = IACTIV 0000000
2952 IF (NN.EQ.1) GO TO 20 0000000
2953 NNM1 = NN - 1 0000000
2954 DO 10 I=1,NNM1 0000000
2955 IF (IP(I).LT.0) NUM = NUM + 1 0000000
2956 IPTR(I+1) = IPTR(I) + LENR(I) 0000000
2957 10 CONTINUE 0000000
2958 C ILAST IS THE LAST ROW IN THE PREVIOUS BLOCK. 0000000
2959 20 ILAST = 0 0000000
2960 C 0000000
2961 C *********************************************** 0000000
2962 C **** LU DECOMPOSITION OF BLOCK NBLOCK **** 0000000
2963 C *********************************************** 0000000
2964 C 0000000
2965 C EACH PASS THROUGH THIS LOOP PERFORMS LU DECOMPOSITION ON ONE 0000000
2966 C OF THE DIAGONAL BLOCKS. 0000000
2967 DO 1000 NBLOCK=1,NUM 0000000
2968 ISTART = ILAST + 1 0000000
2969 DO 30 IROWS=ISTART,NN 0000000
2970 IF (IP(IROWS).LT.0) GO TO 40 0000000
2971 30 CONTINUE 0000000
2972 IROWS = NN 0000000
2973 40 ILAST = IROWS 0000000
2974 C N IS THE NUMBER OF ROWS IN THE CURRENT BLOCK. 0000000
2975 C ISTART IS THE INDEX OF THE FIRST ROW IN THE CURRENT BLOCK. 0000000
2976 C ILAST IS THE INDEX OF THE LAST ROW IN THE CURRENT BLOCK. 0000000
2977 C IACTIV IS THE POSITION OF THE FIRST ENTRY IN THE BLOCK. 0000000
2978 C ITOP IS THE POSITION OF THE LAST ENTRY IN THE BLOCK. 0000000
2979 N = ILAST - ISTART + 1 0000000
2980 IF (N.NE.1) GO TO 90 0000000
2981 C 0000000
2982 C CODE FOR DEALING WITH 1X1 BLOCK. 0000000
2983 LENRL(ILAST) = 0 0000000
2984 ISING = ISTART 0000000
2985 IF (LENR(ILAST).NE.0) GO TO 50 0000000
2986 C BLOCK IS STRUCTURALLY SINGULAR. 0000000
2987 IRANK = IRANK - 1 0000000
2988 ISING = -ISING 0000000
2989 IF (IFLAG.NE.2 .AND. IFLAG.NE.-5) IFLAG = 1 0000000
2990 IF (.NOT.ABORT1) GO TO 80 0000000
2991 IDISP(2) = IACTIV 0000000
2992 IFLAG = -1 0000000
2993 IF (LP.NE.0) WRITE (LP,99999) 0000000
2994 C RETURN 0000000
2995 GO TO 1120 0000000
2996 50 SCALE = DABS(A(IACTIV))
2997 IF (SCALE.EQ.ZERO) GO TO 60 0000000
2998 IF (LBIG) BIG = DMAX1(BIG,SCALE)
2999 GO TO 70 0000000
3000 60 ISING = -ISING 0000000
3001 IRANK = IRANK - 1 0000000
3002 IPTR(ILAST) = 0 0000000
3003 IF (IFLAG.NE.-5) IFLAG = 2 0000000
3004 IF (.NOT.ABORT2) GO TO 70 0000000
3005 IDISP(2) = IACTIV 0000000
3006 IFLAG = -2 0000000
3007 IF (LP.NE.0) WRITE (LP,99998) 0000000
3008 GO TO 1120 0000000
3009 70 A(IBEG) = A(IACTIV) 0000000
3010 ICN(IBEG) = ICN(IACTIV) 0000000
3011 IACTIV = IACTIV + 1 0000000
3012 IPTR(ISTART) = 0 0000000
3013 IBEG = IBEG + 1 0000000
3014 NZROW = NZROW - 1 0000000
3015 80 LASTR(ISTART) = ISTART 0000000
3016 IPC(ISTART) = -ISING 0000000
3017 GO TO 1000 0000000
3018 C 0000000
3019 C NON-TRIVIAL BLOCK. 0000000
3020 90 ITOP = LICN 0000000
3021 IF (ILAST.NE.NN) ITOP = IPTR(ILAST+1) - 1 0000000
3022 C 0000000
3023 C SET UP COLUMN ORIENTED STORAGE. 0000000
3024 DO 100 I=ISTART,ILAST 0000000
3025 LENRL(I) = 0 0000000
3026 LENC(I) = 0 0000000
3027 100 CONTINUE 0000000
3028 IF (ITOP-IACTIV.LT.LIRN) GO TO 110 0000000
3029 MINIRN = ITOP - IACTIV + 1 0000000
3030 PIVOT = ISTART - 1 0000000
3031 GO TO 1100 0000000
3032 C 0000000
3033 C CALCULATE COLUMN COUNTS. 0000000
3034 110 DO 120 II=IACTIV,ITOP 0000000
3035 I = ICN(II) 0000000
3036 LENC(I) = LENC(I) + 1 0000000
3037 120 CONTINUE 0000000
3038 C SET UP COLUMN POINTERS SO THAT IPC(J) POINTS TO POSITION AFTER END 0000000
3039 C OF COLUMN J IN COLUMN FILE. 0000000
3040 IPC(ILAST) = LIRN + 1 0000000
3041 J1 = ISTART + 1 0000000
3042 DO 130 JJ=J1,ILAST 0000000
3043 J = ILAST - JJ + J1 - 1 0000000
3044 IPC(J) = IPC(J+1) - LENC(J+1) 0000000
3045 130 CONTINUE 0000000
3046 DO 150 INDROW=ISTART,ILAST 0000000
3047 J1 = IPTR(INDROW) 0000000
3048 J2 = J1 + LENR(INDROW) - 1 0000000
3049 IF (J1.GT.J2) GO TO 150 0000000
3050 DO 140 JJ=J1,J2 0000000
3051 J = ICN(JJ) 0000000
3052 IPOS = IPC(J) - 1 0000000
3053 IRN(IPOS) = INDROW 0000000
3054 IPC(J) = IPOS 0000000
3055 140 CONTINUE 0000000
3056 150 CONTINUE 0000000
3057 C DISPC IS THE LOWEST INDEXED ACTIVE LOCATION IN THE COLUMN FILE. 0000000
3058 DISPC = IPC(ISTART) 0000000
3059 NZCOL = LIRN - DISPC + 1 0000000
3060 MINIRN = MAX0(NZCOL,MINIRN) 0000000
3061 NZMIN = 1 0000000
3062 C 0000000
3063 C INITIALIZE ARRAY IFIRST. IFIRST(I) = +/- K INDICATES THAT ROW/COL 0000000
3064 C K HAS I NON-ZEROS. IF IFIRST(I) = 0, THERE IS NO ROW OR COLUMN 0000000
3065 C WITH I NON ZEROS. 0000000
3066 DO 160 I=1,N 0000000
3067 IFIRST(I) = 0 0000000
3068 160 CONTINUE 0000000
3069 C 0000000
3070 C COMPUTE ORDERING OF ROW AND COLUMN COUNTS. 0000000
3071 C FIRST RUN THROUGH COLUMNS (FROM COLUMN N TO COLUMN 1). 0000000
3072 DO 180 JJ=ISTART,ILAST 0000000
3073 J = ILAST - JJ + ISTART 0000000
3074 NZ = LENC(J) 0000000
3075 IF (NZ.NE.0) GO TO 170 0000000
3076 IPC(J) = 0 0000000
3077 GO TO 180 0000000
3078 170 IF (NSRCH.LE.NN) GO TO 180 0000000
3079 ISW = IFIRST(NZ) 0000000
3080 IFIRST(NZ) = -J 0000000
3081 LASTC(J) = 0 0000000
3082 NEXTC(J) = -ISW 0000000
3083 ISW1 = IABS(ISW) 0000000
3084 IF (ISW.NE.0) LASTC(ISW1) = J 0000000
3085 180 CONTINUE 0000000
3086 C NOW RUN THROUGH ROWS (AGAIN FROM N TO 1). 0000000
3087 DO 210 II=ISTART,ILAST 0000000
3088 I = ILAST - II + ISTART 0000000
3089 NZ = LENR(I) 0000000
3090 IF (NZ.NE.0) GO TO 190 0000000
3091 IPTR(I) = 0 0000000
3092 LASTR(I) = 0 0000000
3093 GO TO 210 0000000
3094 190 ISW = IFIRST(NZ) 0000000
3095 IFIRST(NZ) = I 0000000
3096 IF (ISW.GT.0) GO TO 200 0000000
3097 NEXTR(I) = 0 0000000
3098 LASTR(I) = ISW 0000000
3099 GO TO 210 0000000
3100 200 NEXTR(I) = ISW 0000000
3101 LASTR(I) = LASTR(ISW) 0000000
3102 LASTR(ISW) = I 0000000
3103 210 CONTINUE 0000000
3104 C 0000000
3105 C 0000000
3106 C ********************************************** 0000000
3107 C **** START OF MAIN ELIMINATION LOOP **** 0000000
3108 C ********************************************** 0000000
3109 DO 980 PIVOT=ISTART,ILAST 0000000
3110 C 0000000
3111 C FIRST FIND THE PIVOT USING MARKOWITZ CRITERION WITH STABILITY 0000000
3112 C CONTROL. 0000000
3113 C JCOST IS THE MARKOWITZ COST OF THE BEST PIVOT SO FAR,.. THIS 0000000
3114 C PIVOT IS IN ROW IPIV AND COLUMN JPIV. 0000000
3115 NZ2 = NZMIN 0000000
3116 JCOST = N*N 0000000
3117 C 0000000
3118 C EXAMINE ROWS/COLUMNS IN ORDER OF ASCENDING COUNT. 0000000
3119 DO 340 L=1,2 0000000
3120 PIVRAT = ZERO 0000000
3121 ISRCH = 1 0000000
3122 LL = L 0000000
3123 C A PASS WITH L EQUAL TO 2 IS ONLY PERFORMED IN THE CASE OF SINGULARITY.0000000
3124 DO 330 NZ=NZ2,N 0000000
3125 IF (JCOST.LE.(NZ-1)**2) GO TO 420 0000000
3126 IJFIR = IFIRST(NZ) 0000000
3127 IF (IJFIR) 230, 220, 240 0000000
3128 220 IF (LL.EQ.1) NZMIN = NZ + 1 0000000
3129 GO TO 330 0000000
3130 230 LL = 2 0000000
3131 IJFIR = -IJFIR 0000000
3132 GO TO 290 0000000
3133 240 LL = 2 0000000
3134 C SCAN ROWS WITH NZ NON-ZEROS. 0000000
3135 DO 270 IDUMMY=1,N 0000000
3136 IF (JCOST.LE.(NZ-1)**2) GO TO 420 0000000
3137 IF (ISRCH.GT.MSRCH) GO TO 420 0000000
3138 IF (IJFIR.EQ.0) GO TO 280 0000000
3139 C ROW IJFIR IS NOW EXAMINED. 0000000
3140 I = IJFIR 0000000
3141 IJFIR = NEXTR(I) 0000000
3142 C FIRST CALCULATE MULTIPLIER THRESHOLD LEVEL. 0000000
3143 AMAX = ZERO 0000000
3144 J1 = IPTR(I) + LENRL(I) 0000000
3145 J2 = IPTR(I) + LENR(I) - 1 0000000
3146 DO 250 JJ=J1,J2 0000000
3147 AMAX = DMAX1(AMAX,DABS(A(JJ)))
3148 250 CONTINUE 0000000
3149 AU = AMAX*U 0000000
3150 ISRCH = ISRCH + 1 0000000
3151 C SCAN ROW FOR POSSIBLE PIVOTS 0000000
3152 DO 260 JJ=J1,J2 0000000
3153 IF (DABS(A(JJ)).LE.AU .AND. L.EQ.1) GO TO 260
3154 J = ICN(JJ) 0000000
3155 KCOST = (NZ-1)*(LENC(J)-1) 0000000
3156 IF (KCOST.GT.JCOST) GO TO 260 0000000
3157 PIVR = ZERO 0000000
3158 IF (AMAX.NE.ZERO) PIVR = DABS(A(JJ))/AMAX
3159 IF (KCOST.EQ.JCOST .AND. (PIVR.LE.PIVRAT .OR. 0000000
3160 * NSRCH.GT.NN+1)) GO TO 260 0000000
3161 C BEST PIVOT SO FAR IS FOUND. 0000000
3162 JCOST = KCOST 0000000
3163 IJPOS = JJ 0000000
3164 IPIV = I 0000000
3165 JPIV = J 0000000
3166 IF (MSRCH.GT.NN+1 .AND. JCOST.LE.(NZ-1)**2) GO TO 420 0000000
3167 PIVRAT = PIVR 0000000
3168 260 CONTINUE 0000000
3169 270 CONTINUE 0000000
3170 C 0000000
3171 C COLUMNS WITH NZ NON-ZEROS NOW EXAMINED. 0000000
3172 280 IJFIR = IFIRST(NZ) 0000000
3173 IJFIR = -LASTR(IJFIR) 0000000
3174 290 IF (JCOST.LE.NZ*(NZ-1)) GO TO 420 0000000
3175 IF (MSRCH.LE.NN) GO TO 330 0000000
3176 DO 320 IDUMMY=1,N 0000000
3177 IF (IJFIR.EQ.0) GO TO 330 0000000
3178 J = IJFIR 0000000
3179 IJFIR = NEXTC(IJFIR) 0000000
3180 I1 = IPC(J) 0000000
3181 I2 = I1 + NZ - 1 0000000
3182 C SCAN COLUMN J. 0000000
3183 DO 310 II=I1,I2 0000000
3184 I = IRN(II) 0000000
3185 KCOST = (NZ-1)*(LENR(I)-LENRL(I)-1) 0000000
3186 IF (KCOST.GE.JCOST) GO TO 310 0000000
3187 C PIVOT HAS BEST MARKOWITZ COUNT SO FAR ... NOW CHECK ITS 0000000
3188 C SUITABILITY ON NUMERIC GROUNDS BY EXAMINING THE OTHER NON-ZEROS 0000000
3189 C IN ITS ROW. 0000000
3190 J1 = IPTR(I) + LENRL(I) 0000000
3191 J2 = IPTR(I) + LENR(I) - 1 0000000
3192 C WE NEED A STABILITY CHECK ON SINGLETON COLUMNS BECAUSE OF POSSIBLE 0000000
3193 C PROBLEMS WITH UNDERDETERMINED SYSTEMS. 0000000
3194 AMAX = ZERO 0000000
3195 DO 300 JJ=J1,J2 0000000
3196 AMAX = DMAX1(AMAX,DABS(A(JJ)))
3197 IF (ICN(JJ).EQ.J) JPOS = JJ 0000000
3198 300 CONTINUE 0000000
3199 IF (DABS(A(JPOS)).LE.AMAX*U .AND. L.EQ.1) GO TO 310
3200 JCOST = KCOST 0000000
3201 IPIV = I 0000000
3202 JPIV = J 0000000
3203 IJPOS = JPOS 0000000
3204 IF (AMAX.NE.ZERO) PIVRAT = DABS(A(JPOS))/AMAX
3205 IF (JCOST.LE.NZ*(NZ-1)) GO TO 420 0000000
3206 310 CONTINUE 0000000
3207 C 0000000
3208 320 CONTINUE 0000000
3209 C 0000000
3210 330 CONTINUE 0000000
3211 C IN THE EVENT OF SINGULARITY, WE MUST MAKE SURE ALL ROWS AND COLUMNS 0000000
3212 C ARE TESTED. 0000000
3213 MSRCH = N 0000000
3214 C 0000000
3215 C MATRIX IS NUMERICALLY OR STRUCTURALLY SINGULAR ... WHICH IT IS WILL 0000000
3216 C BE DIAGNOSED LATER. 0000000
3217 IRANK = IRANK - 1 0000000
3218 340 CONTINUE 0000000
3219 C ASSIGN REST OF ROWS AND COLUMNS TO ORDERING ARRAY. 0000000
3220 C MATRIX IS STRUCTURALLY SINGULAR. 0000000
3221 IF (IFLAG.NE.2 .AND. IFLAG.NE.-5) IFLAG = 1 0000000
3222 IRANK = IRANK - ILAST + PIVOT + 1 0000000
3223 IF (.NOT.ABORT1) GO TO 350 0000000
3224 IDISP(2) = IACTIV 0000000
3225 IFLAG = -1 0000000
3226 IF (LP.NE.0) WRITE (LP,99999) 0000000
3227 GO TO 1120 0000000
3228 350 K = PIVOT - 1 0000000
3229 DO 390 I=ISTART,ILAST 0000000
3230 IF (LASTR(I).NE.0) GO TO 390 0000000
3231 K = K + 1 0000000
3232 LASTR(I) = K 0000000
3233 IF (LENRL(I).EQ.0) GO TO 380 0000000
3234 MINICN = MAX0(MINICN,NZROW+IBEG-1+MOREI+LENRL(I)) 0000000
3235 IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 360 0000000
3236 CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
3237 C CHECK NOW TO SEE IF MA30D/DD HAS CREATED ENOUGH AVAILABLE SPACE. 0000000
3238 IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 360 0000000
3239 C CREATE MORE SPACE BY DESTROYING PREVIOUSLY CREATED LU FACTORS. 0000000
3240 MOREI = MOREI + IBEG - IDISP(1) 0000000
3241 IBEG = IDISP(1) 0000000
3242 IF (LP.NE.0) WRITE (LP,99997) 0000000
3243 IFLAG = -5 0000000
3244 IF (ABORT3) GO TO 1090 0000000
3245 360 J1 = IPTR(I) 0000000
3246 J2 = J1 + LENRL(I) - 1 0000000
3247 IPTR(I) = 0 0000000
3248 DO 370 JJ=J1,J2 0000000
3249 A(IBEG) = A(JJ) 0000000
3250 ICN(IBEG) = ICN(JJ) 0000000
3251 ICN(JJ) = 0 0000000
3252 IBEG = IBEG + 1 0000000
3253 370 CONTINUE 0000000
3254 NZROW = NZROW - LENRL(I) 0000000
3255 380 IF (K.EQ.ILAST) GO TO 400 0000000
3256 390 CONTINUE 0000000
3257 400 K = PIVOT - 1 0000000
3258 DO 410 I=ISTART,ILAST 0000000
3259 IF (IPC(I).NE.0) GO TO 410 0000000
3260 K = K + 1 0000000
3261 IPC(I) = K 0000000
3262 IF (K.EQ.ILAST) GO TO 990 0000000
3263 410 CONTINUE 0000000
3264 C 0000000
3265 C THE PIVOT HAS NOW BEEN FOUND IN POSITION (IPIV,JPIV) IN LOCATION 0000000
3266 C IJPOS IN ROW FILE. 0000000
3267 C UPDATE COLUMN AND ROW ORDERING ARRAYS TO CORRESPOND WITH REMOVAL 0000000
3268 C OF THE ACTIVE PART OF THE MATRIX. 0000000
3269 420 ISING = PIVOT 0000000
3270 IF (A(IJPOS).NE.ZERO) GO TO 430 0000000
3271 C NUMERICAL SINGULARITY IS RECORDED HERE. 0000000
3272 ISING = -ISING 0000000
3273 IF (IFLAG.NE.-5) IFLAG = 2 0000000
3274 IF (.NOT.ABORT2) GO TO 430 0000000
3275 IDISP(2) = IACTIV 0000000
3276 IFLAG = -2 0000000
3277 IF (LP.NE.0) WRITE (LP,99998) 0000000
3278 GO TO 1120 0000000
3279 430 OLDPIV = IPTR(IPIV) + LENRL(IPIV) 0000000
3280 OLDEND = IPTR(IPIV) + LENR(IPIV) - 1 0000000
3281 C CHANGES TO COLUMN ORDERING. 0000000
3282 IF (NSRCH.LE.NN) GO TO 460 0000000
3283 DO 450 JJ=OLDPIV,OLDEND 0000000
3284 J = ICN(JJ) 0000000
3285 LC = LASTC(J) 0000000
3286 NC = NEXTC(J) 0000000
3287 IF (NC.NE.0) LASTC(NC) = LC 0000000
3288 IF (LC.EQ.0) GO TO 440 0000000
3289 NEXTC(LC) = NC 0000000
3290 GO TO 450 0000000
3291 440 NZ = LENC(J) 0000000
3292 ISW = IFIRST(NZ) 0000000
3293 IF (ISW.GT.0) LASTR(ISW) = -NC 0000000
3294 IF (ISW.LT.0) IFIRST(NZ) = -NC 0000000
3295 450 CONTINUE 0000000
3296 C CHANGES TO ROW ORDERING. 0000000
3297 460 I1 = IPC(JPIV) 0000000
3298 I2 = I1 + LENC(JPIV) - 1 0000000
3299 DO 480 II=I1,I2 0000000
3300 I = IRN(II) 0000000
3301 LR = LASTR(I) 0000000
3302 NR = NEXTR(I) 0000000
3303 IF (NR.NE.0) LASTR(NR) = LR 0000000
3304 IF (LR.LE.0) GO TO 470 0000000
3305 NEXTR(LR) = NR 0000000
3306 GO TO 480 0000000
3307 470 NZ = LENR(I) - LENRL(I) 0000000
3308 IF (NR.NE.0) IFIRST(NZ) = NR 0000000
3309 IF (NR.EQ.0) IFIRST(NZ) = LR 0000000
3310 480 CONTINUE 0000000
3311 C 0000000
3312 C MOVE PIVOT TO POSITION LENRL+1 IN PIVOT ROW AND MOVE PIVOT ROW 0000000
3313 C TO THE BEGINNING OF THE AVAILABLE STORAGE. 0000000
3314 C THE L PART AND THE PIVOT IN THE OLD COPY OF THE PIVOT ROW IS 0000000
3315 C NULLIFIED WHILE, IN THE STRICTLY UPPER TRIANGULAR PART, THE 0000000
3316 C COLUMN INDICES, J SAY, ARE OVERWRITTEN BY THE CORRESPONDING 0000000
3317 C ENTRY OF IQ (IQ(J)) AND IQ(J) IS SET TO THE NEGATIVE OF THE 0000000
3318 C DISPLACEMENT OF THE COLUMN INDEX FROM THE PIVOT ENTRY. 0000000
3319 IF (OLDPIV.EQ.IJPOS) GO TO 490 0000000
3320 AU = A(OLDPIV) 0000000
3321 A(OLDPIV) = A(IJPOS) 0000000
3322 A(IJPOS) = AU 0000000
3323 ICN(IJPOS) = ICN(OLDPIV) 0000000
3324 ICN(OLDPIV) = JPIV 0000000
3325 C CHECK TO SEE IF THERE IS SPACE IMMEDIATELY AVAILABLE IN A/ICN TO 0000000
3326 C HOLD NEW COPY OF PIVOT ROW. 0000000
3327 490 MINICN = MAX0(MINICN,NZROW+IBEG-1+MOREI+LENR(IPIV)) 0000000
3328 IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500 0000000
3329 CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
3330 OLDPIV = IPTR(IPIV) + LENRL(IPIV) 0000000
3331 OLDEND = IPTR(IPIV) + LENR(IPIV) - 1 0000000
3332 C CHECK NOW TO SEE IF MA30D/DD HAS CREATED ENOUGH AVAILABLE SPACE. 0000000
3333 IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500 0000000
3334 C CREATE MORE SPACE BY DESTROYING PREVIOUSLY CREATED LU FACTORS. 0000000
3335 MOREI = MOREI + IBEG - IDISP(1) 0000000
3336 IBEG = IDISP(1) 0000000
3337 IF (LP.NE.0) WRITE (LP,99997) 0000000
3338 IFLAG = -5 0000000
3339 IF (ABORT3) GO TO 1090 0000000
3340 IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500 0000000
3341 C THERE IS STILL NOT ENOUGH ROOM IN A/ICN. 0000000
3342 IFLAG = -4 0000000
3343 GO TO 1090 0000000
3344 C COPY PIVOT ROW AND SET UP IQ ARRAY. 0000000
3345 500 IJPOS = 0 0000000
3346 J1 = IPTR(IPIV) 0000000
3347 C 0000000
3348 DO 530 JJ=J1,OLDEND 0000000
3349 A(IBEG) = A(JJ) 0000000
3350 ICN(IBEG) = ICN(JJ) 0000000
3351 IF (IJPOS.NE.0) GO TO 510 0000000
3352 IF (ICN(JJ).EQ.JPIV) IJPOS = IBEG 0000000
3353 ICN(JJ) = 0 0000000
3354 GO TO 520 0000000
3355 510 K = IBEG - IJPOS 0000000
3356 J = ICN(JJ) 0000000
3357 ICN(JJ) = IQ(J) 0000000
3358 IQ(J) = -K 0000000
3359 520 IBEG = IBEG + 1 0000000
3360 530 CONTINUE 0000000
3361 C 0000000
3362 IJP1 = IJPOS + 1 0000000
3363 PIVEND = IBEG - 1 0000000
3364 LENPIV = PIVEND - IJPOS 0000000
3365 NZROW = NZROW - LENRL(IPIV) - 1 0000000
3366 IPTR(IPIV) = OLDPIV + 1 0000000
3367 IF (LENPIV.EQ.0) IPTR(IPIV) = 0 0000000
3368 C 0000000
3369 C REMOVE PIVOT ROW (INCLUDING PIVOT) FROM COLUMN ORIENTED FILE. 0000000
3370 DO 560 JJ=IJPOS,PIVEND 0000000
3371 J = ICN(JJ) 0000000
3372 I1 = IPC(J) 0000000
3373 LENC(J) = LENC(J) - 1 0000000
3374 C I2 IS LAST POSITION IN NEW COLUMN. 0000000
3375 I2 = IPC(J) + LENC(J) - 1 0000000
3376 IF (I2.LT.I1) GO TO 550 0000000
3377 DO 540 II=I1,I2 0000000
3378 IF (IRN(II).NE.IPIV) GO TO 540 0000000
3379 IRN(II) = IRN(I2+1) 0000000
3380 GO TO 550 0000000
3381 540 CONTINUE 0000000
3382 550 IRN(I2+1) = 0 0000000
3383 560 CONTINUE 0000000
3384 NZCOL = NZCOL - LENPIV - 1 0000000
3385 C 0000000
3386 C GO DOWN THE PIVOT COLUMN AND FOR EACH ROW WITH A NON-ZERO ADD 0000000
3387 C THE APPROPRIATE MULTIPLE OF THE PIVOT ROW TO IT. 0000000
3388 C WE LOOP ON THE NUMBER OF NON-ZEROS IN THE PIVOT COLUMN SINCE 0000000
3389 C MA30D/DD MAY CHANGE ITS ACTUAL POSITION. 0000000
3390 C 0000000
3391 NZPC = LENC(JPIV) 0000000
3392 IF (NZPC.EQ.0) GO TO 900 0000000
3393 DO 840 III=1,NZPC 0000000
3394 II = IPC(JPIV) + III - 1 0000000
3395 I = IRN(II) 0000000
3396 C SEARCH ROW I FOR NON-ZERO TO BE ELIMINATED, CALCULATE MULTIPLIER, 0000000
3397 C AND PLACE IT IN POSITION LENRL+1 IN ITS ROW. 0000000
3398 C IDROP IS THE NUMBER OF NON-ZERO ENTRIES DROPPED FROM ROW I 0000000
3399 C BECAUSE THESE FALL BENEATH TOLERANCE LEVEL. 0000000
3400 C 0000000
3401 IDROP = 0 0000000
3402 J1 = IPTR(I) + LENRL(I) 0000000
3403 IEND = IPTR(I) + LENR(I) - 1 0000000
3404 DO 570 JJ=J1,IEND 0000000
3405 IF (ICN(JJ).NE.JPIV) GO TO 570 0000000
3406 C IF PIVOT IS ZERO, REST OF COLUMN IS AND SO MULTIPLIER IS ZERO. 0000000
3407 AU = ZERO 0000000
3408 IF (A(IJPOS).NE.ZERO) AU = -A(JJ)/A(IJPOS) 0000000
3409 IF (LBIG) BIG = DMAX1(BIG,DABS(AU))
3410 A(JJ) = A(J1) 0000000
3411 A(J1) = AU 0000000
3412 ICN(JJ) = ICN(J1) 0000000
3413 ICN(J1) = JPIV 0000000
3414 LENRL(I) = LENRL(I) + 1 0000000
3415 GO TO 580 0000000
3416 570 CONTINUE 0000000
3417 C JUMP IF PIVOT ROW IS A SINGLETON. 0000000
3418 580 IF (LENPIV.EQ.0) GO TO 840 0000000
3419 C NOW PERFORM NECESSARY OPERATIONS ON REST OF NON-PIVOT ROW I. 0000000
3420 ROWI = J1 + 1 0000000
3421 IOP = 0 0000000
3422 C JUMP IF ALL THE PIVOT ROW CAUSES FILL-IN. 0000000
3423 IF (ROWI.GT.IEND) GO TO 650 0000000
3424 C PERFORM OPERATIONS ON CURRENT NON-ZEROS IN ROW I. 0000000
3425 C INNERMOST LOOP. 0000000
3426 DO 590 JJ=ROWI,IEND 0000000
3427 J = ICN(JJ) 0000000
3428 IF (IQ(J).GT.0) GO TO 590 0000000
3429 IOP = IOP + 1 0000000
3430 PIVROW = IJPOS - IQ(J) 0000000
3431 A(JJ) = A(JJ) + AU*A(PIVROW) 0000000
3432 IF (LBIG) BIG = DMAX1(DABS(A(JJ)),BIG)
3433 ICN(PIVROW) = -ICN(PIVROW) 0000000
3434 IF (DABS(A(JJ)).LT.TOL) IDROP = IDROP + 1
3435 590 CONTINUE 0000000
3436 C 0000000
3437 C JUMP IF NO NON-ZEROS IN NON-PIVOT ROW HAVE BEEN REMOVED 0000000
3438 C BECAUSE THESE ARE BENEATH THE DROP-TOLERANCE TOL. 0000000
3439 C 0000000
3440 IF (IDROP.EQ.0) GO TO 650 0000000
3441 C 0000000
3442 C RUN THROUGH NON-PIVOT ROW COMPRESSING ROW SO THAT ONLY 0000000
3443 C NON-ZEROS GREATER THAN TOL ARE STORED. ALL NON-ZEROS 0000000
3444 C LESS THAN TOL ARE ALSO REMOVED FROM THE COLUMN STRUCTURE. 0000000
3445 C 0000000
3446 JNEW = ROWI 0000000
3447 DO 630 JJ=ROWI,IEND 0000000
3448 IF (DABS(A(JJ)).LT.TOL) GO TO 600
3449 A(JNEW) = A(JJ) 0000000
3450 ICN(JNEW) = ICN(JJ) 0000000
3451 JNEW = JNEW + 1 0000000
3452 GO TO 630 0000000
3453 C 0000000
3454 C REMOVE NON-ZERO ENTRY FROM COLUMN STRUCTURE. 0000000
3455 C 0000000
3456 600 J = ICN(JJ) 0000000
3457 I1 = IPC(J) 0000000
3458 I2 = I1 + LENC(J) - 1 0000000
3459 DO 610 II=I1,I2 0000000
3460 IF (IRN(II).EQ.I) GO TO 620 0000000
3461 610 CONTINUE 0000000
3462 620 IRN(II) = IRN(I2) 0000000
3463 IRN(I2) = 0 0000000
3464 LENC(J) = LENC(J) - 1 0000000
3465 630 CONTINUE 0000000
3466 DO 640 JJ=JNEW,IEND 0000000
3467 ICN(JJ) = 0 0000000
3468 640 CONTINUE 0000000
3469 C THE VALUE OF IDROP MIGHT BE DIFFERENT FROM THAT CALCULATED EARLIER 0000000
3470 C BECAUSE, WE MAY NOW HAVE DROPPED SOME NON-ZEROS WHICH WERE NOT 0000000
3471 C MODIFIED BY THE PIVOT ROW. 0000000
3472 IDROP = IEND + 1 - JNEW 0000000
3473 IEND = JNEW - 1 0000000
3474 LENR(I) = LENR(I) - IDROP 0000000
3475 NZROW = NZROW - IDROP 0000000
3476 NZCOL = NZCOL - IDROP 0000000
3477 NDROP = NDROP + IDROP 0000000
3478 650 IFILL = LENPIV - IOP 0000000
3479 C JUMP IS IF THERE IS NO FILL-IN. 0000000
3480 IF (IFILL.EQ.0) GO TO 750 0000000
3481 C NOW FOR THE FILL-IN. 0000000
3482 MINICN = MAX0(MINICN,MOREI+IBEG-1+NZROW+IFILL+LENR(I)) 0000000
3483 C SEE IF THERE IS ROOM FOR FILL-IN. 0000000
3484 C GET MAXIMUM SPACE FOR ROW I IN SITU. 0000000
3485 DO 660 JDIFF=1,IFILL 0000000
3486 JNPOS = IEND + JDIFF 0000000
3487 IF (JNPOS.GT.LICN) GO TO 670 0000000
3488 IF (ICN(JNPOS).NE.0) GO TO 670 0000000
3489 660 CONTINUE 0000000
3490 C THERE IS ROOM FOR ALL THE FILL-IN AFTER THE END OF THE ROW SO IT 0000000
3491 C CAN BE LEFT IN SITU. 0000000
3492 C NEXT AVAILABLE SPACE FOR FILL-IN. 0000000
3493 IEND = IEND + 1 0000000
3494 GO TO 750 0000000
3495 C JMORE SPACES FOR FILL-IN ARE REQUIRED IN FRONT OF ROW. 0000000
3496 670 JMORE = IFILL - JDIFF + 1 0000000
3497 I1 = IPTR(I) 0000000
3498 C WE NOW LOOK IN FRONT OF THE ROW TO SEE IF THERE IS SPACE FOR 0000000
3499 C THE REST OF THE FILL-IN. 0000000
3500 DO 680 JDIFF=1,JMORE 0000000
3501 JNPOS = I1 - JDIFF 0000000
3502 IF (JNPOS.LT.IACTIV) GO TO 690 0000000
3503 IF (ICN(JNPOS).NE.0) GO TO 700 0000000
3504 680 CONTINUE 0000000
3505 690 JNPOS = I1 - JMORE 0000000
3506 GO TO 710 0000000
3507 C WHOLE ROW MUST BE MOVED TO THE BEGINNING OF AVAILABLE STORAGE. 0000000
3508 700 JNPOS = IACTIV - LENR(I) - IFILL 0000000
3509 C JUMP IF THERE IS SPACE IMMEDIATELY AVAILABLE FOR THE SHIFTED ROW. 0000000
3510 710 IF (JNPOS.GE.IBEG) GO TO 730 0000000
3511 CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
3512 I1 = IPTR(I) 0000000
3513 IEND = I1 + LENR(I) - 1 0000000
3514 JNPOS = IACTIV - LENR(I) - IFILL 0000000
3515 IF (JNPOS.GE.IBEG) GO TO 730 0000000
3516 C NO SPACE AVAILABLE SO TRY TO CREATE SOME BY THROWING AWAY PREVIOUS 0000000
3517 C LU DECOMPOSITION. 0000000
3518 MOREI = MOREI + IBEG - IDISP(1) - LENPIV - 1 0000000
3519 IF (LP.NE.0) WRITE (LP,99997) 0000000
3520 IFLAG = -5 0000000
3521 IF (ABORT3) GO TO 1090 0000000
3522 C KEEP RECORD OF CURRENT PIVOT ROW. 0000000
3523 IBEG = IDISP(1) 0000000
3524 ICN(IBEG) = JPIV 0000000
3525 A(IBEG) = A(IJPOS) 0000000
3526 IJPOS = IBEG 0000000
3527 DO 720 JJ=IJP1,PIVEND 0000000
3528 IBEG = IBEG + 1 0000000
3529 A(IBEG) = A(JJ) 0000000
3530 ICN(IBEG) = ICN(JJ) 0000000
3531 720 CONTINUE 0000000
3532 IJP1 = IJPOS + 1 0000000
3533 PIVEND = IBEG 0000000
3534 IBEG = IBEG + 1 0000000
3535 IF (JNPOS.GE.IBEG) GO TO 730 0000000
3536 C THIS STILL DOES NOT GIVE ENOUGH ROOM. 0000000
3537 IFLAG = -4 0000000
3538 GO TO 1090 0000000
3539 730 IACTIV = MIN0(IACTIV,JNPOS) 0000000
3540 C MOVE NON-PIVOT ROW I. 0000000
3541 IPTR(I) = JNPOS 0000000
3542 DO 740 JJ=I1,IEND 0000000
3543 A(JNPOS) = A(JJ) 0000000
3544 ICN(JNPOS) = ICN(JJ) 0000000
3545 JNPOS = JNPOS + 1 0000000
3546 ICN(JJ) = 0 0000000
3547 740 CONTINUE 0000000
3548 C FIRST NEW AVAILABLE SPACE. 0000000
3549 IEND = JNPOS 0000000
3550 750 NZROW = NZROW + IFILL 0000000
3551 C INNERMOST FILL-IN LOOP WHICH ALSO RESETS ICN. 0000000
3552 DO 830 JJ=IJP1,PIVEND 0000000
3553 J = ICN(JJ) 0000000
3554 IF (J.LT.0) GO TO 820 0000000
3555 ANEW = AU*A(JJ) 0000000
3556 AANEW = DABS(ANEW)
3557 IF (AANEW.GE.TOL) GO TO 760 0000000
3558 NDROP = NDROP + 1 0000000
3559 NZROW = NZROW - 1 0000000
3560 MINICN = MINICN - 1 0000000
3561 IFILL = IFILL - 1 0000000
3562 GO TO 830 0000000
3563 760 IF (LBIG) BIG = DMAX1(AANEW,BIG)
3564 A(IEND) = ANEW 0000000
3565 ICN(IEND) = J 0000000
3566 IEND = IEND + 1 0000000
3567 C 0000000
3568 C PUT NEW ENTRY IN COLUMN FILE. 0000000
3569 MINIRN = MAX0(MINIRN,NZCOL+LENC(J)+1) 0000000
3570 JEND = IPC(J) + LENC(J) 0000000
3571 JROOM = NZPC - III + 1 + LENC(J) 0000000
3572 IF (JEND.GT.LIRN) GO TO 770 0000000
3573 IF (IRN(JEND).EQ.0) GO TO 810 0000000
3574 770 IF (JROOM.LT.DISPC) GO TO 780 0000000
3575 C COMPRESS COLUMN FILE TO OBTAIN SPACE FOR NEW COPY OF COLUMN. 0000000
3576 CALL MA30DD(A, IRN, IPC(ISTART), N, DISPC, LIRN, .FALSE.)
3577 IF (JROOM.LT.DISPC) GO TO 780 0000000
3578 JROOM = DISPC - 1 0000000
3579 IF (JROOM.GE.LENC(J)+1) GO TO 780 0000000
3580 C COLUMN FILE IS NOT LARGE ENOUGH. 0000000
3581 GO TO 1100 0000000
3582 C COPY COLUMN TO BEGINNING OF FILE. 0000000
3583 780 JBEG = IPC(J) 0000000
3584 JEND = IPC(J) + LENC(J) - 1 0000000
3585 JZERO = DISPC - 1 0000000
3586 DISPC = DISPC - JROOM 0000000
3587 IDISPC = DISPC 0000000
3588 DO 790 II=JBEG,JEND 0000000
3589 IRN(IDISPC) = IRN(II) 0000000
3590 IRN(II) = 0 0000000
3591 IDISPC = IDISPC + 1 0000000
3592 790 CONTINUE 0000000
3593 IPC(J) = DISPC 0000000
3594 JEND = IDISPC 0000000
3595 DO 800 II=JEND,JZERO 0000000
3596 IRN(II) = 0 0000000
3597 800 CONTINUE 0000000
3598 810 IRN(JEND) = I 0000000
3599 NZCOL = NZCOL + 1 0000000
3600 LENC(J) = LENC(J) + 1 0000000
3601 C END OF ADJUSTMENT TO COLUMN FILE. 0000000
3602 GO TO 830 0000000
3603 C 0000000
3604 820 ICN(JJ) = -J 0000000
3605 830 CONTINUE 0000000
3606 LENR(I) = LENR(I) + IFILL 0000000
3607 C END OF SCAN OF PIVOT COLUMN. 0000000
3608 840 CONTINUE 0000000
3609 C 0000000
3610 C 0000000
3611 C REMOVE PIVOT COLUMN FROM COLUMN ORIENTED STORAGE AND UPDATE ROW 0000000
3612 C ORDERING ARRAYS. 0000000
3613 I1 = IPC(JPIV) 0000000
3614 I2 = IPC(JPIV) + LENC(JPIV) - 1 0000000
3615 NZCOL = NZCOL - LENC(JPIV) 0000000
3616 DO 890 II=I1,I2 0000000
3617 I = IRN(II) 0000000
3618 IRN(II) = 0 0000000
3619 NZ = LENR(I) - LENRL(I) 0000000
3620 IF (NZ.NE.0) GO TO 850 0000000
3621 LASTR(I) = 0 0000000
3622 GO TO 890 0000000
3623 850 IFIR = IFIRST(NZ) 0000000
3624 IFIRST(NZ) = I 0000000
3625 IF (IFIR) 860, 880, 870 0000000
3626 860 LASTR(I) = IFIR 0000000
3627 NEXTR(I) = 0 0000000
3628 GO TO 890 0000000
3629 870 LASTR(I) = LASTR(IFIR) 0000000
3630 NEXTR(I) = IFIR 0000000
3631 LASTR(IFIR) = I 0000000
3632 GO TO 890 0000000
3633 880 LASTR(I) = 0 0000000
3634 NEXTR(I) = 0 0000000
3635 NZMIN = MIN0(NZMIN,NZ) 0000000
3636 890 CONTINUE 0000000
3637 C RESTORE IQ AND NULLIFY U PART OF OLD PIVOT ROW. 0000000
3638 C RECORD THE COLUMN PERMUTATION IN LASTC(JPIV) AND THE ROW 0000000
3639 C PERMUTATION IN LASTR(IPIV). 0000000
3640 900 IPC(JPIV) = -ISING 0000000
3641 LASTR(IPIV) = PIVOT 0000000
3642 IF (LENPIV.EQ.0) GO TO 980 0000000
3643 NZROW = NZROW - LENPIV 0000000
3644 JVAL = IJP1 0000000
3645 JZER = IPTR(IPIV) 0000000
3646 IPTR(IPIV) = 0 0000000
3647 DO 910 JCOUNT=1,LENPIV 0000000
3648 J = ICN(JVAL) 0000000
3649 IQ(J) = ICN(JZER) 0000000
3650 ICN(JZER) = 0 0000000
3651 JVAL = JVAL + 1 0000000
3652 JZER = JZER + 1 0000000
3653 910 CONTINUE 0000000
3654 C ADJUST COLUMN ORDERING ARRAYS. 0000000
3655 DO 970 JJ=IJP1,PIVEND 0000000
3656 J = ICN(JJ) 0000000
3657 NZ = LENC(J) 0000000
3658 IF (NZ.NE.0) GO TO 920 0000000
3659 IPC(J) = 0 0000000
3660 GO TO 970 0000000
3661 920 IF (NSRCH.LE.NN) GO TO 960 0000000
3662 IFIR = IFIRST(NZ) 0000000
3663 LASTC(J) = 0 0000000
3664 IF (IFIR) 930, 940, 950 0000000
3665 930 IFIRST(NZ) = -J 0000000
3666 IFIR = -IFIR 0000000
3667 LASTC(IFIR) = J 0000000
3668 NEXTC(J) = IFIR 0000000
3669 GO TO 970 0000000
3670 940 IFIRST(NZ) = -J 0000000
3671 NEXTC(J) = 0 0000000
3672 GO TO 960 0000000
3673 950 LC = -LASTR(IFIR) 0000000
3674 LASTR(IFIR) = -J 0000000
3675 NEXTC(J) = LC 0000000
3676 IF (LC.NE.0) LASTC(LC) = J 0000000
3677 960 NZMIN = MIN0(NZMIN,NZ) 0000000
3678 970 CONTINUE 0000000
3679 980 CONTINUE 0000000
3680 C ******************************************** 0000000
3681 C **** END OF MAIN ELIMINATION LOOP **** 0000000
3682 C ******************************************** 0000000
3683 C 0000000
3684 C RESET IACTIV TO POINT TO THE BEGINNING OF THE NEXT BLOCK. 0000000
3685 990 IF (ILAST.NE.NN) IACTIV = IPTR(ILAST+1) 0000000
3686 1000 CONTINUE 0000000
3687 C 0000000
3688 C ******************************************** 0000000
3689 C **** END OF DEOMPOSITION OF BLOCK **** 0000000
3690 C ******************************************** 0000000
3691 C 0000000
3692 C RECORD SINGULARITY (IF ANY) IN IQ ARRAY. 0000000
3693 IF (IRANK.EQ.NN) GO TO 1020 0000000
3694 DO 1010 I=1,NN 0000000
3695 IF (IPC(I).LT.0) GO TO 1010 0000000
3696 ISING = IPC(I) 0000000
3697 IQ(ISING) = -IQ(ISING) 0000000
3698 IPC(I) = -ISING 0000000
3699 1010 CONTINUE 0000000
3700 C 0000000
3701 C RUN THROUGH LU DECOMPOSITION CHANGING COLUMN INDICES TO THAT OF NEW 0000000
3702 C ORDER AND PERMUTING LENR AND LENRL ARRAYS ACCORDING TO PIVOT 0000000
3703 C PERMUTATIONS. 0000000
3704 1020 ISTART = IDISP(1) 0000000
3705 IEND = IBEG - 1 0000000
3706 IF (IEND.LT.ISTART) GO TO 1040 0000000
3707 DO 1030 JJ=ISTART,IEND 0000000
3708 JOLD = ICN(JJ) 0000000
3709 ICN(JJ) = -IPC(JOLD) 0000000
3710 1030 CONTINUE 0000000
3711 1040 DO 1050 II=1,NN 0000000
3712 I = LASTR(II) 0000000
3713 NEXTR(I) = LENR(II) 0000000
3714 IPTR(I) = LENRL(II) 0000000
3715 1050 CONTINUE 0000000
3716 DO 1060 I=1,NN 0000000
3717 LENRL(I) = IPTR(I) 0000000
3718 LENR(I) = NEXTR(I) 0000000
3719 1060 CONTINUE 0000000
3720 C 0000000
3721 C UPDATE PERMUTATION ARRAYS IP AND IQ. 0000000
3722 DO 1070 II=1,NN 0000000
3723 I = LASTR(II) 0000000
3724 J = -IPC(II) 0000000
3725 NEXTR(I) = IABS(IP(II)+0) 0000000
3726 IPTR(J) = IABS(IQ(II)+0) 0000000
3727 1070 CONTINUE 0000000
3728 DO 1080 I=1,NN 0000000
3729 IF (IP(I).LT.0) NEXTR(I) = -NEXTR(I) 0000000
3730 IP(I) = NEXTR(I) 0000000
3731 IF (IQ(I).LT.0) IPTR(I) = -IPTR(I) 0000000
3732 IQ(I) = IPTR(I) 0000000
3733 1080 CONTINUE 0000000
3734 IP(NN) = IABS(IP(NN)+0) 0000000
3735 IDISP(2) = IEND 0000000
3736 GO TO 1120 0000000
3737 C 0000000
3738 C *** ERROR RETURNS *** 0000000
3739 1090 IDISP(2) = IACTIV 0000000
3740 IF (LP.EQ.0) GO TO 1120 0000000
3741 WRITE (LP,99996) 0000000
3742 GO TO 1110 0000000
3743 1100 IF (IFLAG.EQ.-5) IFLAG = -6 0000000
3744 IF (IFLAG.NE.-6) IFLAG = -3 0000000
3745 IDISP(2) = IACTIV 0000000
3746 IF (LP.EQ.0) GO TO 1120 0000000
3747 IF (IFLAG.EQ.-3) WRITE (LP,99995) 0000000
3748 IF (IFLAG.EQ.-6) WRITE (LP,99994) 0000000
3749 1110 PIVOT = PIVOT - ISTART + 1 0000000
3750 WRITE (LP,99993) PIVOT, NBLOCK, ISTART, ILAST 0000000
3751 IF (PIVOT.EQ.0) WRITE (LP,99992) MINIRN 0000000
3752 C 0000000
3753 C 0000000
3754 1120 RETURN 0000000
3755 99999 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE MATRIX IS STRUCTUR,0000000
3756 * 13HALLY SINGULAR) 0000000
3757 99998 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE MATRIX IS NUMERICA,0000000
3758 * 12HLLY SINGULAR) 0000000
3759 99997 FORMAT (48H LU DECOMPOSITION DESTROYED TO CREATE MORE SPACE) 0000000
3760 99996 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE LICN NOT BIG ENOUG,0000000
3761 * 1HH) 0000000
3762 99995 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE LIRN NOT BIG ENOUG,0000000
3763 * 1HH) 0000000
3764 99994 FORMAT (51H ERROR RETURN FROM MA30A/AD LIRN AND LICN TOO SMALL) 0000000
3765 99993 FORMAT (10H AT STAGE , I5, 10H IN BLOCK , I5, 16H WITH FIRST ROW ,0000000
3766 * I5, 14H AND LAST ROW , I5) 0000000
3767 99992 FORMAT (34H TO CONTINUE SET LIRN TO AT LEAST , I8) 0000000
3768 END 0000000
3769 SUBROUTINE MA30DD(A, ICN, IPTR, N, IACTIV, ITOP, REALS)
3770 c_270390
3771 EXTERNAL MA30$DATA
3772 c_270390
3773 C THIS SUBROUTINE PERFORMS GARBAGE COLLECTION OPERATIONS ON THE 0000000
3774 C ARRAYS A, ICN AND IRN. 0000000
3775 C IACTIV IS THE FIRST POSITION IN ARRAYS A/ICN FROM WHICH THE COMPRESS 0000000
3776 C STARTS. ON EXIT, IACTIV EQUALS THE POSITION OF THE FIRST ENTRY 0000000
3777 C IN THE COMPRESSED PART OF A/ICN 0000000
3778 C 0000000
3779 DOUBLE PRECISION A(ITOP)
3780 LOGICAL REALS 0000000
3781 INTEGER IPTR(N) 0000000
3782 Change
3783 C INTEGER*2 ICN(ITOP)
3784 INTEGER ICN(ITOP)
3785 Change
3786 C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON.
3787 COMMON /MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN
3788 C 0000000
3789 IF (REALS) ICNCP = ICNCP + 1 0000000
3790 IF (.NOT.REALS) IRNCP = IRNCP + 1 0000000
3791 C SET THE FIRST NON-ZERO ENTRY IN EACH ROW TO THE NEGATIVE OF THE 0000000
3792 C ROW/COL NUMBER AND HOLD THIS ROW/COL INDEX IN THE ROW/COL 0000000
3793 C POINTER. THIS IS SO THAT THE BEGINNING OF EACH ROW/COL CAN 0000000
3794 C BE RECOGNIZED IN THE SUBSEQUENT SCAN. 0000000
3795 DO 10 J=1,N 0000000
3796 K = IPTR(J) 0000000
3797 IF (K.LT.IACTIV) GO TO 10 0000000
3798 IPTR(J) = ICN(K) 0000000
3799 ICN(K) = -J 0000000
3800 10 CONTINUE 0000000
3801 KN = ITOP + 1 0000000
3802 KL = ITOP - IACTIV + 1 0000000
3803 C GO THROUGH ARRAYS IN REVERSE ORDER COMPRESSING TO THE BACK SO 0000000
3804 C THAT THERE ARE NO ZEROS HELD IN POSITIONS IACTIV TO ITOP IN ICN. 0000000
3805 C RESET FIRST ENTRY OF EACH ROW/COL AND POINTER ARRAY IPTR. 0000000
3806 DO 30 K=1,KL 0000000
3807 JPOS = ITOP - K + 1 0000000
3808 IF (ICN(JPOS).EQ.0) GO TO 30 0000000
3809 KN = KN - 1 0000000
3810 IF (REALS) A(KN) = A(JPOS) 0000000
3811 IF (ICN(JPOS).GE.0) GO TO 20 0000000
3812 C FIRST NON-ZERO OF ROW/COL HAS BEEN LOCATED 0000000
3813 J = -ICN(JPOS) 0000000
3814 ICN(JPOS) = IPTR(J) 0000000
3815 IPTR(J) = KN 0000000
3816 20 ICN(KN) = ICN(JPOS) 0000000
3817 30 CONTINUE 0000000
3818 IACTIV = KN 0000000
3819 RETURN 0000000
3820 END 0000000
3821 SUBROUTINE MA30BD(N, ICN, A, LICN, LENR, LENRL, IDISP, IP, IQ, W,
3822 * IW, IFLAG) 0000000
3823 c_270390
3824 EXTERNAL MA30$DATA
3825 c_270390
3826 C MA30B/BD PERFORMS THE LU DECOMPOSITION OF THE DIAGONAL BLOCKS OF A 0000000
3827 C NEW MATRIX PAQ OF THE SAME SPARSITY PATTERN, USING INFORMATION 0000000
3828 C FROM A PREVIOUS CALL TO MA30A/AD. THE ENTRIES OF THE INPUT 0000000
3829 C MATRIX MUST ALREADY BE IN THEIR FINAL POSITIONS IN THE LU 0000000
3830 C DECOMPOSITION STRUCTURE. THIS ROUTINE EXECUTES ABOUT FIVE TIMES 0000000
3831 C FASTER THAN MA30A/AD. 0000000
3832 C 0000000
3833 C WE NOW DESCRIBE THE ARGUMENT LIST FOR MA30B/BD. CONSULT MA30A/AD FOR 0000000
3834 C FURTHER INFORMATION ON THESE PARAMETERS. 0000000
3835 C N IS AN INTEGER VARIABLE SET TO THE ORDER OF THE MATRIX. 0000000
3836 C ICN IS AN INTEGER*2 ARRAY OF LENGTH LICN. IT SHOULD BE UNCHANGED
3837 C SINCE THE LAST CALL TO MA30A/AD. IT IS NOT ALTERED BY MA30B/BD. 0000000
3838 C A IS A REAL/DOUBLE PRECISION ARRAY OF LENGTH LICN THE USER MUST SET 0000000
3839 C ENTRIES IDISP(1) TO IDISP(2) TO CONTAIN THE ENTRIES IN THE 0000000
3840 C DIAGONAL BLOCKS OF THE MATRIX PAQ WHOSE COLUMN NUMBERS ARE HELD 0000000
3841 C IN ICN, USING CORRESPONDING POSITIONS. NOTE THAT SOME ZEROS MAY 0000000
3842 C NEED TO BE HELD EXPLICITLY. ON OUTPUT ENTRIES IDISP(1) TO 0000000
3843 C IDISP(2) OF ARRAY A CONTAIN THE LU DECOMPOSITION OF THE DIAGONAL 0000000
3844 C BLOCKS OF PAQ. ENTRIES A(1) TO A(IDISP(1)-1) ARE NEITHER 0000000
3845 C REQUIRED NOR ALTERED BY MA30B/BD. 0000000
3846 C LICN IS AN INTEGER VARIABLE WHICH MUST BE SET BY THE USER TO THE 0000000
3847 C LENGTH OF ARRAYS A AND ICN. IT IS NOT ALTERED BY MA30B/BD. 0000000
3848 C LENR,LENRL ARE INTEGER*2 ARRAYS OF LENGTH N. THEY SHOULD BE
3849 C UNCHANGED SINCE THE LAST CALL TO MA30A/AD. THEY ARE NOT ALTERED 0000000
3850 C BY MA30B/BD. 0000000
3851 C IDISP IS AN INTEGER ARRAY OF LENGTH 2. IT SHOULD BE UNCHANGED SINCE 0000000
3852 C THE LAST CALL TO MA30A/AD. IT IS NOT ALTERED BY MA30B/BD. 0000000
3853 C IP,IQ ARE INTEGER*2 ARRAYS OF LENGTH N. THEY SHOULD BE UNCHANGED
3854 C SINCE THE LAST CALL TO MA30A/AD. THEY ARE NOT ALTERED BY 0000000
3855 C MA30B/BD. 0000000
3856 C W IS A REAL/DOUBLE PRECISION ARRAY OF LENGTH N WHICH IS USED AS 0000000
3857 C WORKSPACE BY MA30B/BD. 0000000
3858 C IW IS AN INTEGER ARRAY OF LENGTH N WHICH IS USED AS WORKSPACE BY 0000000
3859 C MA30B/BD. 0000000
3860 C IFLAG IS AN INTEGER VARIABLE. ON OUTPUT FROM MA30B/BD, IFLAG HAS 0000000
3861 C THE VALUE ZERO IF THE FACTORIZATION WAS SUCCESSFUL, HAS THE 0000000
3862 C VALUE I IF PIVOT I WAS VERY SMALL AND HAS THE VALUE -I IF AN 0000000
3863 C UNEXPECTED SINGULARITY WAS DETECTED AT STAGE I OF THE 0000000
3864 C DECOMPOSITION. 0000000
3865 C 0000000
3866 DOUBLE PRECISION A(LICN), W(N), AU, EPS, ROWMAX, ZERO, ONE, RMIN,
3867 * TOL, BIG
3868 LOGICAL ABORT1, ABORT2, ABORT3, STAB, LBIG 0000000
3869 INTEGER IW(N), IDISP(2), PIVPOS 0000000
3870 Change
3871 C INTEGER*2 ICN(LICN), LENR(N), LENRL(N), IP(N), IQ(N)
3872 INTEGER ICN(LICN), LENR(N), LENRL(N), IP(N), IQ(N)
3873 Change
3874 C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON.
3875 COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
3876 COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
3877 COMMON /MA30GD/ EPS, RMIN
3878 DATA ZERO /0.0D0/, ONE /1.0D0/
3879 STAB = EPS.LE.ONE 0000000
3880 RMIN = EPS 0000000
3881 ISING = 0 0000000
3882 IFLAG = 0 0000000
3883 DO 10 I=1,N 0000000
3884 W(I) = ZERO 0000000
3885 10 CONTINUE 0000000
3886 C SET UP POINTERS TO THE BEGINNING OF THE ROWS. 0000000
3887 IW(1) = IDISP(1) 0000000
3888 IF (N.EQ.1) GO TO 25 0000000
3889 DO 20 I=2,N 0000000
3890 IW(I) = IW(I-1) + LENR(I-1) 0000000
3891 20 CONTINUE 0000000
3892 C 0000000
3893 C **** START OF MAIN LOOP **** 0000000
3894 C AT STEP I, ROW I OF A IS TRANSFORMED TO ROW I OF L/U BY ADDING 0000000
3895 C APPROPRIATE MULTIPLES OF ROWS 1 TO I-1. 0000000
3896 C .... USING ROW-GAUSS ELIMINATION. 0000000
3897 25 DO 160 I=1,N 0000000
3898 C ISTART IS BEGINNING OF ROW I OF A AND ROW I OF L. 0000000
3899 ISTART = IW(I) 0000000
3900 C IFIN IS END OF ROW I OF A AND ROW I OF U. 0000000
3901 IFIN = ISTART + LENR(I) - 1 0000000
3902 C ILEND IS END OF ROW I OF L. 0000000
3903 ILEND = ISTART + LENRL(I) - 1 0000000
3904 IF (ISTART.GT.ILEND) GO TO 90 0000000
3905 C LOAD ROW I OF A INTO VECTOR W. 0000000
3906 DO 30 JJ=ISTART,IFIN 0000000
3907 J = ICN(JJ) 0000000
3908 W(J) = A(JJ) 0000000
3909 30 CONTINUE 0000000
3910 C 0000000
3911 C ADD MULTIPLES OF APPROPRIATE ROWS OF I TO I-1 TO ROW I. 0000000
3912 DO 70 JJ=ISTART,ILEND 0000000
3913 J = ICN(JJ) 0000000
3914 C IPIVJ IS POSITION OF PIVOT IN ROW J. 0000000
3915 IPIVJ = IW(J) + LENRL(J) 0000000
3916 C FORM MULTIPLIER AU. 0000000
3917 AU = -W(J)/A(IPIVJ) 0000000
3918 IF (LBIG) BIG = DMAX1(DABS(AU),BIG)
3919 W(J) = AU 0000000
3920 C AU * ROW J (U PART) IS ADDED TO ROW I. 0000000
3921 IPIVJ = IPIVJ + 1 0000000
3922 JFIN = IW(J) + LENR(J) - 1 0000000
3923 IF (IPIVJ.GT.JFIN) GO TO 70 0000000
3924 C INNERMOST LOOP. 0000000
3925 IF (LBIG) GO TO 50 0000000
3926 DO 40 JAYJAY=IPIVJ,JFIN 0000000
3927 JAY = ICN(JAYJAY) 0000000
3928 W(JAY) = W(JAY) + AU*A(JAYJAY) 0000000
3929 40 CONTINUE 0000000
3930 GO TO 70 0000000
3931 50 DO 60 JAYJAY=IPIVJ,JFIN 0000000
3932 JAY = ICN(JAYJAY) 0000000
3933 W(JAY) = W(JAY) + AU*A(JAYJAY) 0000000
3934 BIG = DMAX1(DABS(W(JAY)),BIG)
3935 60 CONTINUE 0000000
3936 70 CONTINUE 0000000
3937 C 0000000
3938 C RELOAD W BACK INTO A (NOW L/U) 0000000
3939 DO 80 JJ=ISTART,IFIN 0000000
3940 J = ICN(JJ) 0000000
3941 A(JJ) = W(J) 0000000
3942 W(J) = ZERO 0000000
3943 80 CONTINUE 0000000
3944 C WE NOW PERFORM THE STABILITY CHECKS. 0000000
3945 90 PIVPOS = ILEND + 1 0000000
3946 IF (IQ(I).GT.0) GO TO 140 0000000
3947 C MATRIX HAD SINGULARITY AT THIS POINT IN MA30A/AD. 0000000
3948 C IS IT THE FIRST SUCH PIVOT IN CURRENT BLOCK ? 0000000
3949 IF (ISING.EQ.0) ISING = I 0000000
3950 C DOES CURRENT MATRIX HAVE A SINGULARITY IN THE SAME PLACE ? 0000000
3951 IF (PIVPOS.GT.IFIN) GO TO 100 0000000
3952 IF (A(PIVPOS).NE.ZERO) GO TO 170 0000000
3953 C IT DOES .. SO SET ISING IF IT IS NOT THE END OF THE CURRENT BLOCK 0000000
3954 C CHECK TO SEE THAT APPROPRIATE PART OF L/U IS ZERO OR NULL. 0000000
3955 100 IF (ISTART.GT.IFIN) GO TO 120 0000000
3956 DO 110 JJ=ISTART,IFIN 0000000
3957 IF (ICN(JJ).LT.ISING) GO TO 110 0000000
3958 IF (A(JJ).NE.ZERO) GO TO 170 0000000
3959 110 CONTINUE 0000000
3960 120 IF (PIVPOS.LE.IFIN) A(PIVPOS) = ONE 0000000
3961 IF (IP(I).GT.0 .AND. I.NE.N) GO TO 160 0000000
3962 C END OF CURRENT BLOCK ... RESET ZERO PIVOTS AND ISING. 0000000
3963 DO 130 J=ISING,I 0000000
3964 IF ((LENR(J)-LENRL(J)).EQ.0) GO TO 130 0000000
3965 JJ = IW(J) + LENRL(J) 0000000
3966 A(JJ) = ZERO 0000000
3967 130 CONTINUE 0000000
3968 ISING = 0 0000000
3969 GO TO 160 0000000
3970 C MATRIX HAD NON-ZERO PIVOT IN MA30A/AD AT THIS STAGE. 0000000
3971 140 IF (PIVPOS.GT.IFIN) GO TO 170 0000000
3972 IF (A(PIVPOS).EQ.ZERO) GO TO 170 0000000
3973 IF (.NOT.STAB) GO TO 160 0000000
3974 ROWMAX = ZERO 0000000
3975 DO 150 JJ=PIVPOS,IFIN 0000000
3976 ROWMAX = DMAX1(ROWMAX,DABS(A(JJ)))
3977 150 CONTINUE 0000000
3978 IF (DABS(A(PIVPOS))/ROWMAX.GE.RMIN) GO TO 160
3979 IFLAG = I 0000000
3980 RMIN = DABS(A(PIVPOS))/ROWMAX
3981 C **** END OF MAIN LOOP **** 0000000
3982 160 CONTINUE 0000000
3983 C 0000000
3984 GO TO 180 0000000
3985 C *** ERROR RETURN *** 0000000
3986 170 IF (LP.NE.0) WRITE (LP,99999) I 0000000
3987 IFLAG = -I 0000000
3988 C 0000000
3989 180 RETURN 0000000
3990 99999 FORMAT (54H ERROR RETURN FROM MA30B/BD SINGULARITY DETECTED IN RO,0000000
3991 * 1HW, I8) 0000000
3992 END 0000000
3993 SUBROUTINE MA30CD(N, ICN, A, LICN, LENR, LENRL, LENOFF, IDISP, IP,
3994 * IQ, X, W, MTYPE) 0000000
3995 c_270390
3996 EXTERNAL MA30$DATA
3997 c_270390
3998 C MA30C/CD USES THE FACTORS PRODUCED BY MA30A/AD OR MA30B/BD TO SOLVE 0000000
3999 C AX=B OR A TRANSPOSE X=B WHEN THE MATRIX P1*A*Q1 (PAQ) IS BLOCK 0000000
4000 C LOWER TRIANGULAR (INCLUDING THE CASE OF ONLY ONE DIAGONAL 0000000
4001 C BLOCK). 0000000
4002 C 0000000
4003 C WE NOW DESCRIBE THE ARGUMENT LIST FOR MA30C/CD. 0000000
4004 C N IS AN INTEGER VARIABLE SET TO THE ORDER OF THE MATRIX. IT IS NOT 0000000
4005 C ALTERED BY THE SUBROUTINE. 0000000
4006 C ICN IS AN INTEGER*2 ARRAY OF LENGTH LICN. ENTRIES IDISP(1) TO
4007 C IDISP(2) SHOULD BE UNCHANGED SINCE THE LAST CALL TO MA30A/AD. IF 0000000
4008 C THE MATRIX HAS MORE THAN ONE DIAGONAL BLOCK, THEN COLUMN INDICES 0000000
4009 C CORRESPONDING TO NON-ZEROS IN SUB-DIAGONAL BLOCKS OF PAQ MUST 0000000
4010 C APPEAR IN POSITIONS 1 TO IDISP(1)-1. FOR THE SAME ROW THOSE 0000000
4011 C ENTRIES MUST BE CONTIGUOUS, WITH THOSE IN ROW I PRECEDING THOSE 0000000
4012 C IN ROW I+1 (I=1,...,N-1) AND NO WASTED SPACE BETWEEN ROWS. 0000000
4013 C ENTRIES MAY BE IN ANY ORDER WITHIN EACH ROW. IT IS NOT ALTERED 0000000
4014 C BY MA30C/CD. 0000000
4015 C A IS A REAL/DOUBLE PRECISION ARRAY OF LENGTH LICN. ENTRIES 0000000
4016 C IDISP(1) TO IDISP(2) SHOULD BE UNCHANGED SINCE THE LAST CALL TO 0000000
4017 C MA30A/AD OR MA30B/BD. IF THE MATRIX HAS MORE THAN ONE DIAGONAL 0000000
4018 C BLOCK, THEN THE VALUES OF THE NON-ZEROS IN SUB-DIAGONAL BLOCKS 0000000
4019 C MUST BE IN POSITIONS 1 TO IDISP(1)-1 IN THE ORDER GIVEN BY ICN. 0000000
4020 C IT IS NOT ALTERED BY MA30C/CD. 0000000
4021 C LICN IS AN INTEGER VARIABLE SET TO THE SIZE OF ARRAYS ICN AND A. 0000000
4022 C IT IS NOT ALTERED BY MA30C/CD. 0000000
4023 C LENR,LENRL ARE INTEGER*2 ARRAYS OF LENGTH N WHICH SHOULD BE
4024 C UNCHANGED SINCE THE LAST CALL TO MA30A/AD. THEY ARE NOT ALTERED 0000000
4025 C BY MA30C/CD. 0000000
4026 C LENOFF IS AN INTEGER*2 ARRAY OF LENGTH N. IF THE MATRIX PAQ (OR
4027 C P1*A*Q1) HAS MORE THAN ONE DIAGONAL BLOCK, THEN LENOFF(I), 0000000
4028 C I=1,...,N SHOULD BE SET TO THE NUMBER OF NON-ZEROS IN ROW I OF 0000000
4029 C THE MATRIX PAQ WHICH ARE IN SUB-DIAGONAL BLOCKS. IF THERE IS 0000000
4030 C ONLY ONE DIAGONAL BLOCK THEN LENOFF(1) MAY BE SET TO -1, IN 0000000
4031 C WHICH CASE THE OTHER ENTRIES OF LENOFF ARE NEVER ACCESSED. IT IS 0000000
4032 C NOT ALTERED BY MA30C/CD. 0000000
4033 C IDISP IS AN INTEGER ARRAY OF LENGTH 2 WHICH SHOULD BE UNCHANGED 0000000
4034 C SINCE THE LAST CALL TO MA30A/AD. IT IS NOT ALTERED BY MA30C/CD. 0000000
4035 C IP,IQ ARE INTEGER*2 ARRAYS OF LENGTH N WHICH SHOULD BE UNCHANGED
4036 C SINCE THE LAST CALL TO MA30A/AD. THEY ARE NOT ALTERED BY 0000000
4037 C MA30C/CD. 0000000
4038 C X IS A REAL/DOUBLE PRECISION ARRAY OF LENGTH N. IT MUST BE SET BY 0000000
4039 C THE USER TO THE VALUES OF THE RIGHT HAND SIDE VECTOR B FOR THE 0000000
4040 C EQUATIONS BEING SOLVED. ON EXIT FROM MA30C/CD IT WILL BE EQUAL 0000000
4041 C TO THE SOLUTION X REQUIRED. 0000000
4042 C W IS A REAL/DOUBLE PRECISION ARRAY OF LENGTH N WHICH IS USED AS 0000000
4043 C WORKSPACE BY MA30C/CD. 0000000
4044 C MTYPE IS AN INTEGER VARIABLE WHICH MUST BE SET BY THE USER. IF 0000000
4045 C MTYPE=1, THEN THE SOLUTION TO THE SYSTEM AX=B IS RETURNED; ANY 0000000
4046 C OTHER VALUE FOR MTYPE WILL RETURN THE SOLUTION TO THE SYSTEM A 0000000
4047 C TRANSPOSE X=B. IT IS NOT ALTERED BY MA30C/CD. 0000000
4048 C 0000000
4049 DOUBLE PRECISION A(LICN), X(N), W(N), WII, WI, RESID, ZERO
4050 LOGICAL NEG, NOBLOC 0000000
4051 INTEGER IDISP(2) 0000000
4052 Change
4053 C INTEGER*2 ICN(LICN), LENR(N), LENRL(N), LENOFF(N), IP(N), IQ(N)
4054 INTEGER ICN(LICN), LENR(N), LENRL(N), LENOFF(N), IP(N), IQ(N)
4055 Change
4056 C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON.
4057 COMMON /MA30HD/ RESID
4058 DATA ZERO /0.0D0/
4059 C 0000000
4060 C THE FINAL VALUE OF RESID IS THE MAXIMUM RESIDUAL FOR AN INCONSISTENT 0000000
4061 C SET OF EQUATIONS. 0000000
4062 RESID = ZERO 0000000
4063 C NOBLOC IS .TRUE. IF SUBROUTINE BLOCK HAS BEEN USED PREVIOUSLY AND 0000000
4064 C IS .FALSE. OTHERWISE. THE VALUE .FALSE. MEANS THAT LENOFF 0000000
4065 C WILL NOT BE SUBSEQUENTLY ACCESSED. 0000000
4066 NOBLOC = LENOFF(1).LT.0 0000000
4067 IF (MTYPE.NE.1) GO TO 140 0000000
4068 C 0000000
4069 C WE NOW SOLVE A * X = B. 0000000
4070 C NEG IS USED TO INDICATE WHEN THE LAST ROW IN A BLOCK HAS BEEN 0000000
4071 C REACHED. IT IS THEN SET TO TRUE WHEREAFTER BACKSUBSTITUTION IS 0000000
4072 C PERFORMED ON THE BLOCK. 0000000
4073 NEG = .FALSE. 0000000
4074 C IP(N) IS NEGATED SO THAT THE LAST ROW OF THE LAST BLOCK CAN BE 0000000
4075 C RECOGNISED. IT IS RESET TO ITS POSITIVE VALUE ON EXIT. 0000000
4076 IP(N) = -IP(N) 0000000
4077 C PREORDER VECTOR ... W(I) = X(IP(I)) 0000000
4078 DO 10 II=1,N 0000000
4079 I = IP(II) 0000000
4080 I = IABS(I) 0000000
4081 W(II) = X(I) 0000000
4082 10 CONTINUE 0000000
4083 C LT HOLDS THE POSITION OF THE FIRST NON-ZERO IN THE CURRENT ROW OF THE 0000000
4084 C OFF-DIAGONAL BLOCKS. 0000000
4085 LT = 1 0000000
4086 C IFIRST HOLDS THE INDEX OF THE FIRST ROW IN THE CURRENT BLOCK. 0000000
4087 IFIRST = 1 0000000
4088 C IBLOCK HOLDS THE POSITION OF THE FIRST NON-ZERO IN THE CURRENT ROW 0000000
4089 C OF THE LU DECOMPOSITION OF THE DIAGONAL BLOCKS. 0000000
4090 IBLOCK = IDISP(1) 0000000
4091 C IF I IS NOT THE LAST ROW OF A BLOCK, THEN A PASS THROUGH THIS LOOP 0000000
4092 C ADDS THE INNER PRODUCT OF ROW I OF THE OFF-DIAGONAL BLOCKS AND W 0000000
4093 C TO W AND PERFORMS FORWARD ELIMINATION USING ROW I OF THE LU 0000000
4094 C DECOMPOSITION. IF I IS THE LAST ROW OF A BLOCK THEN, AFTER 0000000
4095 C PERFORMING THESE AFOREMENTIONED OPERATIONS, BACKSUBSTITUTION IS 0000000
4096 C PERFORMED USING THE ROWS OF THE BLOCK. 0000000
4097 DO 120 I=1,N 0000000
4098 WI = W(I) 0000000
4099 IF (NOBLOC) GO TO 30 0000000
4100 IF (LENOFF(I).EQ.0) GO TO 30 0000000
4101 C OPERATIONS USING LOWER TRIANGULAR BLOCKS. 0000000
4102 C LTEND IS THE END OF ROW I IN THE OFF-DIAGONAL BLOCKS. 0000000
4103 LTEND = LT + LENOFF(I) - 1 0000000
4104 DO 20 JJ=LT,LTEND 0000000
4105 J = ICN(JJ) 0000000
4106 WI = WI - A(JJ)*W(J) 0000000
4107 20 CONTINUE 0000000
4108 C LT IS SET THE BEGINNING OF THE NEXT OFF-DIAGONAL ROW. 0000000
4109 LT = LTEND + 1 0000000
4110 C SET NEG TO .TRUE. IF WE ARE ON THE LAST ROW OF THE BLOCK. 0000000
4111 30 IF (IP(I).LT.0) NEG = .TRUE. 0000000
4112 IF (LENRL(I).EQ.0) GO TO 50 0000000
4113 C FORWARD ELIMINATION PHASE. 0000000
4114 C IEND IS THE END OF THE L PART OF ROW I IN THE LU DECOMPOSITION. 0000000
4115 IEND = IBLOCK + LENRL(I) - 1 0000000
4116 DO 40 JJ=IBLOCK,IEND 0000000
4117 J = ICN(JJ) 0000000
4118 WI = WI + A(JJ)*W(J) 0000000
4119 40 CONTINUE 0000000
4120 C IBLOCK IS ADJUSTED TO POINT TO THE START OF THE NEXT ROW. 0000000
4121 50 IBLOCK = IBLOCK + LENR(I) 0000000
4122 W(I) = WI 0000000
4123 IF (.NOT.NEG) GO TO 120 0000000
4124 C BACK SUBSTITUTION PHASE. 0000000
4125 C J1 IS POSITION IN A/ICN AFTER END OF BLOCK BEGINNING IN ROW IFIRST 0000000
4126 C AND ENDING IN ROW I. 0000000
4127 J1 = IBLOCK 0000000
4128 C ARE THERE ANY SINGULARITIES IN THIS BLOCK? IF NOT, CONTINUE WITH 0000000
4129 C THE BACKSUBSTITUTION. 0000000
4130 IB = I 0000000
4131 IF (IQ(I).GT.0) GO TO 70 0000000
4132 DO 60 III=IFIRST,I 0000000
4133 IB = I - III + IFIRST 0000000
4134 IF (IQ(IB).GT.0) GO TO 70 0000000
4135 J1 = J1 - LENR(IB) 0000000
4136 RESID = DMAX1(RESID,DABS(W(IB)))
4137 W(IB) = ZERO 0000000
4138 60 CONTINUE 0000000
4139 C ENTIRE BLOCK IS SINGULAR. 0000000
4140 GO TO 110 0000000
4141 C EACH PASS THROUGH THIS LOOP PERFORMS THE BACK-SUBSTITUTION 0000000
4142 C OPERATIONS FOR A SINGLE ROW, STARTING AT THE END OF THE BLOCK AND 0000000
4143 C WORKING THROUGH IT IN REVERSE ORDER. 0000000
4144 70 DO 100 III=IFIRST,IB 0000000
4145 II = IB - III + IFIRST 0000000
4146 C J2 IS END OF ROW II. 0000000
4147 J2 = J1 - 1 0000000
4148 C J1 IS BEGINNING OF ROW II. 0000000
4149 J1 = J1 - LENR(II) 0000000
4150 C JPIV IS THE POSITION OF THE PIVOT IN ROW II. 0000000
4151 JPIV = J1 + LENRL(II) 0000000
4152 JPIVP1 = JPIV + 1 0000000
4153 C JUMP IF ROW II OF U HAS NO NON-ZEROS. 0000000
4154 IF (J2.LT.JPIVP1) GO TO 90 0000000
4155 WII = W(II) 0000000
4156 DO 80 JJ=JPIVP1,J2 0000000
4157 J = ICN(JJ) 0000000
4158 WII = WII - A(JJ)*W(J) 0000000
4159 80 CONTINUE 0000000
4160 W(II) = WII 0000000
4161 90 W(II) = W(II)/A(JPIV) 0000000
4162 100 CONTINUE 0000000
4163 110 IFIRST = I + 1 0000000
4164 NEG = .FALSE. 0000000
4165 120 CONTINUE 0000000
4166 C 0000000
4167 C REORDER SOLUTION VECTOR ... X(I) = W(IQINVERSE(I)) 0000000
4168 DO 130 II=1,N 0000000
4169 I = IQ(II) 0000000
4170 I = IABS(I) 0000000
4171 X(I) = W(II) 0000000
4172 130 CONTINUE 0000000
4173 IP(N) = -IP(N) 0000000
4174 GO TO 320 0000000
4175 C 0000000
4176 C 0000000
4177 C WE NOW SOLVE ATRANSPOSE * X = B. 0000000
4178 C PREORDER VECTOR ... W(I)=X(IQ(I)) 0000000
4179 140 DO 150 II=1,N 0000000
4180 I = IQ(II) 0000000
4181 I = IABS(I) 0000000
4182 W(II) = X(I) 0000000
4183 150 CONTINUE 0000000
4184 C LJ1 POINTS TO THE BEGINNING THE CURRENT ROW IN THE OFF-DIAGONAL 0000000
4185 C BLOCKS. 0000000
4186 LJ1 = IDISP(1) 0000000
4187 C IBLOCK IS INITIALIZED TO POINT TO THE BEGINNING OF THE BLOCK AFTER 0000000
4188 C THE LAST ONE ] 0000000
4189 IBLOCK = IDISP(2) + 1 0000000
4190 C ILAST IS THE LAST ROW IN THE CURRENT BLOCK. 0000000
4191 ILAST = N 0000000
4192 C IBLEND POINTS TO THE POSITION AFTER THE LAST NON-ZERO IN THE 0000000
4193 C CURRENT BLOCK. 0000000
4194 IBLEND = IBLOCK 0000000
4195 C EACH PASS THROUGH THIS LOOP OPERATES WITH ONE DIAGONAL BLOCK AND 0000000
4196 C THE OFF-DIAGONAL PART OF THE MATRIX CORRESPONDING TO THE ROWS 0000000
4197 C OF THIS BLOCK. THE BLOCKS ARE TAKEN IN REVERSE ORDER AND THE 0000000
4198 C NUMBER OF TIMES THE LOOP IS ENTERED IS MIN(N,NO. BLOCKS+1). 0000000
4199 DO 290 NUMBLK=1,N 0000000
4200 IF (ILAST.EQ.0) GO TO 300 0000000
4201 IBLOCK = IBLOCK - LENR(ILAST) 0000000
4202 C THIS LOOP FINDS THE INDEX OF THE FIRST ROW IN THE CURRENT BLOCK.. 0000000
4203 C IT IS FIRST AND IBLOCK IS SET TO THE POSITION OF THE BEGINNING 0000000
4204 C OF THIS FIRST ROW. 0000000
4205 DO 160 K=1,N 0000000
4206 II = ILAST - K 0000000
4207 IF (II.EQ.0) GO TO 170 0000000
4208 IF (IP(II).LT.0) GO TO 170 0000000
4209 IBLOCK = IBLOCK - LENR(II) 0000000
4210 160 CONTINUE 0000000
4211 170 IFIRST = II + 1 0000000
4212 C J1 POINTS TO THE POSITION OF THE BEGINNING OF ROW I (LT PART) OR PIVOT0000000
4213 J1 = IBLOCK 0000000
4214 C FORWARD ELIMINATION. 0000000
4215 C EACH PASS THROUGH THIS LOOP PERFORMS THE OPERATIONS FOR ONE ROW OF THE0000000
4216 C BLOCK. IF THE CORRESPONDING ENTRY OF W IS ZERO THEN THE 0000000
4217 C OPERATIONS CAN BE AVOIDED. 0000000
4218 DO 210 I=IFIRST,ILAST 0000000
4219 IF (W(I).EQ.ZERO) GO TO 200 0000000
4220 C JUMP IF ROW I SINGULAR. 0000000
4221 IF (IQ(I).LT.0) GO TO 220 0000000
4222 C J2 FIRST POINTS TO THE PIVOT IN ROW I AND THEN IS MADE TO POINT TO THE0000000
4223 C FIRST NON-ZERO IN THE U TRANSPOSE PART OF THE ROW. 0000000
4224 J2 = J1 + LENRL(I) 0000000
4225 WI = W(I)/A(J2) 0000000
4226 IF (LENR(I)-LENRL(I).EQ.1) GO TO 190 0000000
4227 J2 = J2 + 1 0000000
4228 C J3 POINTS TO THE END OF ROW I. 0000000
4229 J3 = J1 + LENR(I) - 1 0000000
4230 DO 180 JJ=J2,J3 0000000
4231 J = ICN(JJ) 0000000
4232 W(J) = W(J) - A(JJ)*WI 0000000
4233 180 CONTINUE 0000000
4234 190 W(I) = WI 0000000
4235 200 J1 = J1 + LENR(I) 0000000
4236 210 CONTINUE 0000000
4237 GO TO 240 0000000
4238 C DEALS WITH REST OF BLOCK WHICH IS SINGULAR. 0000000
4239 220 DO 230 II=I,ILAST 0000000
4240 RESID = DMAX1(RESID,DABS(W(II)))
4241 W(II) = ZERO 0000000
4242 230 CONTINUE 0000000
4243 C BACK SUBSTITUTION. 0000000
4244 C THIS LOOP DOES THE BACK SUBSTITUTION ON THE ROWS OF THE BLOCK IN 0000000
4245 C THE REVERSE ORDER DOING IT SIMULTANEOUSLY ON THE L TRANSPOSE PART 0000000
4246 C OF THE DIAGONAL BLOCKS AND THE OFF-DIAGONAL BLOCKS. 0000000
4247 240 J1 = IBLEND 0000000
4248 DO 280 IBACK=IFIRST,ILAST 0000000
4249 I = ILAST - IBACK + IFIRST 0000000
4250 C J1 POINTS TO THE BEGINNING OF ROW I. 0000000
4251 J1 = J1 - LENR(I) 0000000
4252 IF (LENRL(I).EQ.0) GO TO 260 0000000
4253 C J2 POINTS TO THE END OF THE L TRANSPOSE PART OF ROW I. 0000000
4254 J2 = J1 + LENRL(I) - 1 0000000
4255 DO 250 JJ=J1,J2 0000000
4256 J = ICN(JJ) 0000000
4257 W(J) = W(J) + A(JJ)*W(I) 0000000
4258 250 CONTINUE 0000000
4259 260 IF (NOBLOC) GO TO 280 0000000
4260 C OPERATIONS USING LOWER TRIANGULAR BLOCKS. 0000000
4261 IF (LENOFF(I).EQ.0) GO TO 280 0000000
4262 C LJ2 POINTS TO THE END OF ROW I OF THE OFF-DIAGONAL BLOCKS. 0000000
4263 LJ2 = LJ1 - 1 0000000
4264 C LJ1 POINTS TO THE BEGINNING OF ROW I OF THE OFF-DIAGONAL BLOCKS. 0000000
4265 LJ1 = LJ1 - LENOFF(I) 0000000
4266 DO 270 JJ=LJ1,LJ2 0000000
4267 J = ICN(JJ) 0000000
4268 W(J) = W(J) - A(JJ)*W(I) 0000000
4269 270 CONTINUE 0000000
4270 280 CONTINUE 0000000
4271 IBLEND = J1 0000000
4272 ILAST = IFIRST - 1 0000000
4273 290 CONTINUE 0000000
4274 C REORDER SOLUTION VECTOR ... X(I)=W(IPINVERSE(I)) 0000000
4275 300 DO 310 II=1,N 0000000
4276 I = IP(II) 0000000
4277 I = IABS(I) 0000000
4278 X(I) = W(II) 0000000
4279 310 CONTINUE 0000000
4280 C 0000000
4281 320 RETURN 0000000
4282 END 0000000
4283 BLOCK DATA ma30$data
4284 C_270390ak BLOCK DATA 0000000
4285 C ALTHOUGH ALL COMMON BLOCK VARIABLES DO NOT HAVE DEFAULT VALUES, 0000000
4286 C WE COMMENT ON ALL THE COMMON BLOCK VARIABLES HERE. 0000000
4287 C 0000000
4288 C COMMON BLOCK MA30E/ED HOLDS CONTROL PARAMETERS .... 0000000
4289 C COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3 0000000
4290 C THE INTEGER LP IS THE UNIT NUMBER TO WHICH THE ERROR MESSAGES ARE 0000000
4291 C SENT. LP HAS A DEFAULT VALUE OF 6. THIS DEFAULT VALUE CAN BE 0000000
4292 C RESET BY THE USER, IF DESIRED. A VALUE OF 0 SUPPRESSES ALL 0000000
4293 C MESSAGES. 0000000
4294 C THE LOGICAL VARIABLES ABORT1,ABORT2,ABORT3 ARE USED TO CONTROL THE 0000000
4295 C CONDITIONS UNDER WHICH THE SUBROUTINE WILL TERMINATE. 0000000
4296 C IF ABORT1 IS .TRUE. THEN THE SUBROUTINE WILL EXIT IMMEDIATELY ON 0000000
4297 C DETECTING STRUCTURAL SINGULARITY. 0000000
4298 C IF ABORT2 IS .TRUE. THEN THE SUBROUTINE WILL EXIT IMMEDIATELY ON 0000000
4299 C DETECTING NUMERICAL SINGULARITY. 0000000
4300 C IF ABORT3 IS .TRUE. THEN THE SUBROUTINE WILL EXIT IMMEDIATELY WHEN 0000000
4301 C THE AVAILABLE SPACE IN A/ICN IS FILLED UP BY THE PREVIOUSLY 0000000
4302 C DECOMPOSED, ACTIVE, AND UNDECOMPOSED PARTS OF THE MATRIX. 0000000
4303 C THE DEFAULT VALUES FOR ABORT1,ABORT2,ABORT3 ARE SET TO .TRUE.,.TRUE. 0000000
4304 C AND .FALSE. RESPECTIVELY. 0000000
4305 C 0000000
4306 C THE VARIABLES IN THE COMMON BLOCK MA30F/FD ARE USED TO PROVIDE THE 0000000
4307 C USER WITH INFORMATION ON THE DECOMPOSITION. 0000000
4308 C COMMON /MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN 0000000
4309 C IRNCP AND ICNCP ARE INTEGER VARIABLES USED TO MONITOR THE ADEQUACY 0000000
4310 C OF THE ALLOCATED SPACE IN ARRAYS IRN AND A/ICN RESPECTIVELY, BY 0000000
4311 C TAKING ACCOUNT OF THE NUMBER OF DATA MANAGEMENT COMPRESSES 0000000
4312 C REQUIRED ON THESE ARRAYS. IF IRNCP OR ICNCP IS FAIRLY LARGE (SAY 0000000
4313 C GREATER THAN N/10), IT MAY BE ADVANTAGEOUS TO INCREASE THE SIZE 0000000
4314 C OF THE CORRESPONDING ARRAY(S). IRNCP AND ICNCP ARE INITIALIZED 0000000
4315 C TO ZERO ON ENTRY TO MA30A/AD AND ARE INCREMENTED EACH TIME THE 0000000
4316 C COMPRESSING ROUTINE MA30D/DD IS ENTERED. 0000000
4317 C ICNCP IS THE NUMBER OF COMPRESSES ON A/ICN. 0000000
4318 C IRNCP IS THE NUMBER OF COMPRESSES ON IRN. 0000000
4319 C IRANK IS AN INTEGER VARIABLE WHICH GIVES AN ESTIMATE (ACTUALLY AN 0000000
4320 C UPPER BOUND) OF THE RANK OF THE MATRIX. ON AN EXIT WITH IFLAG 0000000
4321 C EQUAL TO 0, THIS WILL BE EQUAL TO N. 0000000
4322 C MINIRN IS AN INTEGER VARIABLE WHICH, AFTER A SUCCESSFUL CALL TO 0000000
4323 C MA30A/AD, INDICATES THE MINIMUM LENGTH TO WHICH IRN CAN BE 0000000
4324 C REDUCED WHILE STILL PERMITTING A SUCCESSFUL DECOMPOSITION OF THE 0000000
4325 C SAME MATRIX. IF, HOWEVER, THE USER WERE TO DECREASE THE LENGTH 0000000
4326 C OF IRN TO THAT SIZE, THE NUMBER OF COMPRESSES (IRNCP) MAY BE 0000000
4327 C VERY HIGH AND QUITE COSTLY. IF LIRN IS NOT LARGE ENOUGH TO BEGIN 0000000
4328 C THE DECOMPOSITION ON A DIAGONAL BLOCK, MINIRN WILL BE EQUAL TO 0000000
4329 C THE VALUE REQUIRED TO CONTINUE THE DECOMPOSITION AND IFLAG WILL 0000000
4330 C BE SET TO -3 OR -6. A VALUE OF LIRN SLIGHTLY GREATER THAN THIS 0000000
4331 C (SAY ABOUT N/2) WILL USUALLY PROVIDE ENOUGH SPACE TO COMPLETE 0000000
4332 C THE DECOMPOSITION ON THAT BLOCK. IN THE EVENT OF ANY OTHER 0000000
4333 C FAILURE MINIRN GIVES THE MINIMUM SIZE OF IRN REQUIRED FOR A 0000000
4334 C SUCCESSFUL DECOMPOSITION UP TO THAT POINT. 0000000
4335 C MINICN IS AN INTEGER VARIABLE WHICH AFTER A SUCCESSFUL CALL TO 0000000
4336 C MA30A/AD, INDICATES THE MINIMUM SIZE OF LICN REQUIRED TO ENABLE 0000000
4337 C A SUCCESSFUL DECOMPOSITION. IN THE EVENT OF FAILURE WITH IFLAG= 0000000
4338 C -5, MINICN WILL, IF ABORT3 IS LEFT SET TO .FALSE., INDICATE THE 0000000
4339 C MINIMUM LENGTH THAT WOULD BE SUFFICIENT TO PREVENT THIS ERROR IN 0000000
4340 C A SUBSEQUENT RUN ON AN IDENTICAL MATRIX. AGAIN THE USER MAY 0000000
4341 C PREFER TO USE A VALUE OF ICN SLIGHTLY GREATER THAN MINICN FOR 0000000
4342 C SUBSEQUENT RUNS TO AVOID TOO MANY CONPRESSES (ICNCP). IN THE 0000000
4343 C EVENT OF FAILURE WITH IFLAG EQUAL TO ANY NEGATIVE VALUE EXCEPT 0000000
4344 C -4, MINICN WILL GIVE THE MINIMUM LENGTH TO WHICH LICN COULD BE 0000000
4345 C REDUCED TO ENABLE A SUCCESSFUL DECOMPOSITION TO THE POINT AT 0000000
4346 C WHICH FAILURE OCCURRED. NOTICE THAT, ON A SUCCESSFUL ENTRY 0000000
4347 C IDISP(2) GIVES THE AMOUNT OF SPACE IN A/ICN REQUIRED FOR THE 0000000
4348 C DECOMPOSITION WHILE MINICN WILL USUALLY BE SLIGHTLY GREATER 0000000
4349 C BECAUSE OF THE NEED FOR "ELBOW ROOM". IF THE USER IS VERY 0000000
4350 C UNSURE HOW LARGE TO MAKE LICN, THE VARIABLE MINICN CAN BE USED 0000000
4351 C TO PROVIDE THAT INFORMATION. A PRELIMINARY RUN SHOULD BE 0000000
4352 C PERFORMED WITH ABORT3 LEFT SET TO .FALSE. AND LICN ABOUT 3/2 0000000
4353 C TIMES AS BIG AS THE NUMBER OF NON-ZEROS IN THE ORIGINAL MATRIX. 0000000
4354 C UNLESS THE INITIAL PROBLEM IS VERY SPARSE (WHEN THE RUN WILL BE 0000000
4355 C SUCCESSFUL) OR FILLS IN EXTREMELY BADLY (GIVING AN ERROR RETURN 0000000
4356 C WITH IFLAG EQUAL TO -4), AN ERROR RETURN WITH IFLAG EQUAL TO -5 0000000
4357 C SHOULD RESULT AND MINICN WILL GIVE THE AMOUNT OF SPACE REQUIRED 0000000
4358 C FOR A SUCCESSFUL DECOMPOSITION. 0000000
4359 C 0000000
4360 C COMMON BLOCK MA30G/GD IS USED BY THE MA30B/BD ENTRY ONLY. 0000000
4361 C COMMON /MA30GD/ EPS, RMIN 0000000
4362 C EPS IS A REAL/DOUBLE PRECISION VARIABLE. IT IS USED TO TEST FOR 0000000
4363 C SMALL PIVOTS. ITS DEFAULT VALUE IS 1.0E-4 (1.0D-4 IN D VERSION). 0000000
4364 C IF THE USER SETS EPS TO ANY VALUE GREATER THAN 1.0, THEN NO 0000000
4365 C CHECK IS MADE ON THE SIZE OF THE PIVOTS. ALTHOUGH THE ABSENCE OF 0000000
4366 C SUCH A CHECK WOULD FAIL TO WARN THE USER OF BAD INSTABILITY, ITS 0000000
4367 C ABSENCE WILL ENABLE MA30B/BD TO RUN SLIGHTLY FASTER. AN A 0000000
4368 C POSTERIORI CHECK ON THE STABILITY OF THE FACTORIZATION CAN BE 0000000
4369 C OBTAINED FROM MC24A/AD. 0000000
4370 C RMIN IS A REAL/DOUBLE PRECISION VARIABLE WHICH GIVES THE USER SOME 0000000
4371 C INFORMATION ABOUT THE STABILITY OF THE DECOMPOSITION. AT EACH 0000000
4372 C STAGE OF THE LU DECOMPOSITION THE MAGNITUDE OF THE PIVOT APIV 0000000
4373 C IS COMPARED WITH THE LARGEST OFF-DIAGONAL ENTRY CURRENTLY IN ITS 0000000
4374 C ROW (ROW OF U), ROWMAX SAY. IF THE RATIO 0000000
4375 C MIN (APIV/ROWMAX) 0000000
4376 C WHERE THE MINIMUM IS TAKEN OVER ALL THE ROWS, IS LESS THAN EPS 0000000
4377 C THEN RMIN IS SET TO THIS MINIMUM VALUE AND IFLAG IS RETURNED 0000000
4378 C WITH THE VALUE +I WHERE I IS THE ROW IN WHICH THIS MINIMUM 0000000
4379 C OCCURS. IF THE USER SETS EPS GREATER THAN ONE, THEN THIS TEST 0000000
4380 C IS NOT PERFORMED. IN THIS CASE, AND WHEN THERE ARE NO SMALL 0000000
4381 C PIVOTS RMIN WILL BE SET EQUAL TO EPS. 0000000
4382 C 0000000
4383 C COMMON BLOCK MA30H/HD IS USED BY MA30C/CD ONLY. 0000000
4384 C COMMON /MA30HD/ RESID 0000000
4385 C RESID IS A REAL/DOUBLE PRECISION VARIABLE. IN THE CASE OF SINGULAR 0000000
4386 C OR RECTANGULAR MATRICES ITS FINAL VALUE WILL BE EQUAL TO THE 0000000
4387 C MAXIMUM RESIDUAL FOR THE UNSATISFIED EQUATIONS; OTHERWISE ITS 0000000
4388 C VALUE WILL BE SET TO ZERO. 0000000
4389 C 0000000
4390 C COMMON BLOCK MA30I/ID CONTROLS THE USE OF DROP TOLERANCES, THE 0000000
4391 C MODIFIED PIVOT OPTION AND THE THE CALCULATION OF THE LARGEST 0000000
4392 C ENTRY IN THE FACTORIZATION PROCESS. THIS COMMON BLOCK WAS ADDED 0000000
4393 C TO THE MA30 PACKAGE IN FEBRUARY, 1983. 0000000
4394 C COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG 0000000
4395 C TOL IS A REAL/DOUBLE PRECISION VARIABLE. IF IT IS SET TO A POSITIVE 0000000
4396 C VALUE, THEN MA30A/AD WILL DROP FROM THE FACTORS ANY NON-ZERO 0000000
4397 C WHOSE MODULUS IS LESS THAN TOL. THE FACTORIZATION WILL THEN 0000000
4398 C REQUIRE LESS STORAGE BUT WILL BE INACCURATE. AFTER A RUN OF 0000000
4399 C MA30A/AD WHERE ENTRIES HAVE BEEN DROPPED, MA30B/BD SHOULD NOT 0000000
4400 C BE CALLED. THE DEFAULT VALUE FOR TOL IS 0.0. 0000000
4401 C BIG IS A REAL/DOUBLE PRECISION VARIABLE. IF LBIG HAS BEEN SET TO 0000000
4402 C .TRUE., BIG WILL BE SET TO THE LARGEST ENTRY ENCOUNTERED DURING 0000000
4403 C THE FACTORIZATION. 0000000
4404 C NDROP IS AN INTEGER VARIABLE. IF TOL HAS BEEN SET POSITIVE, ON EXIT 0000000
4405 C FROM MA30A/AD, NDROP WILL HOLD THE NUMBER OF ENTRIES DROPPED 0000000
4406 C FROM THE DATA STRUCTURE. 0000000
4407 C NSRCH IS AN INTEGER VARIABLE. IF NSRCH IS SET TO A VALUE LESS THAN 0000000
4408 C OR EQUAL TO N, THEN A DIFFERENT PIVOT OPTION WILL BE EMPLOYED BY 0000000
4409 C MA30A/AD. THIS MAY RESULT IN DIFFERENT FILL-IN AND EXECUTION 0000000
4410 C TIME FOR MA30A/AD. IF NSRCH IS LESS THAN OR EQUAL TO N, THE 0000000
4411 C WORKSPACE ARRAYS LASTC AND NEXTC ARE NOT REFERENCED BY MA30A/AD. 0000000
4412 C THE DEFAULT VALUE FOR NSRCH IS 32768. 0000000
4413 C LBIG IS A LOGICAL VARIABLE. IF LBIG IS SET TO .TRUE., THE VALUE OF 0000000
4414 C THE LARGEST ENTRY ENCOUNTERED IN THE FACTORIZATION BY MA30A/AD 0000000
4415 C IS RETURNED IN BIG. SETTING LBIG TO .TRUE. WILL MARGINALLY 0000000
4416 C INCREASE THE FACTORIZATION TIME FOR MA30A/AD AND WILL INCREASE 0000000
4417 C THAT FOR MA30B/BD BY ABOUT 20%. THE DEFAULT VALUE FOR LBIG IS 0000000
4418 C .FALSE. 0000000
4419 C 0000000
4420 DOUBLE PRECISION EPS, RMIN, TOL, BIG
4421 LOGICAL ABORT1, ABORT2, ABORT3, LBIG
4422 COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
4423 COMMON /MA30GD/ EPS, RMIN
4424 COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
4425 DATA EPS /1.0D-4/, TOL /0.0D0/, BIG /0.0D0/
4426 DATA LP /6/, NSRCH /32768/
4427 DATA LBIG /.FALSE./
4428 DATA ABORT1 /.TRUE./, ABORT2 /.TRUE./, ABORT3 /.FALSE./
4429 END
4430 SUBROUTINE XERRWV (MSG, NMES, NERR, IERT, NI, I1, I2, NR, R1, R2)
4431 INTEGER MSG, NMES, NERR, IERT, NI, I1, I2, NR,
4432 1 I, LUN, LUNIT, MESFLG, NCPW, NCH, NWDS
4433 DOUBLE PRECISION R1, R2
4434 DIMENSION MSG(NMES)
4435 C-----------------------------------------------------------------------
4436 C SUBROUTINES XERRWV, XSETF, AND XSETUN, AS GIVEN HERE, CONSTITUTE
4437 C A SIMPLIFIED VERSION OF THE SLATEC ERROR HANDLING PACKAGE.
4438 C WRITTEN BY A. C. HINDMARSH AT LLNL. VERSION OF AUGUST 13, 1981.
4439 C THIS VERSION IS IN DOUBLE PRECISION.
4440 C
4441 C ALL ARGUMENTS ARE INPUT ARGUMENTS.
4442 C
4443 C MSG = THE MESSAGE (HOLLERITH LITTERAL OR INTEGER ARRAY).
4444 C NMES = THE LENGTH OF MSG (NUMBER OF CHARACTERS).
4445 C NERR = THE ERROR NUMBER (NOT USED).
4446 C IERT = THE ERROR TYPE..
4447 C 1 MEANS RECOVERABLE (CONTROL RETURNS TO CALLER).
4448 C 2 MEANS FATAL (RUN IS ABORTED--SEE NOTE BELOW).
4449 C NI = NUMBER OF INTEGERS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE.
4450 C I1,I2 = INTEGERS TO BE PRINTED, DEPENDING ON NI.
4451 C NR = NUMBER OF REALS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE.
4452 C R1,R2 = REALS TO BE PRINTED, DEPENDING ON NR.
4453 C
4454 C NOTE.. THIS ROUTINE IS MACHINE-DEPENDENT AND SPECIALIZED FOR USE
4455 C IN LIMITED CONTEXT, IN THE FOLLOWING WAYS..
4456 C 1. THE NUMBER OF HOLLERITH CHARACTERS STORED PER WORD, DENOTED
4457 C BY NCPW BELOW, IS A DATA-LOADED CONSTANT.
4458 C 2. THE VALUE OF NMES IS ASSUMED TO BE AT MOST 60.
4459 C (MULTI-LINE MESSAGES ARE GENERATED BY REPEATED CALLS.)
4460 C 3. IF IERT = 2, CONTROL PASSES TO THE STATEMENT STOP
4461 C TO ABORT THE RUN. THIS STATEMENT MAY BE MACHINE-DEPENDENT.
4462 C 4. R1 AND R2 ARE ASSUMED TO BE IN DOUBLE PRECISION AND ARE PRINTED
4463 C IN D21.13 FORMAT.
4464 C 5. THE COMMON BLOCK /EH0001/ BELOW IS DATA-LOADED (A MACHINE-
4465 C DEPENDENT FEATURE) WITH DEFAULT VALUES.
4466 C THIS BLOCK IS NEEDED FOR PROPER RETENTION OF PARAMETERS USED BY
4467 C THIS ROUTINE WHICH THE USER CAN RESET BY CALLING XSETF OR XSETUN.
4468 C THE VARIABLES IN THIS BLOCK ARE AS FOLLOWS..
4469 C MESFLG = PRINT CONTROL FLAG..
4470 C 1 MEANS PRINT ALL MESSAGES (THE DEFAULT).
4471 C 0 MEANS NO PRINTING.
4472 C LUNIT = LOGICAL UNIT NUMBER FOR MESSAGES.
4473 C THE DEFAULT IS 6 (MACHINE-DEPENDENT).
4474 C-----------------------------------------------------------------------
4475 C THE FOLLOWING ARE INSTRUCTIONS FOR INSTALLING THIS ROUTINE
4476 C IN DIFFERENT MACHINE ENVIRONMENTS.
4477 C
4478 C TO CHANGE THE DEFAULT OUTPUT UNIT, CHANGE THE DATA STATEMENT
4479 C IN THE BLOCK DATA SUBPROGRAM BELOW.
4480 C
4481 C FOR A DIFFERENT NUMBER OF CHARACTERS PER WORD, CHANGE THE
4482 C DATA STATEMENT SETTING NCPW BELOW, AND FORMAT 10. ALTERNATIVES FOR
4483 C VARIOUS COMPUTERS ARE SHOWN IN COMMENT CARDS.
4484 C
4485 C FOR A DIFFERENT RUN-ABORT COMMAND, CHANGE THE STATEMENT FOLLOWING
4486 C STATEMENT 100 AT THE END.
4487 C-----------------------------------------------------------------------
4488 COMMON /EH0001/ MESFLG, LUNIT
4489 C-----------------------------------------------------------------------
4490 C THE FOLLOWING DATA-LOADED VALUE OF NCPW IS VALID FOR THE CDC-6600
4491 C AND CDC-7600 COMPUTERS.
4492 C DATA NCPW/10/
4493 C THE FOLLOWING IS VALID FOR THE CRAY-1 COMPUTER.
4494 C DATA NCPW/8/
4495 C THE FOLLOWING IS VALID FOR THE BURROUGHS 6700 AND 7800 COMPUTERS.
4496 C DATA NCPW/6/
4497 C THE FOLLOWING IS VALID FOR THE PDP-10 COMPUTER.
4498 C DATA NCPW/5/
4499 C THE FOLLOWING IS VALID FOR THE VAX COMPUTER WITH 4 BYTES PER INTEGER,
4500 C AND FOR THE IBM-360, IBM-370, IBM-303X, AND IBM-43XX COMPUTERS.
4501 DATA NCPW/4/
4502 C THE FOLLOWING IS VALID FOR THE PDP-11, OR VAX WITH 2-BYTE INTEGERS.
4503 C DATA NCPW/2/
4504 C-----------------------------------------------------------------------
4505 IF (MESFLG .EQ. 0) GO TO 100
4506 C GET LOGICAL UNIT NUMBER. ---------------------------------------------
4507 LUN = LUNIT
4508 C GET NUMBER OF WORDS IN MESSAGE. --------------------------------------
4509 NCH = MIN0(NMES,60)
4510 NWDS = NCH/NCPW
4511 IF (NCH .NE. NWDS*NCPW) NWDS = NWDS + 1
4512 C WRITE THE MESSAGE. ---------------------------------------------------
4513 WRITE (LUN, 10) (MSG(I),I=1,NWDS)
4514 C-----------------------------------------------------------------------
4515 C THE FOLLOWING FORMAT STATEMENT IS TO HAVE THE FORM
4516 C 10 FORMAT(1X,MMANN)
4517 C WHERE NN = NCPW AND MM IS THE SMALLEST INTEGER .GE. 60/NCPW.
4518 C THE FOLLOWING IS VALID FOR NCPW = 10.
4519 C 10 FORMAT(1X,6A10)
4520 C THE FOLLOWING IS VALID FOR NCPW = 8.
4521 C 10 FORMAT(1X,8A8)
4522 C THE FOLLOWING IS VALID FOR NCPW = 6.
4523 C 10 FORMAT(1X,10A6)
4524 C THE FOLLOWING IS VALID FOR NCPW = 5.
4525 C 10 FORMAT(1X,12A5)
4526 C THE FOLLOWING IS VALID FOR NCPW = 4.
4527 10 FORMAT(1X,15A4)
4528 C THE FOLLOWING IS VALID FOR NCPW = 2.
4529 C 10 FORMAT(1X,30A2)
4530 C-----------------------------------------------------------------------
4531 IF (NI .EQ. 1) WRITE (LUN, 20) I1
4532 20 FORMAT(6X,23HIN ABOVE MESSAGE, I1 =,I10)
4533 IF (NI .EQ. 2) WRITE (LUN, 30) I1,I2
4534 30 FORMAT(6X,23HIN ABOVE MESSAGE, I1 =,I10,3X,4HI2 =,I10)
4535 IF (NR .EQ. 1) WRITE (LUN, 40) R1
4536 40 FORMAT(6X,23HIN ABOVE MESSAGE, R1 =,D21.13)
4537 IF (NR .EQ. 2) WRITE (LUN, 50) R1,R2
4538 50 FORMAT(6X,15HIN ABOVE, R1 =,D21.13,3X,4HR2 =,D21.13)
4539 C ABORT THE RUN IF IERT = 2. -------------------------------------------
4540 100 IF (IERT .NE. 2) RETURN
4541 STOP
4542 C----------------------- END OF SUBROUTINE XERRWV ----------------------
4543 END
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22