MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
opkda2.f
1 *DECK DGEFA
2  SUBROUTINE dgefa (A, LDA, N, IPVT, INFO)
3 C***BEGIN PROLOGUE DGEFA
4 C***PURPOSE Factor a matrix using Gaussian elimination.
5 C***CATEGORY D2A1
6 C***TYPE DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C)
7 C***KEYWORDS GENERAL MATRIX, LINEAR ALGEBRA, LINPACK,
8 C MATRIX FACTORIZATION
9 C***AUTHOR Moler, C. B., (U. of New Mexico)
10 C***DESCRIPTION
11 C
12 C DGEFA factors a double precision matrix by Gaussian elimination.
13 C
14 C DGEFA is usually called by DGECO, but it can be called
15 C directly with a saving in time if RCOND is not needed.
16 C (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .
17 C
18 C On Entry
19 C
20 C A DOUBLE PRECISION(LDA, N)
21 C the matrix to be factored.
22 C
23 C LDA INTEGER
24 C the leading dimension of the array A .
25 C
26 C N INTEGER
27 C the order of the matrix A .
28 C
29 C On Return
30 C
31 C A an upper triangular matrix and the multipliers
32 C which were used to obtain it.
33 C The factorization can be written A = L*U where
34 C L is a product of permutation and unit lower
35 C triangular matrices and U is upper triangular.
36 C
37 C IPVT INTEGER(N)
38 C an integer vector of pivot indices.
39 C
40 C INFO INTEGER
41 C = 0 normal value.
42 C = K if U(K,K) .EQ. 0.0 . This is not an error
43 C condition for this subroutine, but it does
44 C indicate that DGESL or DGEDI will divide by zero
45 C if called. Use RCOND in DGECO for a reliable
46 C indication of singularity.
47 C
48 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
49 C Stewart, LINPACK Users' Guide, SIAM, 1979.
50 C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX
51 C***REVISION HISTORY (YYMMDD)
52 C 780814 DATE WRITTEN
53 C 890831 Modified array declarations. (WRB)
54 C 890831 REVISION DATE from Version 3.2
55 C 891214 Prologue converted to Version 4.0 format. (BAB)
56 C 900326 Removed duplicate information from DESCRIPTION section.
57 C (WRB)
58 C 920501 Reformatted the REFERENCES section. (WRB)
59 C***END PROLOGUE DGEFA
60  INTEGER lda,n,ipvt(*),info
61  DOUBLE PRECISION a(lda,*)
62 C
63  DOUBLE PRECISION t
64  INTEGER idamax,j,k,kp1,l,nm1
65 C
66 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
67 C
68 C***FIRST EXECUTABLE STATEMENT DGEFA
69  info = 0
70  nm1 = n - 1
71  IF (nm1 .LT. 1) GO TO 70
72  DO 60 k = 1, nm1
73  kp1 = k + 1
74 C
75 C FIND L = PIVOT INDEX
76 C
77  l = idamax(n-k+1,a(k,k),1) + k - 1
78  ipvt(k) = l
79 C
80 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
81 C
82  IF (a(l,k) .EQ. 0.0d0) GO TO 40
83 C
84 C INTERCHANGE IF NECESSARY
85 C
86  IF (l .EQ. k) GO TO 10
87  t = a(l,k)
88  a(l,k) = a(k,k)
89  a(k,k) = t
90  10 CONTINUE
91 C
92 C COMPUTE MULTIPLIERS
93 C
94  t = -1.0d0/a(k,k)
95  CALL dscal(n-k,t,a(k+1,k),1)
96 C
97 C ROW ELIMINATION WITH COLUMN INDEXING
98 C
99  DO 30 j = kp1, n
100  t = a(l,j)
101  IF (l .EQ. k) GO TO 20
102  a(l,j) = a(k,j)
103  a(k,j) = t
104  20 CONTINUE
105  CALL daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
106  30 CONTINUE
107  GO TO 50
108  40 CONTINUE
109  info = k
110  50 CONTINUE
111  60 CONTINUE
112  70 CONTINUE
113  ipvt(n) = n
114  IF (a(n,n) .EQ. 0.0d0) info = n
115  RETURN
116  END
117 *DECK DGESL
118  SUBROUTINE dgesl (A, LDA, N, IPVT, B, JOB)
119 C***BEGIN PROLOGUE DGESL
120 C***PURPOSE Solve the real system A*X=B or TRANS(A)*X=B using the
121 C factors computed by DGECO or DGEFA.
122 C***CATEGORY D2A1
123 C***TYPE DOUBLE PRECISION (SGESL-S, DGESL-D, CGESL-C)
124 C***KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
125 C***AUTHOR Moler, C. B., (U. of New Mexico)
126 C***DESCRIPTION
127 C
128 C DGESL solves the double precision system
129 C A * X = B or TRANS(A) * X = B
130 C using the factors computed by DGECO or DGEFA.
131 C
132 C On Entry
133 C
134 C A DOUBLE PRECISION(LDA, N)
135 C the output from DGECO or DGEFA.
136 C
137 C LDA INTEGER
138 C the leading dimension of the array A .
139 C
140 C N INTEGER
141 C the order of the matrix A .
142 C
143 C IPVT INTEGER(N)
144 C the pivot vector from DGECO or DGEFA.
145 C
146 C B DOUBLE PRECISION(N)
147 C the right hand side vector.
148 C
149 C JOB INTEGER
150 C = 0 to solve A*X = B ,
151 C = nonzero to solve TRANS(A)*X = B where
152 C TRANS(A) is the transpose.
153 C
154 C On Return
155 C
156 C B the solution vector X .
157 C
158 C Error Condition
159 C
160 C A division by zero will occur if the input factor contains a
161 C zero on the diagonal. Technically this indicates singularity
162 C but it is often caused by improper arguments or improper
163 C setting of LDA . It will not occur if the subroutines are
164 C called correctly and if DGECO has set RCOND .GT. 0.0
165 C or DGEFA has set INFO .EQ. 0 .
166 C
167 C To compute INVERSE(A) * C where C is a matrix
168 C with P columns
169 C CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
170 C IF (RCOND is too small) GO TO ...
171 C DO 10 J = 1, P
172 C CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
173 C 10 CONTINUE
174 C
175 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
176 C Stewart, LINPACK Users' Guide, SIAM, 1979.
177 C***ROUTINES CALLED DAXPY, DDOT
178 C***REVISION HISTORY (YYMMDD)
179 C 780814 DATE WRITTEN
180 C 890831 Modified array declarations. (WRB)
181 C 890831 REVISION DATE from Version 3.2
182 C 891214 Prologue converted to Version 4.0 format. (BAB)
183 C 900326 Removed duplicate information from DESCRIPTION section.
184 C (WRB)
185 C 920501 Reformatted the REFERENCES section. (WRB)
186 C***END PROLOGUE DGESL
187  INTEGER lda,n,ipvt(*),job
188  DOUBLE PRECISION a(lda,*),b(*)
189 C
190  DOUBLE PRECISION ddot,t
191  INTEGER k,kb,l,nm1
192 C***FIRST EXECUTABLE STATEMENT DGESL
193  nm1 = n - 1
194  IF (job .NE. 0) GO TO 50
195 C
196 C JOB = 0 , SOLVE A * X = B
197 C FIRST SOLVE L*Y = B
198 C
199  IF (nm1 .LT. 1) GO TO 30
200  DO 20 k = 1, nm1
201  l = ipvt(k)
202  t = b(l)
203  IF (l .EQ. k) GO TO 10
204  b(l) = b(k)
205  b(k) = t
206  10 CONTINUE
207  CALL daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
208  20 CONTINUE
209  30 CONTINUE
210 C
211 C NOW SOLVE U*X = Y
212 C
213  DO 40 kb = 1, n
214  k = n + 1 - kb
215  b(k) = b(k)/a(k,k)
216  t = -b(k)
217  CALL daxpy(k-1,t,a(1,k),1,b(1),1)
218  40 CONTINUE
219  GO TO 100
220  50 CONTINUE
221 C
222 C JOB = NONZERO, SOLVE TRANS(A) * X = B
223 C FIRST SOLVE TRANS(U)*Y = B
224 C
225  DO 60 k = 1, n
226  t = ddot(k-1,a(1,k),1,b(1),1)
227  b(k) = (b(k) - t)/a(k,k)
228  60 CONTINUE
229 C
230 C NOW SOLVE TRANS(L)*X = Y
231 C
232  IF (nm1 .LT. 1) GO TO 90
233  DO 80 kb = 1, nm1
234  k = n - kb
235  b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
236  l = ipvt(k)
237  IF (l .EQ. k) GO TO 70
238  t = b(l)
239  b(l) = b(k)
240  b(k) = t
241  70 CONTINUE
242  80 CONTINUE
243  90 CONTINUE
244  100 CONTINUE
245  RETURN
246  END
247 *DECK DGBFA
248  SUBROUTINE dgbfa (ABD, LDA, N, ML, MU, IPVT, INFO)
249 C***BEGIN PROLOGUE DGBFA
250 C***PURPOSE Factor a band matrix using Gaussian elimination.
251 C***CATEGORY D2A2
252 C***TYPE DOUBLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
253 C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
254 C***AUTHOR Moler, C. B., (U. of New Mexico)
255 C***DESCRIPTION
256 C
257 C DGBFA factors a double precision band matrix by elimination.
258 C
259 C DGBFA is usually called by DGBCO, but it can be called
260 C directly with a saving in time if RCOND is not needed.
261 C
262 C On Entry
263 C
264 C ABD DOUBLE PRECISION(LDA, N)
265 C contains the matrix in band storage. The columns
266 C of the matrix are stored in the columns of ABD and
267 C the diagonals of the matrix are stored in rows
268 C ML+1 through 2*ML+MU+1 of ABD .
269 C See the comments below for details.
270 C
271 C LDA INTEGER
272 C the leading dimension of the array ABD .
273 C LDA must be .GE. 2*ML + MU + 1 .
274 C
275 C N INTEGER
276 C the order of the original matrix.
277 C
278 C ML INTEGER
279 C number of diagonals below the main diagonal.
280 C 0 .LE. ML .LT. N .
281 C
282 C MU INTEGER
283 C number of diagonals above the main diagonal.
284 C 0 .LE. MU .LT. N .
285 C More efficient if ML .LE. MU .
286 C On Return
287 C
288 C ABD an upper triangular matrix in band storage and
289 C the multipliers which were used to obtain it.
290 C The factorization can be written A = L*U where
291 C L is a product of permutation and unit lower
292 C triangular matrices and U is upper triangular.
293 C
294 C IPVT INTEGER(N)
295 C an integer vector of pivot indices.
296 C
297 C INFO INTEGER
298 C = 0 normal value.
299 C = K if U(K,K) .EQ. 0.0 . This is not an error
300 C condition for this subroutine, but it does
301 C indicate that DGBSL will divide by zero if
302 C called. Use RCOND in DGBCO for a reliable
303 C indication of singularity.
304 C
305 C Band Storage
306 C
307 C If A is a band matrix, the following program segment
308 C will set up the input.
309 C
310 C ML = (band width below the diagonal)
311 C MU = (band width above the diagonal)
312 C M = ML + MU + 1
313 C DO 20 J = 1, N
314 C I1 = MAX(1, J-MU)
315 C I2 = MIN(N, J+ML)
316 C DO 10 I = I1, I2
317 C K = I - J + M
318 C ABD(K,J) = A(I,J)
319 C 10 CONTINUE
320 C 20 CONTINUE
321 C
322 C This uses rows ML+1 through 2*ML+MU+1 of ABD .
323 C In addition, the first ML rows in ABD are used for
324 C elements generated during the triangularization.
325 C The total number of rows needed in ABD is 2*ML+MU+1 .
326 C The ML+MU by ML+MU upper left triangle and the
327 C ML by ML lower right triangle are not referenced.
328 C
329 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
330 C Stewart, LINPACK Users' Guide, SIAM, 1979.
331 C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX
332 C***REVISION HISTORY (YYMMDD)
333 C 780814 DATE WRITTEN
334 C 890531 Changed all specific intrinsics to generic. (WRB)
335 C 890831 Modified array declarations. (WRB)
336 C 890831 REVISION DATE from Version 3.2
337 C 891214 Prologue converted to Version 4.0 format. (BAB)
338 C 900326 Removed duplicate information from DESCRIPTION section.
339 C (WRB)
340 C 920501 Reformatted the REFERENCES section. (WRB)
341 C***END PROLOGUE DGBFA
342  INTEGER lda,n,ml,mu,ipvt(*),info
343  DOUBLE PRECISION abd(lda,*)
344 C
345  DOUBLE PRECISION t
346  INTEGER i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
347 C
348 C***FIRST EXECUTABLE STATEMENT DGBFA
349  m = ml + mu + 1
350  info = 0
351 C
352 C ZERO INITIAL FILL-IN COLUMNS
353 C
354  j0 = mu + 2
355  j1 = min(n,m) - 1
356  IF (j1 .LT. j0) GO TO 30
357  DO 20 jz = j0, j1
358  i0 = m + 1 - jz
359  DO 10 i = i0, ml
360  abd(i,jz) = 0.0d0
361  10 CONTINUE
362  20 CONTINUE
363  30 CONTINUE
364  jz = j1
365  ju = 0
366 C
367 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
368 C
369  nm1 = n - 1
370  IF (nm1 .LT. 1) GO TO 130
371  DO 120 k = 1, nm1
372  kp1 = k + 1
373 C
374 C ZERO NEXT FILL-IN COLUMN
375 C
376  jz = jz + 1
377  IF (jz .GT. n) GO TO 50
378  IF (ml .LT. 1) GO TO 50
379  DO 40 i = 1, ml
380  abd(i,jz) = 0.0d0
381  40 CONTINUE
382  50 CONTINUE
383 C
384 C FIND L = PIVOT INDEX
385 C
386  lm = min(ml,n-k)
387  l = idamax(lm+1,abd(m,k),1) + m - 1
388  ipvt(k) = l + k - m
389 C
390 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
391 C
392  IF (abd(l,k) .EQ. 0.0d0) GO TO 100
393 C
394 C INTERCHANGE IF NECESSARY
395 C
396  IF (l .EQ. m) GO TO 60
397  t = abd(l,k)
398  abd(l,k) = abd(m,k)
399  abd(m,k) = t
400  60 CONTINUE
401 C
402 C COMPUTE MULTIPLIERS
403 C
404  t = -1.0d0/abd(m,k)
405  CALL dscal(lm,t,abd(m+1,k),1)
406 C
407 C ROW ELIMINATION WITH COLUMN INDEXING
408 C
409  ju = min(max(ju,mu+ipvt(k)),n)
410  mm = m
411  IF (ju .LT. kp1) GO TO 90
412  DO 80 j = kp1, ju
413  l = l - 1
414  mm = mm - 1
415  t = abd(l,j)
416  IF (l .EQ. mm) GO TO 70
417  abd(l,j) = abd(mm,j)
418  abd(mm,j) = t
419  70 CONTINUE
420  CALL daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
421  80 CONTINUE
422  90 CONTINUE
423  GO TO 110
424  100 CONTINUE
425  info = k
426  110 CONTINUE
427  120 CONTINUE
428  130 CONTINUE
429  ipvt(n) = n
430  IF (abd(m,n) .EQ. 0.0d0) info = n
431  RETURN
432  END
433 *DECK DGBSL
434  SUBROUTINE dgbsl (ABD, LDA, N, ML, MU, IPVT, B, JOB)
435 C***BEGIN PROLOGUE DGBSL
436 C***PURPOSE Solve the real band system A*X=B or TRANS(A)*X=B using
437 C the factors computed by DGBCO or DGBFA.
438 C***CATEGORY D2A2
439 C***TYPE DOUBLE PRECISION (SGBSL-S, DGBSL-D, CGBSL-C)
440 C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
441 C***AUTHOR Moler, C. B., (U. of New Mexico)
442 C***DESCRIPTION
443 C
444 C DGBSL solves the double precision band system
445 C A * X = B or TRANS(A) * X = B
446 C using the factors computed by DGBCO or DGBFA.
447 C
448 C On Entry
449 C
450 C ABD DOUBLE PRECISION(LDA, N)
451 C the output from DGBCO or DGBFA.
452 C
453 C LDA INTEGER
454 C the leading dimension of the array ABD .
455 C
456 C N INTEGER
457 C the order of the original matrix.
458 C
459 C ML INTEGER
460 C number of diagonals below the main diagonal.
461 C
462 C MU INTEGER
463 C number of diagonals above the main diagonal.
464 C
465 C IPVT INTEGER(N)
466 C the pivot vector from DGBCO or DGBFA.
467 C
468 C B DOUBLE PRECISION(N)
469 C the right hand side vector.
470 C
471 C JOB INTEGER
472 C = 0 to solve A*X = B ,
473 C = nonzero to solve TRANS(A)*X = B , where
474 C TRANS(A) is the transpose.
475 C
476 C On Return
477 C
478 C B the solution vector X .
479 C
480 C Error Condition
481 C
482 C A division by zero will occur if the input factor contains a
483 C zero on the diagonal. Technically this indicates singularity
484 C but it is often caused by improper arguments or improper
485 C setting of LDA . It will not occur if the subroutines are
486 C called correctly and if DGBCO has set RCOND .GT. 0.0
487 C or DGBFA has set INFO .EQ. 0 .
488 C
489 C To compute INVERSE(A) * C where C is a matrix
490 C with P columns
491 C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
492 C IF (RCOND is too small) GO TO ...
493 C DO 10 J = 1, P
494 C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
495 C 10 CONTINUE
496 C
497 C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
498 C Stewart, LINPACK Users' Guide, SIAM, 1979.
499 C***ROUTINES CALLED DAXPY, DDOT
500 C***REVISION HISTORY (YYMMDD)
501 C 780814 DATE WRITTEN
502 C 890531 Changed all specific intrinsics to generic. (WRB)
503 C 890831 Modified array declarations. (WRB)
504 C 890831 REVISION DATE from Version 3.2
505 C 891214 Prologue converted to Version 4.0 format. (BAB)
506 C 900326 Removed duplicate information from DESCRIPTION section.
507 C (WRB)
508 C 920501 Reformatted the REFERENCES section. (WRB)
509 C***END PROLOGUE DGBSL
510  INTEGER lda,n,ml,mu,ipvt(*),job
511  DOUBLE PRECISION abd(lda,*),b(*)
512 C
513  DOUBLE PRECISION ddot,t
514  INTEGER k,kb,l,la,lb,lm,m,nm1
515 C***FIRST EXECUTABLE STATEMENT DGBSL
516  m = mu + ml + 1
517  nm1 = n - 1
518  IF (job .NE. 0) GO TO 50
519 C
520 C JOB = 0 , SOLVE A * X = B
521 C FIRST SOLVE L*Y = B
522 C
523  IF (ml .EQ. 0) GO TO 30
524  IF (nm1 .LT. 1) GO TO 30
525  DO 20 k = 1, nm1
526  lm = min(ml,n-k)
527  l = ipvt(k)
528  t = b(l)
529  IF (l .EQ. k) GO TO 10
530  b(l) = b(k)
531  b(k) = t
532  10 CONTINUE
533  CALL daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
534  20 CONTINUE
535  30 CONTINUE
536 C
537 C NOW SOLVE U*X = Y
538 C
539  DO 40 kb = 1, n
540  k = n + 1 - kb
541  b(k) = b(k)/abd(m,k)
542  lm = min(k,m) - 1
543  la = m - lm
544  lb = k - lm
545  t = -b(k)
546  CALL daxpy(lm,t,abd(la,k),1,b(lb),1)
547  40 CONTINUE
548  GO TO 100
549  50 CONTINUE
550 C
551 C JOB = NONZERO, SOLVE TRANS(A) * X = B
552 C FIRST SOLVE TRANS(U)*Y = B
553 C
554  DO 60 k = 1, n
555  lm = min(k,m) - 1
556  la = m - lm
557  lb = k - lm
558  t = ddot(lm,abd(la,k),1,b(lb),1)
559  b(k) = (b(k) - t)/abd(m,k)
560  60 CONTINUE
561 C
562 C NOW SOLVE TRANS(L)*X = Y
563 C
564  IF (ml .EQ. 0) GO TO 90
565  IF (nm1 .LT. 1) GO TO 90
566  DO 80 kb = 1, nm1
567  k = n - kb
568  lm = min(ml,n-k)
569  b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
570  l = ipvt(k)
571  IF (l .EQ. k) GO TO 70
572  t = b(l)
573  b(l) = b(k)
574  b(k) = t
575  70 CONTINUE
576  80 CONTINUE
577  90 CONTINUE
578  100 CONTINUE
579  RETURN
580  END
581 *DECK DAXPY
582  SUBROUTINE daxpy (N, DA, DX, INCX, DY, INCY)
583 C***BEGIN PROLOGUE DAXPY
584 C***PURPOSE Compute a constant times a vector plus a vector.
585 C***CATEGORY D1A7
586 C***TYPE DOUBLE PRECISION (SAXPY-S, DAXPY-D, CAXPY-C)
587 C***KEYWORDS BLAS, LINEAR ALGEBRA, TRIAD, VECTOR
588 C***AUTHOR Lawson, C. L., (JPL)
589 C Hanson, R. J., (SNLA)
590 C Kincaid, D. R., (U. of Texas)
591 C Krogh, F. T., (JPL)
592 C***DESCRIPTION
593 C
594 C B L A S Subprogram
595 C Description of Parameters
596 C
597 C --Input--
598 C N number of elements in input vector(s)
599 C DA double precision scalar multiplier
600 C DX double precision vector with N elements
601 C INCX storage spacing between elements of DX
602 C DY double precision vector with N elements
603 C INCY storage spacing between elements of DY
604 C
605 C --Output--
606 C DY double precision result (unchanged if N .LE. 0)
607 C
608 C Overwrite double precision DY with double precision DA*DX + DY.
609 C For I = 0 to N-1, replace DY(LY+I*INCY) with DA*DX(LX+I*INCX) +
610 C DY(LY+I*INCY),
611 C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
612 C defined in a similar way using INCY.
613 C
614 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
615 C Krogh, Basic linear algebra subprograms for Fortran
616 C usage, Algorithm No. 539, Transactions on Mathematical
617 C Software 5, 3 (September 1979), pp. 308-323.
618 C***ROUTINES CALLED (NONE)
619 C***REVISION HISTORY (YYMMDD)
620 C 791001 DATE WRITTEN
621 C 890831 Modified array declarations. (WRB)
622 C 890831 REVISION DATE from Version 3.2
623 C 891214 Prologue converted to Version 4.0 format. (BAB)
624 C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
625 C 920501 Reformatted the REFERENCES section. (WRB)
626 C***END PROLOGUE DAXPY
627  DOUBLE PRECISION dx(*), dy(*), da
628 C***FIRST EXECUTABLE STATEMENT DAXPY
629  IF (n.LE.0 .OR. da.EQ.0.0d0) RETURN
630  IF (incx .EQ. incy) IF (incx-1) 5,20,60
631 C
632 C Code for unequal or nonpositive increments.
633 C
634  5 ix = 1
635  iy = 1
636  IF (incx .LT. 0) ix = (-n+1)*incx + 1
637  IF (incy .LT. 0) iy = (-n+1)*incy + 1
638  DO 10 i = 1,n
639  dy(iy) = dy(iy) + da*dx(ix)
640  ix = ix + incx
641  iy = iy + incy
642  10 CONTINUE
643  RETURN
644 C
645 C Code for both increments equal to 1.
646 C
647 C Clean-up loop so remaining vector length is a multiple of 4.
648 C
649  20 m = mod(n,4)
650  IF (m .EQ. 0) GO TO 40
651  DO 30 i = 1,m
652  dy(i) = dy(i) + da*dx(i)
653  30 CONTINUE
654  IF (n .LT. 4) RETURN
655  40 mp1 = m + 1
656  DO 50 i = mp1,n,4
657  dy(i) = dy(i) + da*dx(i)
658  dy(i+1) = dy(i+1) + da*dx(i+1)
659  dy(i+2) = dy(i+2) + da*dx(i+2)
660  dy(i+3) = dy(i+3) + da*dx(i+3)
661  50 CONTINUE
662  RETURN
663 C
664 C Code for equal, positive, non-unit increments.
665 C
666  60 ns = n*incx
667  DO 70 i = 1,ns,incx
668  dy(i) = da*dx(i) + dy(i)
669  70 CONTINUE
670  RETURN
671  END
672 *DECK DCOPY
673  SUBROUTINE dcopy (N, DX, INCX, DY, INCY)
674 C***BEGIN PROLOGUE DCOPY
675 C***PURPOSE Copy a vector.
676 C***CATEGORY D1A5
677 C***TYPE DOUBLE PRECISION (SCOPY-S, DCOPY-D, CCOPY-C, ICOPY-I)
678 C***KEYWORDS BLAS, COPY, LINEAR ALGEBRA, VECTOR
679 C***AUTHOR Lawson, C. L., (JPL)
680 C Hanson, R. J., (SNLA)
681 C Kincaid, D. R., (U. of Texas)
682 C Krogh, F. T., (JPL)
683 C***DESCRIPTION
684 C
685 C B L A S Subprogram
686 C Description of Parameters
687 C
688 C --Input--
689 C N number of elements in input vector(s)
690 C DX double precision vector with N elements
691 C INCX storage spacing between elements of DX
692 C DY double precision vector with N elements
693 C INCY storage spacing between elements of DY
694 C
695 C --Output--
696 C DY copy of vector DX (unchanged if N .LE. 0)
697 C
698 C Copy double precision DX to double precision DY.
699 C For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY),
700 C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
701 C defined in a similar way using INCY.
702 C
703 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
704 C Krogh, Basic linear algebra subprograms for Fortran
705 C usage, Algorithm No. 539, Transactions on Mathematical
706 C Software 5, 3 (September 1979), pp. 308-323.
707 C***ROUTINES CALLED (NONE)
708 C***REVISION HISTORY (YYMMDD)
709 C 791001 DATE WRITTEN
710 C 890831 Modified array declarations. (WRB)
711 C 890831 REVISION DATE from Version 3.2
712 C 891214 Prologue converted to Version 4.0 format. (BAB)
713 C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
714 C 920501 Reformatted the REFERENCES section. (WRB)
715 C***END PROLOGUE DCOPY
716  DOUBLE PRECISION dx(*), dy(*)
717 C***FIRST EXECUTABLE STATEMENT DCOPY
718  IF (n .LE. 0) RETURN
719  IF (incx .EQ. incy) IF (incx-1) 5,20,60
720 C
721 C Code for unequal or nonpositive increments.
722 C
723  5 ix = 1
724  iy = 1
725  IF (incx .LT. 0) ix = (-n+1)*incx + 1
726  IF (incy .LT. 0) iy = (-n+1)*incy + 1
727  DO 10 i = 1,n
728  dy(iy) = dx(ix)
729  ix = ix + incx
730  iy = iy + incy
731  10 CONTINUE
732  RETURN
733 C
734 C Code for both increments equal to 1.
735 C
736 C Clean-up loop so remaining vector length is a multiple of 7.
737 C
738  20 m = mod(n,7)
739  IF (m .EQ. 0) GO TO 40
740  DO 30 i = 1,m
741  dy(i) = dx(i)
742  30 CONTINUE
743  IF (n .LT. 7) RETURN
744  40 mp1 = m + 1
745  DO 50 i = mp1,n,7
746  dy(i) = dx(i)
747  dy(i+1) = dx(i+1)
748  dy(i+2) = dx(i+2)
749  dy(i+3) = dx(i+3)
750  dy(i+4) = dx(i+4)
751  dy(i+5) = dx(i+5)
752  dy(i+6) = dx(i+6)
753  50 CONTINUE
754  RETURN
755 C
756 C Code for equal, positive, non-unit increments.
757 C
758  60 ns = n*incx
759  DO 70 i = 1,ns,incx
760  dy(i) = dx(i)
761  70 CONTINUE
762  RETURN
763  END
764 *DECK DDOT
765  DOUBLE PRECISION FUNCTION ddot (N, DX, INCX, DY, INCY)
766 C***BEGIN PROLOGUE DDOT
767 C***PURPOSE Compute the inner product of two vectors.
768 C***CATEGORY D1A4
769 C***TYPE DOUBLE PRECISION (SDOT-S, DDOT-D, CDOTU-C)
770 C***KEYWORDS BLAS, INNER PRODUCT, LINEAR ALGEBRA, VECTOR
771 C***AUTHOR Lawson, C. L., (JPL)
772 C Hanson, R. J., (SNLA)
773 C Kincaid, D. R., (U. of Texas)
774 C Krogh, F. T., (JPL)
775 C***DESCRIPTION
776 C
777 C B L A S Subprogram
778 C Description of Parameters
779 C
780 C --Input--
781 C N number of elements in input vector(s)
782 C DX double precision vector with N elements
783 C INCX storage spacing between elements of DX
784 C DY double precision vector with N elements
785 C INCY storage spacing between elements of DY
786 C
787 C --Output--
788 C DDOT double precision dot product (zero if N .LE. 0)
789 C
790 C Returns the dot product of double precision DX and DY.
791 C DDOT = sum for I = 0 to N-1 of DX(LX+I*INCX) * DY(LY+I*INCY),
792 C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
793 C defined in a similar way using INCY.
794 C
795 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
796 C Krogh, Basic linear algebra subprograms for Fortran
797 C usage, Algorithm No. 539, Transactions on Mathematical
798 C Software 5, 3 (September 1979), pp. 308-323.
799 C***ROUTINES CALLED (NONE)
800 C***REVISION HISTORY (YYMMDD)
801 C 791001 DATE WRITTEN
802 C 890831 Modified array declarations. (WRB)
803 C 890831 REVISION DATE from Version 3.2
804 C 891214 Prologue converted to Version 4.0 format. (BAB)
805 C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
806 C 920501 Reformatted the REFERENCES section. (WRB)
807 C***END PROLOGUE DDOT
808  DOUBLE PRECISION dx(*), dy(*)
809 C***FIRST EXECUTABLE STATEMENT DDOT
810  ddot = 0.0d0
811  IF (n .LE. 0) RETURN
812  IF (incx .EQ. incy) IF (incx-1) 5,20,60
813 C
814 C Code for unequal or nonpositive increments.
815 C
816  5 ix = 1
817  iy = 1
818  IF (incx .LT. 0) ix = (-n+1)*incx + 1
819  IF (incy .LT. 0) iy = (-n+1)*incy + 1
820  DO 10 i = 1,n
821  ddot = ddot + dx(ix)*dy(iy)
822  ix = ix + incx
823  iy = iy + incy
824  10 CONTINUE
825  RETURN
826 C
827 C Code for both increments equal to 1.
828 C
829 C Clean-up loop so remaining vector length is a multiple of 5.
830 C
831  20 m = mod(n,5)
832  IF (m .EQ. 0) GO TO 40
833  DO 30 i = 1,m
834  ddot = ddot + dx(i)*dy(i)
835  30 CONTINUE
836  IF (n .LT. 5) RETURN
837  40 mp1 = m + 1
838  DO 50 i = mp1,n,5
839  ddot = ddot + dx(i)*dy(i) + dx(i+1)*dy(i+1) + dx(i+2)*dy(i+2) +
840  1 dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4)
841  50 CONTINUE
842  RETURN
843 C
844 C Code for equal, positive, non-unit increments.
845 C
846  60 ns = n*incx
847  DO 70 i = 1,ns,incx
848  ddot = ddot + dx(i)*dy(i)
849  70 CONTINUE
850  RETURN
851  END
852 *DECK DNRM2
853  DOUBLE PRECISION FUNCTION dnrm2 (N, DX, INCX)
854 C***BEGIN PROLOGUE DNRM2
855 C***PURPOSE Compute the Euclidean length (L2 norm) of a vector.
856 C***CATEGORY D1A3B
857 C***TYPE DOUBLE PRECISION (SNRM2-S, DNRM2-D, SCNRM2-C)
858 C***KEYWORDS BLAS, EUCLIDEAN LENGTH, EUCLIDEAN NORM, L2,
859 C LINEAR ALGEBRA, UNITARY, VECTOR
860 C***AUTHOR Lawson, C. L., (JPL)
861 C Hanson, R. J., (SNLA)
862 C Kincaid, D. R., (U. of Texas)
863 C Krogh, F. T., (JPL)
864 C***DESCRIPTION
865 C
866 C B L A S Subprogram
867 C Description of parameters
868 C
869 C --Input--
870 C N number of elements in input vector(s)
871 C DX double precision vector with N elements
872 C INCX storage spacing between elements of DX
873 C
874 C --Output--
875 C DNRM2 double precision result (zero if N .LE. 0)
876 C
877 C Euclidean norm of the N-vector stored in DX with storage
878 C increment INCX.
879 C If N .LE. 0, return with result = 0.
880 C If N .GE. 1, then INCX must be .GE. 1
881 C
882 C Four phase method using two built-in constants that are
883 C hopefully applicable to all machines.
884 C CUTLO = maximum of SQRT(U/EPS) over all known machines.
885 C CUTHI = minimum of SQRT(V) over all known machines.
886 C where
887 C EPS = smallest no. such that EPS + 1. .GT. 1.
888 C U = smallest positive no. (underflow limit)
889 C V = largest no. (overflow limit)
890 C
891 C Brief outline of algorithm.
892 C
893 C Phase 1 scans zero components.
894 C move to phase 2 when a component is nonzero and .LE. CUTLO
895 C move to phase 3 when a component is .GT. CUTLO
896 C move to phase 4 when a component is .GE. CUTHI/M
897 C where M = N for X() real and M = 2*N for complex.
898 C
899 C Values for CUTLO and CUTHI.
900 C From the environmental parameters listed in the IMSL converter
901 C document the limiting values are as follows:
902 C CUTLO, S.P. U/EPS = 2**(-102) for Honeywell. Close seconds are
903 C Univac and DEC at 2**(-103)
904 C Thus CUTLO = 2**(-51) = 4.44089E-16
905 C CUTHI, S.P. V = 2**127 for Univac, Honeywell, and DEC.
906 C Thus CUTHI = 2**(63.5) = 1.30438E19
907 C CUTLO, D.P. U/EPS = 2**(-67) for Honeywell and DEC.
908 C Thus CUTLO = 2**(-33.5) = 8.23181D-11
909 C CUTHI, D.P. same as S.P. CUTHI = 1.30438D19
910 C DATA CUTLO, CUTHI /8.232D-11, 1.304D19/
911 C DATA CUTLO, CUTHI /4.441E-16, 1.304E19/
912 C
913 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
914 C Krogh, Basic linear algebra subprograms for Fortran
915 C usage, Algorithm No. 539, Transactions on Mathematical
916 C Software 5, 3 (September 1979), pp. 308-323.
917 C***ROUTINES CALLED (NONE)
918 C***REVISION HISTORY (YYMMDD)
919 C 791001 DATE WRITTEN
920 C 890531 Changed all specific intrinsics to generic. (WRB)
921 C 890831 Modified array declarations. (WRB)
922 C 890831 REVISION DATE from Version 3.2
923 C 891214 Prologue converted to Version 4.0 format. (BAB)
924 C 920501 Reformatted the REFERENCES section. (WRB)
925 C***END PROLOGUE DNRM2
926  INTEGER next
927  DOUBLE PRECISION dx(*), cutlo, cuthi, hitest, sum, xmax, zero,
928  + one
929  SAVE cutlo, cuthi, zero, one
930  DATA zero, one /0.0d0, 1.0d0/
931 C
932  DATA cutlo, cuthi /8.232d-11, 1.304d19/
933 C***FIRST EXECUTABLE STATEMENT DNRM2
934  IF (n .GT. 0) GO TO 10
935  dnrm2 = zero
936  GO TO 300
937 C
938  10 assign 30 to next
939  sum = zero
940  nn = n * incx
941 C
942 C BEGIN MAIN LOOP
943 C
944  i = 1
945  20 GO TO next,(30, 50, 70, 110)
946  30 IF (abs(dx(i)) .GT. cutlo) GO TO 85
947  assign 50 to next
948  xmax = zero
949 C
950 C PHASE 1. SUM IS ZERO
951 C
952  50 IF (dx(i) .EQ. zero) GO TO 200
953  IF (abs(dx(i)) .GT. cutlo) GO TO 85
954 C
955 C PREPARE FOR PHASE 2.
956 C
957  assign 70 to next
958  GO TO 105
959 C
960 C PREPARE FOR PHASE 4.
961 C
962  100 i = j
963  assign 110 to next
964  sum = (sum / dx(i)) / dx(i)
965  105 xmax = abs(dx(i))
966  GO TO 115
967 C
968 C PHASE 2. SUM IS SMALL.
969 C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
970 C
971  70 IF (abs(dx(i)) .GT. cutlo) GO TO 75
972 C
973 C COMMON CODE FOR PHASES 2 AND 4.
974 C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
975 C
976  110 IF (abs(dx(i)) .LE. xmax) GO TO 115
977  sum = one + sum * (xmax / dx(i))**2
978  xmax = abs(dx(i))
979  GO TO 200
980 C
981  115 sum = sum + (dx(i)/xmax)**2
982  GO TO 200
983 C
984 C PREPARE FOR PHASE 3.
985 C
986  75 sum = (sum * xmax) * xmax
987 C
988 C FOR REAL OR D.P. SET HITEST = CUTHI/N
989 C FOR COMPLEX SET HITEST = CUTHI/(2*N)
990 C
991  85 hitest = cuthi / n
992 C
993 C PHASE 3. SUM IS MID-RANGE. NO SCALING.
994 C
995  DO 95 j = i,nn,incx
996  IF (abs(dx(j)) .GE. hitest) GO TO 100
997  95 sum = sum + dx(j)**2
998  dnrm2 = sqrt(sum)
999  GO TO 300
1000 C
1001  200 CONTINUE
1002  i = i + incx
1003  IF (i .LE. nn) GO TO 20
1004 C
1005 C END OF MAIN LOOP.
1006 C
1007 C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
1008 C
1009  dnrm2 = xmax * sqrt(sum)
1010  300 CONTINUE
1011  RETURN
1012  END
1013 *DECK DSCAL
1014  SUBROUTINE dscal (N, DA, DX, INCX)
1015 C***BEGIN PROLOGUE DSCAL
1016 C***PURPOSE Multiply a vector by a constant.
1017 C***CATEGORY D1A6
1018 C***TYPE DOUBLE PRECISION (SSCAL-S, DSCAL-D, CSCAL-C)
1019 C***KEYWORDS BLAS, LINEAR ALGEBRA, SCALE, VECTOR
1020 C***AUTHOR Lawson, C. L., (JPL)
1021 C Hanson, R. J., (SNLA)
1022 C Kincaid, D. R., (U. of Texas)
1023 C Krogh, F. T., (JPL)
1024 C***DESCRIPTION
1025 C
1026 C B L A S Subprogram
1027 C Description of Parameters
1028 C
1029 C --Input--
1030 C N number of elements in input vector(s)
1031 C DA double precision scale factor
1032 C DX double precision vector with N elements
1033 C INCX storage spacing between elements of DX
1034 C
1035 C --Output--
1036 C DX double precision result (unchanged if N.LE.0)
1037 C
1038 C Replace double precision DX by double precision DA*DX.
1039 C For I = 0 to N-1, replace DX(IX+I*INCX) with DA * DX(IX+I*INCX),
1040 C where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX.
1041 C
1042 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
1043 C Krogh, Basic linear algebra subprograms for Fortran
1044 C usage, Algorithm No. 539, Transactions on Mathematical
1045 C Software 5, 3 (September 1979), pp. 308-323.
1046 C***ROUTINES CALLED (NONE)
1047 C***REVISION HISTORY (YYMMDD)
1048 C 791001 DATE WRITTEN
1049 C 890831 Modified array declarations. (WRB)
1050 C 890831 REVISION DATE from Version 3.2
1051 C 891214 Prologue converted to Version 4.0 format. (BAB)
1052 C 900821 Modified to correct problem with a negative increment.
1053 C (WRB)
1054 C 920501 Reformatted the REFERENCES section. (WRB)
1055 C***END PROLOGUE DSCAL
1056  DOUBLE PRECISION da, dx(*)
1057  INTEGER i, incx, ix, m, mp1, n
1058 C***FIRST EXECUTABLE STATEMENT DSCAL
1059  IF (n .LE. 0) RETURN
1060  IF (incx .EQ. 1) GOTO 20
1061 C
1062 C Code for increment not equal to 1.
1063 C
1064  ix = 1
1065  IF (incx .LT. 0) ix = (-n+1)*incx + 1
1066  DO 10 i = 1,n
1067  dx(ix) = da*dx(ix)
1068  ix = ix + incx
1069  10 CONTINUE
1070  RETURN
1071 C
1072 C Code for increment equal to 1.
1073 C
1074 C Clean-up loop so remaining vector length is a multiple of 5.
1075 C
1076  20 m = mod(n,5)
1077  IF (m .EQ. 0) GOTO 40
1078  DO 30 i = 1,m
1079  dx(i) = da*dx(i)
1080  30 CONTINUE
1081  IF (n .LT. 5) RETURN
1082  40 mp1 = m + 1
1083  DO 50 i = mp1,n,5
1084  dx(i) = da*dx(i)
1085  dx(i+1) = da*dx(i+1)
1086  dx(i+2) = da*dx(i+2)
1087  dx(i+3) = da*dx(i+3)
1088  dx(i+4) = da*dx(i+4)
1089  50 CONTINUE
1090  RETURN
1091  END
1092 *DECK IDAMAX
1093  INTEGER FUNCTION idamax (N, DX, INCX)
1094 C***BEGIN PROLOGUE IDAMAX
1095 C***PURPOSE Find the smallest index of that component of a vector
1096 C having the maximum magnitude.
1097 C***CATEGORY D1A2
1098 C***TYPE DOUBLE PRECISION (ISAMAX-S, IDAMAX-D, ICAMAX-C)
1099 C***KEYWORDS BLAS, LINEAR ALGEBRA, MAXIMUM COMPONENT, VECTOR
1100 C***AUTHOR Lawson, C. L., (JPL)
1101 C Hanson, R. J., (SNLA)
1102 C Kincaid, D. R., (U. of Texas)
1103 C Krogh, F. T., (JPL)
1104 C***DESCRIPTION
1105 C
1106 C B L A S Subprogram
1107 C Description of Parameters
1108 C
1109 C --Input--
1110 C N number of elements in input vector(s)
1111 C DX double precision vector with N elements
1112 C INCX storage spacing between elements of DX
1113 C
1114 C --Output--
1115 C IDAMAX smallest index (zero if N .LE. 0)
1116 C
1117 C Find smallest index of maximum magnitude of double precision DX.
1118 C IDAMAX = first I, I = 1 to N, to maximize ABS(DX(IX+(I-1)*INCX)),
1119 C where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX.
1120 C
1121 C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
1122 C Krogh, Basic linear algebra subprograms for Fortran
1123 C usage, Algorithm No. 539, Transactions on Mathematical
1124 C Software 5, 3 (September 1979), pp. 308-323.
1125 C***ROUTINES CALLED (NONE)
1126 C***REVISION HISTORY (YYMMDD)
1127 C 791001 DATE WRITTEN
1128 C 890531 Changed all specific intrinsics to generic. (WRB)
1129 C 890531 REVISION DATE from Version 3.2
1130 C 891214 Prologue converted to Version 4.0 format. (BAB)
1131 C 900821 Modified to correct problem with a negative increment.
1132 C (WRB)
1133 C 920501 Reformatted the REFERENCES section. (WRB)
1134 C***END PROLOGUE IDAMAX
1135  DOUBLE PRECISION dx(*), dmax, xmag
1136  INTEGER i, incx, ix, n
1137 C***FIRST EXECUTABLE STATEMENT IDAMAX
1138  idamax = 0
1139  IF (n .LE. 0) RETURN
1140  idamax = 1
1141  IF (n .EQ. 1) RETURN
1142 C
1143  IF (incx .EQ. 1) GOTO 20
1144 C
1145 C Code for increments not equal to 1.
1146 C
1147  ix = 1
1148  IF (incx .LT. 0) ix = (-n+1)*incx + 1
1149  dmax = abs(dx(ix))
1150  ix = ix + incx
1151  DO 10 i = 2,n
1152  xmag = abs(dx(ix))
1153  IF (xmag .GT. dmax) THEN
1154  idamax = i
1155  dmax = xmag
1156  ENDIF
1157  ix = ix + incx
1158  10 CONTINUE
1159  RETURN
1160 C
1161 C Code for increments equal to 1.
1162 C
1163  20 dmax = abs(dx(1))
1164  DO 30 i = 2,n
1165  xmag = abs(dx(i))
1166  IF (xmag .GT. dmax) THEN
1167  idamax = i
1168  dmax = xmag
1169  ENDIF
1170  30 CONTINUE
1171  RETURN
1172  END
1173 *DECK XERRWD
1174  SUBROUTINE xerrwd (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
1175 C***BEGIN PROLOGUE XERRWD
1176 C***SUBSIDIARY
1177 C***PURPOSE Write error message with values.
1178 C***CATEGORY R3C
1179 C***TYPE DOUBLE PRECISION (XERRWV-S, XERRWD-D)
1180 C***AUTHOR Hindmarsh, Alan C., (LLNL)
1181 C***DESCRIPTION
1182 C
1183 C Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV,
1184 C as given here, constitute a simplified version of the SLATEC error
1185 C handling package.
1186 C
1187 C All arguments are input arguments.
1188 C
1189 C MSG = The message (character array).
1190 C NMES = The length of MSG (number of characters).
1191 C NERR = The error number (not used).
1192 C LEVEL = The error level..
1193 C 0 or 1 means recoverable (control returns to caller).
1194 C 2 means fatal (run is aborted--see note below).
1195 C NI = Number of integers (0, 1, or 2) to be printed with message.
1196 C I1,I2 = Integers to be printed, depending on NI.
1197 C NR = Number of reals (0, 1, or 2) to be printed with message.
1198 C R1,R2 = Reals to be printed, depending on NR.
1199 C
1200 C Note.. this routine is machine-dependent and specialized for use
1201 C in limited context, in the following ways..
1202 C 1. The argument MSG is assumed to be of type CHARACTER, and
1203 C the message is printed with a format of (1X,A).
1204 C 2. The message is assumed to take only one line.
1205 C Multi-line messages are generated by repeated calls.
1206 C 3. If LEVEL = 2, control passes to the statement STOP
1207 C to abort the run. This statement may be machine-dependent.
1208 C 4. R1 and R2 are assumed to be in double precision and are printed
1209 C in D21.13 format.
1210 C
1211 C***ROUTINES CALLED IXSAV
1212 C***REVISION HISTORY (YYMMDD)
1213 C 920831 DATE WRITTEN
1214 C 921118 Replaced MFLGSV/LUNSAV by IXSAV. (ACH)
1215 C 930329 Modified prologue to SLATEC format. (FNF)
1216 C 930407 Changed MSG from CHARACTER*1 array to variable. (FNF)
1217 C 930922 Minor cosmetic change. (FNF)
1218 C***END PROLOGUE XERRWD
1219 C
1220 C*Internal Notes:
1221 C
1222 C For a different default logical unit number, IXSAV (or a subsidiary
1223 C routine that it calls) will need to be modified.
1224 C For a different run-abort command, change the statement following
1225 C statement 100 at the end.
1226 C-----------------------------------------------------------------------
1227 C Subroutines called by XERRWD.. None
1228 C Function routine called by XERRWD.. IXSAV
1229 C-----------------------------------------------------------------------
1230 C**End
1231 C
1232 C Declare arguments.
1233 C
1234  DOUBLE PRECISION r1, r2
1235  INTEGER nmes, nerr, level, ni, i1, i2, nr
1236  CHARACTER*(*) msg
1237 C
1238 C Declare local variables.
1239 C
1240  INTEGER lunit, ixsav, mesflg
1241 C
1242 C Get logical unit number and message print flag.
1243 C
1244 C***FIRST EXECUTABLE STATEMENT XERRWD
1245  lunit = ixsav(1, 0, .false.)
1246  mesflg = ixsav(2, 0, .false.)
1247  IF (mesflg .EQ. 0) GO TO 100
1248 C
1249 C Write the message.
1250 C
1251  WRITE (lunit,10) msg
1252  10 FORMAT(1x,a)
1253  IF (ni .EQ. 1) WRITE (lunit, 20) i1
1254  20 FORMAT(6x,'In above message, I1 =',i10)
1255  IF (ni .EQ. 2) WRITE (lunit, 30) i1,i2
1256  30 FORMAT(6x,'In above message, I1 =',i10,3x,'I2 =',i10)
1257  IF (nr .EQ. 1) WRITE (lunit, 40) r1
1258  40 FORMAT(6x,'In above message, R1 =',d21.13)
1259  IF (nr .EQ. 2) WRITE (lunit, 50) r1,r2
1260  50 FORMAT(6x,'In above, R1 =',d21.13,3x,'R2 =',d21.13)
1261 C
1262 C Abort the run if LEVEL = 2.
1263 C
1264  100 IF (level .NE. 2) RETURN
1265  stop
1266 C----------------------- End of Subroutine XERRWD ----------------------
1267  END
1268 *DECK XSETF
1269  SUBROUTINE xsetf (MFLAG)
1270 C***BEGIN PROLOGUE XSETF
1271 C***PURPOSE Reset the error print control flag.
1272 C***CATEGORY R3A
1273 C***TYPE ALL (XSETF-A)
1274 C***KEYWORDS ERROR CONTROL
1275 C***AUTHOR Hindmarsh, Alan C., (LLNL)
1276 C***DESCRIPTION
1277 C
1278 C XSETF sets the error print control flag to MFLAG:
1279 C MFLAG=1 means print all messages (the default).
1280 C MFLAG=0 means no printing.
1281 C
1282 C***SEE ALSO XERRWD, XERRWV
1283 C***REFERENCES (NONE)
1284 C***ROUTINES CALLED IXSAV
1285 C***REVISION HISTORY (YYMMDD)
1286 C 921118 DATE WRITTEN
1287 C 930329 Added SLATEC format prologue. (FNF)
1288 C 930407 Corrected SEE ALSO section. (FNF)
1289 C 930922 Made user-callable, and other cosmetic changes. (FNF)
1290 C***END PROLOGUE XSETF
1291 C
1292 C Subroutines called by XSETF.. None
1293 C Function routine called by XSETF.. IXSAV
1294 C-----------------------------------------------------------------------
1295 C**End
1296  INTEGER mflag, junk, ixsav
1297 C
1298 C***FIRST EXECUTABLE STATEMENT XSETF
1299  IF (mflag .EQ. 0 .OR. mflag .EQ. 1) junk = ixsav(2,mflag,.true.)
1300  RETURN
1301 C----------------------- End of Subroutine XSETF -----------------------
1302  END
1303 *DECK XSETUN
1304  SUBROUTINE xsetun (LUN)
1305 C***BEGIN PROLOGUE XSETUN
1306 C***PURPOSE Reset the logical unit number for error messages.
1307 C***CATEGORY R3B
1308 C***TYPE ALL (XSETUN-A)
1309 C***KEYWORDS ERROR CONTROL
1310 C***DESCRIPTION
1311 C
1312 C XSETUN sets the logical unit number for error messages to LUN.
1313 C
1314 C***AUTHOR Hindmarsh, Alan C., (LLNL)
1315 C***SEE ALSO XERRWD, XERRWV
1316 C***REFERENCES (NONE)
1317 C***ROUTINES CALLED IXSAV
1318 C***REVISION HISTORY (YYMMDD)
1319 C 921118 DATE WRITTEN
1320 C 930329 Added SLATEC format prologue. (FNF)
1321 C 930407 Corrected SEE ALSO section. (FNF)
1322 C 930922 Made user-callable, and other cosmetic changes. (FNF)
1323 C***END PROLOGUE XSETUN
1324 C
1325 C Subroutines called by XSETUN.. None
1326 C Function routine called by XSETUN.. IXSAV
1327 C-----------------------------------------------------------------------
1328 C**End
1329  INTEGER lun, junk, ixsav
1330 C
1331 C***FIRST EXECUTABLE STATEMENT XSETUN
1332  IF (lun .GT. 0) junk = ixsav(1,lun,.true.)
1333  RETURN
1334 C----------------------- End of Subroutine XSETUN ----------------------
1335  END
1336 *DECK IXSAV
1337  INTEGER FUNCTION ixsav (IPAR, IVALUE, ISET)
1338 C***BEGIN PROLOGUE IXSAV
1339 C***SUBSIDIARY
1340 C***PURPOSE Save and recall error message control parameters.
1341 C***CATEGORY R3C
1342 C***TYPE ALL (IXSAV-A)
1343 C***AUTHOR Hindmarsh, Alan C., (LLNL)
1344 C***DESCRIPTION
1345 C
1346 C IXSAV saves and recalls one of two error message parameters:
1347 C LUNIT, the logical unit number to which messages are printed, and
1348 C MESFLG, the message print flag.
1349 C This is a modification of the SLATEC library routine J4SAVE.
1350 C
1351 C Saved local variables..
1352 C LUNIT = Logical unit number for messages. The default is obtained
1353 C by a call to IUMACH (may be machine-dependent).
1354 C MESFLG = Print control flag..
1355 C 1 means print all messages (the default).
1356 C 0 means no printing.
1357 C
1358 C On input..
1359 C IPAR = Parameter indicator (1 for LUNIT, 2 for MESFLG).
1360 C IVALUE = The value to be set for the parameter, if ISET = .TRUE.
1361 C ISET = Logical flag to indicate whether to read or write.
1362 C If ISET = .TRUE., the parameter will be given
1363 C the value IVALUE. If ISET = .FALSE., the parameter
1364 C will be unchanged, and IVALUE is a dummy argument.
1365 C
1366 C On return..
1367 C IXSAV = The (old) value of the parameter.
1368 C
1369 C***SEE ALSO XERRWD, XERRWV
1370 C***ROUTINES CALLED IUMACH
1371 C***REVISION HISTORY (YYMMDD)
1372 C 921118 DATE WRITTEN
1373 C 930329 Modified prologue to SLATEC format. (FNF)
1374 C 930915 Added IUMACH call to get default output unit. (ACH)
1375 C 930922 Minor cosmetic changes. (FNF)
1376 C 010425 Type declaration for IUMACH added. (ACH)
1377 C***END PROLOGUE IXSAV
1378 C
1379 C Subroutines called by IXSAV.. None
1380 C Function routine called by IXSAV.. IUMACH
1381 C-----------------------------------------------------------------------
1382 C**End
1383  LOGICAL iset
1384  INTEGER ipar, ivalue
1385 C-----------------------------------------------------------------------
1386  INTEGER iumach, lunit, mesflg
1387 C-----------------------------------------------------------------------
1388 C The following Fortran-77 declaration is to cause the values of the
1389 C listed (local) variables to be saved between calls to this routine.
1390 C-----------------------------------------------------------------------
1391  SAVE lunit, mesflg
1392  DATA lunit/-1/, mesflg/1/
1393 C
1394 C***FIRST EXECUTABLE STATEMENT IXSAV
1395  IF (ipar .EQ. 1) THEN
1396  IF (lunit .EQ. -1) lunit = iumach()
1397  ixsav = lunit
1398  IF (iset) lunit = ivalue
1399  ENDIF
1400 C
1401  IF (ipar .EQ. 2) THEN
1402  ixsav = mesflg
1403  IF (iset) mesflg = ivalue
1404  ENDIF
1405 C
1406  RETURN
1407 C----------------------- End of Function IXSAV -------------------------
1408  END
1409 *DECK IUMACH
1410  INTEGER FUNCTION iumach()
1411 C***BEGIN PROLOGUE IUMACH
1412 C***PURPOSE Provide standard output unit number.
1413 C***CATEGORY R1
1414 C***TYPE INTEGER (IUMACH-I)
1415 C***KEYWORDS MACHINE CONSTANTS
1416 C***AUTHOR Hindmarsh, Alan C., (LLNL)
1417 C***DESCRIPTION
1418 C *Usage:
1419 C INTEGER LOUT, IUMACH
1420 C LOUT = IUMACH()
1421 C
1422 C *Function Return Values:
1423 C LOUT : the standard logical unit for Fortran output.
1424 C
1425 C***REFERENCES (NONE)
1426 C***ROUTINES CALLED (NONE)
1427 C***REVISION HISTORY (YYMMDD)
1428 C 930915 DATE WRITTEN
1429 C 930922 Made user-callable, and other cosmetic changes. (FNF)
1430 C***END PROLOGUE IUMACH
1431 C
1432 C*Internal Notes:
1433 C The built-in value of 6 is standard on a wide range of Fortran
1434 C systems. This may be machine-dependent.
1435 C**End
1436 C***FIRST EXECUTABLE STATEMENT IUMACH
1437  iumach = 6
1438 C
1439  RETURN
1440 C----------------------- End of Function IUMACH ------------------------
1441  END