*DECK DGEFA
SUBROUTINE DGEFA (A, LDA, N, IPVT, INFO)
C***BEGIN PROLOGUE DGEFA
C***PURPOSE Factor a matrix using Gaussian elimination.
C***CATEGORY D2A1
C***TYPE DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C)
C***KEYWORDS GENERAL MATRIX, LINEAR ALGEBRA, LINPACK,
C MATRIX FACTORIZATION
C***AUTHOR Moler, C. B., (U. of New Mexico)
C***DESCRIPTION
C
C DGEFA factors a double precision matrix by Gaussian elimination.
C
C DGEFA is usually called by DGECO, but it can be called
C directly with a saving in time if RCOND is not needed.
C (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .
C
C On Entry
C
C A DOUBLE PRECISION(LDA, N)
C the matrix to be factored.
C
C LDA INTEGER
C the leading dimension of the array A .
C
C N INTEGER
C the order of the matrix A .
C
C On Return
C
C A an upper triangular matrix and the multipliers
C which were used to obtain it.
C The factorization can be written A = L*U where
C L is a product of permutation and unit lower
C triangular matrices and U is upper triangular.
C
C IPVT INTEGER(N)
C an integer vector of pivot indices.
C
C INFO INTEGER
C = 0 normal value.
C = K if U(K,K) .EQ. 0.0 . This is not an error
C condition for this subroutine, but it does
C indicate that DGESL or DGEDI will divide by zero
C if called. Use RCOND in DGECO for a reliable
C indication of singularity.
C
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX
C***REVISION HISTORY (YYMMDD)
C 780814 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DGEFA
INTEGER LDA,N,IPVT(*),INFO
DOUBLE PRECISION A(LDA,*)
C
DOUBLE PRECISION T
INTEGER IDAMAX,J,K,KP1,L,NM1
C
C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
C***FIRST EXECUTABLE STATEMENT DGEFA
INFO = 0
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 70
DO 60 K = 1, NM1
KP1 = K + 1
C
C FIND L = PIVOT INDEX
C
L = IDAMAX(N-K+1,A(K,K),1) + K - 1
IPVT(K) = L
C
C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
IF (A(L,K) .EQ. 0.0D0) GO TO 40
C
C INTERCHANGE IF NECESSARY
C
IF (L .EQ. K) GO TO 10
T = A(L,K)
A(L,K) = A(K,K)
A(K,K) = T
10 CONTINUE
C
C COMPUTE MULTIPLIERS
C
T = -1.0D0/A(K,K)
CALL DSCAL(N-K,T,A(K+1,K),1)
C
C ROW ELIMINATION WITH COLUMN INDEXING
C
DO 30 J = KP1, N
T = A(L,J)
IF (L .EQ. K) GO TO 20
A(L,J) = A(K,J)
A(K,J) = T
20 CONTINUE
CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
30 CONTINUE
GO TO 50
40 CONTINUE
INFO = K
50 CONTINUE
60 CONTINUE
70 CONTINUE
IPVT(N) = N
IF (A(N,N) .EQ. 0.0D0) INFO = N
RETURN
END
*DECK DGESL
SUBROUTINE DGESL (A, LDA, N, IPVT, B, JOB)
C***BEGIN PROLOGUE DGESL
C***PURPOSE Solve the real system A*X=B or TRANS(A)*X=B using the
C factors computed by DGECO or DGEFA.
C***CATEGORY D2A1
C***TYPE DOUBLE PRECISION (SGESL-S, DGESL-D, CGESL-C)
C***KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
C***AUTHOR Moler, C. B., (U. of New Mexico)
C***DESCRIPTION
C
C DGESL solves the double precision system
C A * X = B or TRANS(A) * X = B
C using the factors computed by DGECO or DGEFA.
C
C On Entry
C
C A DOUBLE PRECISION(LDA, N)
C the output from DGECO or DGEFA.
C
C LDA INTEGER
C the leading dimension of the array A .
C
C N INTEGER
C the order of the matrix A .
C
C IPVT INTEGER(N)
C the pivot vector from DGECO or DGEFA.
C
C B DOUBLE PRECISION(N)
C the right hand side vector.
C
C JOB INTEGER
C = 0 to solve A*X = B ,
C = nonzero to solve TRANS(A)*X = B where
C TRANS(A) is the transpose.
C
C On Return
C
C B the solution vector X .
C
C Error Condition
C
C A division by zero will occur if the input factor contains a
C zero on the diagonal. Technically this indicates singularity
C but it is often caused by improper arguments or improper
C setting of LDA . It will not occur if the subroutines are
C called correctly and if DGECO has set RCOND .GT. 0.0
C or DGEFA has set INFO .EQ. 0 .
C
C To compute INVERSE(A) * C where C is a matrix
C with P columns
C CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
C IF (RCOND is too small) GO TO ...
C DO 10 J = 1, P
C CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
C 10 CONTINUE
C
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED DAXPY, DDOT
C***REVISION HISTORY (YYMMDD)
C 780814 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DGESL
INTEGER LDA,N,IPVT(*),JOB
DOUBLE PRECISION A(LDA,*),B(*)
C
DOUBLE PRECISION DDOT,T
INTEGER K,KB,L,NM1
C***FIRST EXECUTABLE STATEMENT DGESL
NM1 = N - 1
IF (JOB .NE. 0) GO TO 50
C
C JOB = 0 , SOLVE A * X = B
C FIRST SOLVE L*Y = B
C
IF (NM1 .LT. 1) GO TO 30
DO 20 K = 1, NM1
L = IPVT(K)
T = B(L)
IF (L .EQ. K) GO TO 10
B(L) = B(K)
B(K) = T
10 CONTINUE
CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
20 CONTINUE
30 CONTINUE
C
C NOW SOLVE U*X = Y
C
DO 40 KB = 1, N
K = N + 1 - KB
B(K) = B(K)/A(K,K)
T = -B(K)
CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
40 CONTINUE
GO TO 100
50 CONTINUE
C
C JOB = NONZERO, SOLVE TRANS(A) * X = B
C FIRST SOLVE TRANS(U)*Y = B
C
DO 60 K = 1, N
T = DDOT(K-1,A(1,K),1,B(1),1)
B(K) = (B(K) - T)/A(K,K)
60 CONTINUE
C
C NOW SOLVE TRANS(L)*X = Y
C
IF (NM1 .LT. 1) GO TO 90
DO 80 KB = 1, NM1
K = N - KB
B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
L = IPVT(K)
IF (L .EQ. K) GO TO 70
T = B(L)
B(L) = B(K)
B(K) = T
70 CONTINUE
80 CONTINUE
90 CONTINUE
100 CONTINUE
RETURN
END
*DECK DGBFA
SUBROUTINE DGBFA (ABD, LDA, N, ML, MU, IPVT, INFO)
C***BEGIN PROLOGUE DGBFA
C***PURPOSE Factor a band matrix using Gaussian elimination.
C***CATEGORY D2A2
C***TYPE DOUBLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
C***AUTHOR Moler, C. B., (U. of New Mexico)
C***DESCRIPTION
C
C DGBFA factors a double precision band matrix by elimination.
C
C DGBFA is usually called by DGBCO, but it can be called
C directly with a saving in time if RCOND is not needed.
C
C On Entry
C
C ABD DOUBLE PRECISION(LDA, N)
C contains the matrix in band storage. The columns
C of the matrix are stored in the columns of ABD and
C the diagonals of the matrix are stored in rows
C ML+1 through 2*ML+MU+1 of ABD .
C See the comments below for details.
C
C LDA INTEGER
C the leading dimension of the array ABD .
C LDA must be .GE. 2*ML + MU + 1 .
C
C N INTEGER
C the order of the original matrix.
C
C ML INTEGER
C number of diagonals below the main diagonal.
C 0 .LE. ML .LT. N .
C
C MU INTEGER
C number of diagonals above the main diagonal.
C 0 .LE. MU .LT. N .
C More efficient if ML .LE. MU .
C On Return
C
C ABD an upper triangular matrix in band storage and
C the multipliers which were used to obtain it.
C The factorization can be written A = L*U where
C L is a product of permutation and unit lower
C triangular matrices and U is upper triangular.
C
C IPVT INTEGER(N)
C an integer vector of pivot indices.
C
C INFO INTEGER
C = 0 normal value.
C = K if U(K,K) .EQ. 0.0 . This is not an error
C condition for this subroutine, but it does
C indicate that DGBSL will divide by zero if
C called. Use RCOND in DGBCO for a reliable
C indication of singularity.
C
C Band Storage
C
C If A is a band matrix, the following program segment
C will set up the input.
C
C ML = (band width below the diagonal)
C MU = (band width above the diagonal)
C M = ML + MU + 1
C DO 20 J = 1, N
C I1 = MAX(1, J-MU)
C I2 = MIN(N, J+ML)
C DO 10 I = I1, I2
C K = I - J + M
C ABD(K,J) = A(I,J)
C 10 CONTINUE
C 20 CONTINUE
C
C This uses rows ML+1 through 2*ML+MU+1 of ABD .
C In addition, the first ML rows in ABD are used for
C elements generated during the triangularization.
C The total number of rows needed in ABD is 2*ML+MU+1 .
C The ML+MU by ML+MU upper left triangle and the
C ML by ML lower right triangle are not referenced.
C
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX
C***REVISION HISTORY (YYMMDD)
C 780814 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DGBFA
INTEGER LDA,N,ML,MU,IPVT(*),INFO
DOUBLE PRECISION ABD(LDA,*)
C
DOUBLE PRECISION T
INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
C
C***FIRST EXECUTABLE STATEMENT DGBFA
M = ML + MU + 1
INFO = 0
C
C ZERO INITIAL FILL-IN COLUMNS
C
J0 = MU + 2
J1 = MIN(N,M) - 1
IF (J1 .LT. J0) GO TO 30
DO 20 JZ = J0, J1
I0 = M + 1 - JZ
DO 10 I = I0, ML
ABD(I,JZ) = 0.0D0
10 CONTINUE
20 CONTINUE
30 CONTINUE
JZ = J1
JU = 0
C
C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 130
DO 120 K = 1, NM1
KP1 = K + 1
C
C ZERO NEXT FILL-IN COLUMN
C
JZ = JZ + 1
IF (JZ .GT. N) GO TO 50
IF (ML .LT. 1) GO TO 50
DO 40 I = 1, ML
ABD(I,JZ) = 0.0D0
40 CONTINUE
50 CONTINUE
C
C FIND L = PIVOT INDEX
C
LM = MIN(ML,N-K)
L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
IPVT(K) = L + K - M
C
C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
C
C INTERCHANGE IF NECESSARY
C
IF (L .EQ. M) GO TO 60
T = ABD(L,K)
ABD(L,K) = ABD(M,K)
ABD(M,K) = T
60 CONTINUE
C
C COMPUTE MULTIPLIERS
C
T = -1.0D0/ABD(M,K)
CALL DSCAL(LM,T,ABD(M+1,K),1)
C
C ROW ELIMINATION WITH COLUMN INDEXING
C
JU = MIN(MAX(JU,MU+IPVT(K)),N)
MM = M
IF (JU .LT. KP1) GO TO 90
DO 80 J = KP1, JU
L = L - 1
MM = MM - 1
T = ABD(L,J)
IF (L .EQ. MM) GO TO 70
ABD(L,J) = ABD(MM,J)
ABD(MM,J) = T
70 CONTINUE
CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
80 CONTINUE
90 CONTINUE
GO TO 110
100 CONTINUE
INFO = K
110 CONTINUE
120 CONTINUE
130 CONTINUE
IPVT(N) = N
IF (ABD(M,N) .EQ. 0.0D0) INFO = N
RETURN
END
*DECK DGBSL
SUBROUTINE DGBSL (ABD, LDA, N, ML, MU, IPVT, B, JOB)
C***BEGIN PROLOGUE DGBSL
C***PURPOSE Solve the real band system A*X=B or TRANS(A)*X=B using
C the factors computed by DGBCO or DGBFA.
C***CATEGORY D2A2
C***TYPE DOUBLE PRECISION (SGBSL-S, DGBSL-D, CGBSL-C)
C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
C***AUTHOR Moler, C. B., (U. of New Mexico)
C***DESCRIPTION
C
C DGBSL solves the double precision band system
C A * X = B or TRANS(A) * X = B
C using the factors computed by DGBCO or DGBFA.
C
C On Entry
C
C ABD DOUBLE PRECISION(LDA, N)
C the output from DGBCO or DGBFA.
C
C LDA INTEGER
C the leading dimension of the array ABD .
C
C N INTEGER
C the order of the original matrix.
C
C ML INTEGER
C number of diagonals below the main diagonal.
C
C MU INTEGER
C number of diagonals above the main diagonal.
C
C IPVT INTEGER(N)
C the pivot vector from DGBCO or DGBFA.
C
C B DOUBLE PRECISION(N)
C the right hand side vector.
C
C JOB INTEGER
C = 0 to solve A*X = B ,
C = nonzero to solve TRANS(A)*X = B , where
C TRANS(A) is the transpose.
C
C On Return
C
C B the solution vector X .
C
C Error Condition
C
C A division by zero will occur if the input factor contains a
C zero on the diagonal. Technically this indicates singularity
C but it is often caused by improper arguments or improper
C setting of LDA . It will not occur if the subroutines are
C called correctly and if DGBCO has set RCOND .GT. 0.0
C or DGBFA has set INFO .EQ. 0 .
C
C To compute INVERSE(A) * C where C is a matrix
C with P columns
C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
C IF (RCOND is too small) GO TO ...
C DO 10 J = 1, P
C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
C 10 CONTINUE
C
C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED DAXPY, DDOT
C***REVISION HISTORY (YYMMDD)
C 780814 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DGBSL
INTEGER LDA,N,ML,MU,IPVT(*),JOB
DOUBLE PRECISION ABD(LDA,*),B(*)
C
DOUBLE PRECISION DDOT,T
INTEGER K,KB,L,LA,LB,LM,M,NM1
C***FIRST EXECUTABLE STATEMENT DGBSL
M = MU + ML + 1
NM1 = N - 1
IF (JOB .NE. 0) GO TO 50
C
C JOB = 0 , SOLVE A * X = B
C FIRST SOLVE L*Y = B
C
IF (ML .EQ. 0) GO TO 30
IF (NM1 .LT. 1) GO TO 30
DO 20 K = 1, NM1
LM = MIN(ML,N-K)
L = IPVT(K)
T = B(L)
IF (L .EQ. K) GO TO 10
B(L) = B(K)
B(K) = T
10 CONTINUE
CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
20 CONTINUE
30 CONTINUE
C
C NOW SOLVE U*X = Y
C
DO 40 KB = 1, N
K = N + 1 - KB
B(K) = B(K)/ABD(M,K)
LM = MIN(K,M) - 1
LA = M - LM
LB = K - LM
T = -B(K)
CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
40 CONTINUE
GO TO 100
50 CONTINUE
C
C JOB = NONZERO, SOLVE TRANS(A) * X = B
C FIRST SOLVE TRANS(U)*Y = B
C
DO 60 K = 1, N
LM = MIN(K,M) - 1
LA = M - LM
LB = K - LM
T = DDOT(LM,ABD(LA,K),1,B(LB),1)
B(K) = (B(K) - T)/ABD(M,K)
60 CONTINUE
C
C NOW SOLVE TRANS(L)*X = Y
C
IF (ML .EQ. 0) GO TO 90
IF (NM1 .LT. 1) GO TO 90
DO 80 KB = 1, NM1
K = N - KB
LM = MIN(ML,N-K)
B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
L = IPVT(K)
IF (L .EQ. K) GO TO 70
T = B(L)
B(L) = B(K)
B(K) = T
70 CONTINUE
80 CONTINUE
90 CONTINUE
100 CONTINUE
RETURN
END
*DECK DAXPY
SUBROUTINE DAXPY (N, DA, DX, INCX, DY, INCY)
C***BEGIN PROLOGUE DAXPY
C***PURPOSE Compute a constant times a vector plus a vector.
C***CATEGORY D1A7
C***TYPE DOUBLE PRECISION (SAXPY-S, DAXPY-D, CAXPY-C)
C***KEYWORDS BLAS, LINEAR ALGEBRA, TRIAD, VECTOR
C***AUTHOR Lawson, C. L., (JPL)
C Hanson, R. J., (SNLA)
C Kincaid, D. R., (U. of Texas)
C Krogh, F. T., (JPL)
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C DA double precision scalar multiplier
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C DY double precision vector with N elements
C INCY storage spacing between elements of DY
C
C --Output--
C DY double precision result (unchanged if N .LE. 0)
C
C Overwrite double precision DY with double precision DA*DX + DY.
C For I = 0 to N-1, replace DY(LY+I*INCY) with DA*DX(LX+I*INCX) +
C DY(LY+I*INCY),
C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
C defined in a similar way using INCY.
C
C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C Krogh, Basic linear algebra subprograms for Fortran
C usage, Algorithm No. 539, Transactions on Mathematical
C Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791001 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DAXPY
DOUBLE PRECISION DX(*), DY(*), DA
C***FIRST EXECUTABLE STATEMENT DAXPY
IF (N.LE.0 .OR. DA.EQ.0.0D0) RETURN
IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
C
C Code for unequal or nonpositive increments.
C
5 IX = 1
IY = 1
IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DY(IY) = DY(IY) + DA*DX(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C Code for both increments equal to 1.
C
C Clean-up loop so remaining vector length is a multiple of 4.
C
20 M = MOD(N,4)
IF (M .EQ. 0) GO TO 40
DO 30 I = 1,M
DY(I) = DY(I) + DA*DX(I)
30 CONTINUE
IF (N .LT. 4) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,4
DY(I) = DY(I) + DA*DX(I)
DY(I+1) = DY(I+1) + DA*DX(I+1)
DY(I+2) = DY(I+2) + DA*DX(I+2)
DY(I+3) = DY(I+3) + DA*DX(I+3)
50 CONTINUE
RETURN
C
C Code for equal, positive, non-unit increments.
C
60 NS = N*INCX
DO 70 I = 1,NS,INCX
DY(I) = DA*DX(I) + DY(I)
70 CONTINUE
RETURN
END
*DECK DCOPY
SUBROUTINE DCOPY (N, DX, INCX, DY, INCY)
C***BEGIN PROLOGUE DCOPY
C***PURPOSE Copy a vector.
C***CATEGORY D1A5
C***TYPE DOUBLE PRECISION (SCOPY-S, DCOPY-D, CCOPY-C, ICOPY-I)
C***KEYWORDS BLAS, COPY, LINEAR ALGEBRA, VECTOR
C***AUTHOR Lawson, C. L., (JPL)
C Hanson, R. J., (SNLA)
C Kincaid, D. R., (U. of Texas)
C Krogh, F. T., (JPL)
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C DY double precision vector with N elements
C INCY storage spacing between elements of DY
C
C --Output--
C DY copy of vector DX (unchanged if N .LE. 0)
C
C Copy double precision DX to double precision DY.
C For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY),
C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
C defined in a similar way using INCY.
C
C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C Krogh, Basic linear algebra subprograms for Fortran
C usage, Algorithm No. 539, Transactions on Mathematical
C Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791001 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DCOPY
DOUBLE PRECISION DX(*), DY(*)
C***FIRST EXECUTABLE STATEMENT DCOPY
IF (N .LE. 0) RETURN
IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
C
C Code for unequal or nonpositive increments.
C
5 IX = 1
IY = 1
IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DY(IY) = DX(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C Code for both increments equal to 1.
C
C Clean-up loop so remaining vector length is a multiple of 7.
C
20 M = MOD(N,7)
IF (M .EQ. 0) GO TO 40
DO 30 I = 1,M
DY(I) = DX(I)
30 CONTINUE
IF (N .LT. 7) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,7
DY(I) = DX(I)
DY(I+1) = DX(I+1)
DY(I+2) = DX(I+2)
DY(I+3) = DX(I+3)
DY(I+4) = DX(I+4)
DY(I+5) = DX(I+5)
DY(I+6) = DX(I+6)
50 CONTINUE
RETURN
C
C Code for equal, positive, non-unit increments.
C
60 NS = N*INCX
DO 70 I = 1,NS,INCX
DY(I) = DX(I)
70 CONTINUE
RETURN
END
*DECK DDOT
DOUBLE PRECISION FUNCTION DDOT (N, DX, INCX, DY, INCY)
C***BEGIN PROLOGUE DDOT
C***PURPOSE Compute the inner product of two vectors.
C***CATEGORY D1A4
C***TYPE DOUBLE PRECISION (SDOT-S, DDOT-D, CDOTU-C)
C***KEYWORDS BLAS, INNER PRODUCT, LINEAR ALGEBRA, VECTOR
C***AUTHOR Lawson, C. L., (JPL)
C Hanson, R. J., (SNLA)
C Kincaid, D. R., (U. of Texas)
C Krogh, F. T., (JPL)
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C DY double precision vector with N elements
C INCY storage spacing between elements of DY
C
C --Output--
C DDOT double precision dot product (zero if N .LE. 0)
C
C Returns the dot product of double precision DX and DY.
C DDOT = sum for I = 0 to N-1 of DX(LX+I*INCX) * DY(LY+I*INCY),
C where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
C defined in a similar way using INCY.
C
C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C Krogh, Basic linear algebra subprograms for Fortran
C usage, Algorithm No. 539, Transactions on Mathematical
C Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791001 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920310 Corrected definition of LX in DESCRIPTION. (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DDOT
DOUBLE PRECISION DX(*), DY(*)
C***FIRST EXECUTABLE STATEMENT DDOT
DDOT = 0.0D0
IF (N .LE. 0) RETURN
IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
C
C Code for unequal or nonpositive increments.
C
5 IX = 1
IY = 1
IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
DO 10 I = 1,N
DDOT = DDOT + DX(IX)*DY(IY)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
C
C Code for both increments equal to 1.
C
C Clean-up loop so remaining vector length is a multiple of 5.
C
20 M = MOD(N,5)
IF (M .EQ. 0) GO TO 40
DO 30 I = 1,M
DDOT = DDOT + DX(I)*DY(I)
30 CONTINUE
IF (N .LT. 5) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,5
DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + DX(I+2)*DY(I+2) +
1 DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
50 CONTINUE
RETURN
C
C Code for equal, positive, non-unit increments.
C
60 NS = N*INCX
DO 70 I = 1,NS,INCX
DDOT = DDOT + DX(I)*DY(I)
70 CONTINUE
RETURN
END
*DECK DNRM2
DOUBLE PRECISION FUNCTION DNRM2 (N, DX, INCX)
C***BEGIN PROLOGUE DNRM2
C***PURPOSE Compute the Euclidean length (L2 norm) of a vector.
C***CATEGORY D1A3B
C***TYPE DOUBLE PRECISION (SNRM2-S, DNRM2-D, SCNRM2-C)
C***KEYWORDS BLAS, EUCLIDEAN LENGTH, EUCLIDEAN NORM, L2,
C LINEAR ALGEBRA, UNITARY, VECTOR
C***AUTHOR Lawson, C. L., (JPL)
C Hanson, R. J., (SNLA)
C Kincaid, D. R., (U. of Texas)
C Krogh, F. T., (JPL)
C***DESCRIPTION
C
C B L A S Subprogram
C Description of parameters
C
C --Input--
C N number of elements in input vector(s)
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C
C --Output--
C DNRM2 double precision result (zero if N .LE. 0)
C
C Euclidean norm of the N-vector stored in DX with storage
C increment INCX.
C If N .LE. 0, return with result = 0.
C If N .GE. 1, then INCX must be .GE. 1
C
C Four phase method using two built-in constants that are
C hopefully applicable to all machines.
C CUTLO = maximum of SQRT(U/EPS) over all known machines.
C CUTHI = minimum of SQRT(V) over all known machines.
C where
C EPS = smallest no. such that EPS + 1. .GT. 1.
C U = smallest positive no. (underflow limit)
C V = largest no. (overflow limit)
C
C Brief outline of algorithm.
C
C Phase 1 scans zero components.
C move to phase 2 when a component is nonzero and .LE. CUTLO
C move to phase 3 when a component is .GT. CUTLO
C move to phase 4 when a component is .GE. CUTHI/M
C where M = N for X() real and M = 2*N for complex.
C
C Values for CUTLO and CUTHI.
C From the environmental parameters listed in the IMSL converter
C document the limiting values are as follows:
C CUTLO, S.P. U/EPS = 2**(-102) for Honeywell. Close seconds are
C Univac and DEC at 2**(-103)
C Thus CUTLO = 2**(-51) = 4.44089E-16
C CUTHI, S.P. V = 2**127 for Univac, Honeywell, and DEC.
C Thus CUTHI = 2**(63.5) = 1.30438E19
C CUTLO, D.P. U/EPS = 2**(-67) for Honeywell and DEC.
C Thus CUTLO = 2**(-33.5) = 8.23181D-11
C CUTHI, D.P. same as S.P. CUTHI = 1.30438D19
C DATA CUTLO, CUTHI /8.232D-11, 1.304D19/
C DATA CUTLO, CUTHI /4.441E-16, 1.304E19/
C
C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C Krogh, Basic linear algebra subprograms for Fortran
C usage, Algorithm No. 539, Transactions on Mathematical
C Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791001 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DNRM2
INTEGER NEXT
DOUBLE PRECISION DX(*), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO,
+ ONE
SAVE CUTLO, CUTHI, ZERO, ONE
DATA ZERO, ONE /0.0D0, 1.0D0/
C
DATA CUTLO, CUTHI /8.232D-11, 1.304D19/
C***FIRST EXECUTABLE STATEMENT DNRM2
IF (N .GT. 0) GO TO 10
DNRM2 = ZERO
GO TO 300
C
10 ASSIGN 30 TO NEXT
SUM = ZERO
NN = N * INCX
C
C BEGIN MAIN LOOP
C
I = 1
20 GO TO NEXT,(30, 50, 70, 110)
30 IF (ABS(DX(I)) .GT. CUTLO) GO TO 85
ASSIGN 50 TO NEXT
XMAX = ZERO
C
C PHASE 1. SUM IS ZERO
C
50 IF (DX(I) .EQ. ZERO) GO TO 200
IF (ABS(DX(I)) .GT. CUTLO) GO TO 85
C
C PREPARE FOR PHASE 2.
C
ASSIGN 70 TO NEXT
GO TO 105
C
C PREPARE FOR PHASE 4.
C
100 I = J
ASSIGN 110 TO NEXT
SUM = (SUM / DX(I)) / DX(I)
105 XMAX = ABS(DX(I))
GO TO 115
C
C PHASE 2. SUM IS SMALL.
C SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
70 IF (ABS(DX(I)) .GT. CUTLO) GO TO 75
C
C COMMON CODE FOR PHASES 2 AND 4.
C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
C
110 IF (ABS(DX(I)) .LE. XMAX) GO TO 115
SUM = ONE + SUM * (XMAX / DX(I))**2
XMAX = ABS(DX(I))
GO TO 200
C
115 SUM = SUM + (DX(I)/XMAX)**2
GO TO 200
C
C PREPARE FOR PHASE 3.
C
75 SUM = (SUM * XMAX) * XMAX
C
C FOR REAL OR D.P. SET HITEST = CUTHI/N
C FOR COMPLEX SET HITEST = CUTHI/(2*N)
C
85 HITEST = CUTHI / N
C
C PHASE 3. SUM IS MID-RANGE. NO SCALING.
C
DO 95 J = I,NN,INCX
IF (ABS(DX(J)) .GE. HITEST) GO TO 100
95 SUM = SUM + DX(J)**2
DNRM2 = SQRT(SUM)
GO TO 300
C
200 CONTINUE
I = I + INCX
IF (I .LE. NN) GO TO 20
C
C END OF MAIN LOOP.
C
C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
DNRM2 = XMAX * SQRT(SUM)
300 CONTINUE
RETURN
END
*DECK DSCAL
SUBROUTINE DSCAL (N, DA, DX, INCX)
C***BEGIN PROLOGUE DSCAL
C***PURPOSE Multiply a vector by a constant.
C***CATEGORY D1A6
C***TYPE DOUBLE PRECISION (SSCAL-S, DSCAL-D, CSCAL-C)
C***KEYWORDS BLAS, LINEAR ALGEBRA, SCALE, VECTOR
C***AUTHOR Lawson, C. L., (JPL)
C Hanson, R. J., (SNLA)
C Kincaid, D. R., (U. of Texas)
C Krogh, F. T., (JPL)
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C DA double precision scale factor
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C
C --Output--
C DX double precision result (unchanged if N.LE.0)
C
C Replace double precision DX by double precision DA*DX.
C For I = 0 to N-1, replace DX(IX+I*INCX) with DA * DX(IX+I*INCX),
C where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX.
C
C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C Krogh, Basic linear algebra subprograms for Fortran
C usage, Algorithm No. 539, Transactions on Mathematical
C Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791001 DATE WRITTEN
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900821 Modified to correct problem with a negative increment.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DSCAL
DOUBLE PRECISION DA, DX(*)
INTEGER I, INCX, IX, M, MP1, N
C***FIRST EXECUTABLE STATEMENT DSCAL
IF (N .LE. 0) RETURN
IF (INCX .EQ. 1) GOTO 20
C
C Code for increment not equal to 1.
C
IX = 1
IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
DO 10 I = 1,N
DX(IX) = DA*DX(IX)
IX = IX + INCX
10 CONTINUE
RETURN
C
C Code for increment equal to 1.
C
C Clean-up loop so remaining vector length is a multiple of 5.
C
20 M = MOD(N,5)
IF (M .EQ. 0) GOTO 40
DO 30 I = 1,M
DX(I) = DA*DX(I)
30 CONTINUE
IF (N .LT. 5) RETURN
40 MP1 = M + 1
DO 50 I = MP1,N,5
DX(I) = DA*DX(I)
DX(I+1) = DA*DX(I+1)
DX(I+2) = DA*DX(I+2)
DX(I+3) = DA*DX(I+3)
DX(I+4) = DA*DX(I+4)
50 CONTINUE
RETURN
END
*DECK IDAMAX
INTEGER FUNCTION IDAMAX (N, DX, INCX)
C***BEGIN PROLOGUE IDAMAX
C***PURPOSE Find the smallest index of that component of a vector
C having the maximum magnitude.
C***CATEGORY D1A2
C***TYPE DOUBLE PRECISION (ISAMAX-S, IDAMAX-D, ICAMAX-C)
C***KEYWORDS BLAS, LINEAR ALGEBRA, MAXIMUM COMPONENT, VECTOR
C***AUTHOR Lawson, C. L., (JPL)
C Hanson, R. J., (SNLA)
C Kincaid, D. R., (U. of Texas)
C Krogh, F. T., (JPL)
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C DX double precision vector with N elements
C INCX storage spacing between elements of DX
C
C --Output--
C IDAMAX smallest index (zero if N .LE. 0)
C
C Find smallest index of maximum magnitude of double precision DX.
C IDAMAX = first I, I = 1 to N, to maximize ABS(DX(IX+(I-1)*INCX)),
C where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX.
C
C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C Krogh, Basic linear algebra subprograms for Fortran
C usage, Algorithm No. 539, Transactions on Mathematical
C Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791001 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890531 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900821 Modified to correct problem with a negative increment.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE IDAMAX
DOUBLE PRECISION DX(*), DMAX, XMAG
INTEGER I, INCX, IX, N
C***FIRST EXECUTABLE STATEMENT IDAMAX
IDAMAX = 0
IF (N .LE. 0) RETURN
IDAMAX = 1
IF (N .EQ. 1) RETURN
C
IF (INCX .EQ. 1) GOTO 20
C
C Code for increments not equal to 1.
C
IX = 1
IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
DMAX = ABS(DX(IX))
IX = IX + INCX
DO 10 I = 2,N
XMAG = ABS(DX(IX))
IF (XMAG .GT. DMAX) THEN
IDAMAX = I
DMAX = XMAG
ENDIF
IX = IX + INCX
10 CONTINUE
RETURN
C
C Code for increments equal to 1.
C
20 DMAX = ABS(DX(1))
DO 30 I = 2,N
XMAG = ABS(DX(I))
IF (XMAG .GT. DMAX) THEN
IDAMAX = I
DMAX = XMAG
ENDIF
30 CONTINUE
RETURN
END
*DECK XERRWD
SUBROUTINE XERRWD (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
C***BEGIN PROLOGUE XERRWD
C***SUBSIDIARY
C***PURPOSE Write error message with values.
C***CATEGORY R3C
C***TYPE DOUBLE PRECISION (XERRWV-S, XERRWD-D)
C***AUTHOR Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV,
C as given here, constitute a simplified version of the SLATEC error
C handling package.
C
C All arguments are input arguments.
C
C MSG = The message (character array).
C NMES = The length of MSG (number of characters).
C NERR = The error number (not used).
C LEVEL = The error level..
C 0 or 1 means recoverable (control returns to caller).
C 2 means fatal (run is aborted--see note below).
C NI = Number of integers (0, 1, or 2) to be printed with message.
C I1,I2 = Integers to be printed, depending on NI.
C NR = Number of reals (0, 1, or 2) to be printed with message.
C R1,R2 = Reals to be printed, depending on NR.
C
C Note.. this routine is machine-dependent and specialized for use
C in limited context, in the following ways..
C 1. The argument MSG is assumed to be of type CHARACTER, and
C the message is printed with a format of (1X,A).
C 2. The message is assumed to take only one line.
C Multi-line messages are generated by repeated calls.
C 3. If LEVEL = 2, control passes to the statement STOP
C to abort the run. This statement may be machine-dependent.
C 4. R1 and R2 are assumed to be in double precision and are printed
C in D21.13 format.
C
C***ROUTINES CALLED IXSAV
C***REVISION HISTORY (YYMMDD)
C 920831 DATE WRITTEN
C 921118 Replaced MFLGSV/LUNSAV by IXSAV. (ACH)
C 930329 Modified prologue to SLATEC format. (FNF)
C 930407 Changed MSG from CHARACTER*1 array to variable. (FNF)
C 930922 Minor cosmetic change. (FNF)
C***END PROLOGUE XERRWD
C
C*Internal Notes:
C
C For a different default logical unit number, IXSAV (or a subsidiary
C routine that it calls) will need to be modified.
C For a different run-abort command, change the statement following
C statement 100 at the end.
C-----------------------------------------------------------------------
C Subroutines called by XERRWD.. None
C Function routine called by XERRWD.. IXSAV
C-----------------------------------------------------------------------
C**End
C
C Declare arguments.
C
DOUBLE PRECISION R1, R2
INTEGER NMES, NERR, LEVEL, NI, I1, I2, NR
CHARACTER*(*) MSG
C
C Declare local variables.
C
INTEGER LUNIT, IXSAV, MESFLG
C
C Get logical unit number and message print flag.
C
C***FIRST EXECUTABLE STATEMENT XERRWD
LUNIT = IXSAV (1, 0, .FALSE.)
MESFLG = IXSAV (2, 0, .FALSE.)
IF (MESFLG .EQ. 0) GO TO 100
C
C Write the message.
C
WRITE (LUNIT,10) MSG
10 FORMAT(1X,A)
IF (NI .EQ. 1) WRITE (LUNIT, 20) I1
20 FORMAT(6X,'In above message, I1 =',I10)
IF (NI .EQ. 2) WRITE (LUNIT, 30) I1,I2
30 FORMAT(6X,'In above message, I1 =',I10,3X,'I2 =',I10)
IF (NR .EQ. 1) WRITE (LUNIT, 40) R1
40 FORMAT(6X,'In above message, R1 =',D21.13)
IF (NR .EQ. 2) WRITE (LUNIT, 50) R1,R2
50 FORMAT(6X,'In above, R1 =',D21.13,3X,'R2 =',D21.13)
C
C Abort the run if LEVEL = 2.
C
100 IF (LEVEL .NE. 2) RETURN
STOP
C----------------------- End of Subroutine XERRWD ----------------------
END
*DECK XSETF
SUBROUTINE XSETF (MFLAG)
C***BEGIN PROLOGUE XSETF
C***PURPOSE Reset the error print control flag.
C***CATEGORY R3A
C***TYPE ALL (XSETF-A)
C***KEYWORDS ERROR CONTROL
C***AUTHOR Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C XSETF sets the error print control flag to MFLAG:
C MFLAG=1 means print all messages (the default).
C MFLAG=0 means no printing.
C
C***SEE ALSO XERRWD, XERRWV
C***REFERENCES (NONE)
C***ROUTINES CALLED IXSAV
C***REVISION HISTORY (YYMMDD)
C 921118 DATE WRITTEN
C 930329 Added SLATEC format prologue. (FNF)
C 930407 Corrected SEE ALSO section. (FNF)
C 930922 Made user-callable, and other cosmetic changes. (FNF)
C***END PROLOGUE XSETF
C
C Subroutines called by XSETF.. None
C Function routine called by XSETF.. IXSAV
C-----------------------------------------------------------------------
C**End
INTEGER MFLAG, JUNK, IXSAV
C
C***FIRST EXECUTABLE STATEMENT XSETF
IF (MFLAG .EQ. 0 .OR. MFLAG .EQ. 1) JUNK = IXSAV (2,MFLAG,.TRUE.)
RETURN
C----------------------- End of Subroutine XSETF -----------------------
END
*DECK XSETUN
SUBROUTINE XSETUN (LUN)
C***BEGIN PROLOGUE XSETUN
C***PURPOSE Reset the logical unit number for error messages.
C***CATEGORY R3B
C***TYPE ALL (XSETUN-A)
C***KEYWORDS ERROR CONTROL
C***DESCRIPTION
C
C XSETUN sets the logical unit number for error messages to LUN.
C
C***AUTHOR Hindmarsh, Alan C., (LLNL)
C***SEE ALSO XERRWD, XERRWV
C***REFERENCES (NONE)
C***ROUTINES CALLED IXSAV
C***REVISION HISTORY (YYMMDD)
C 921118 DATE WRITTEN
C 930329 Added SLATEC format prologue. (FNF)
C 930407 Corrected SEE ALSO section. (FNF)
C 930922 Made user-callable, and other cosmetic changes. (FNF)
C***END PROLOGUE XSETUN
C
C Subroutines called by XSETUN.. None
C Function routine called by XSETUN.. IXSAV
C-----------------------------------------------------------------------
C**End
INTEGER LUN, JUNK, IXSAV
C
C***FIRST EXECUTABLE STATEMENT XSETUN
IF (LUN .GT. 0) JUNK = IXSAV (1,LUN,.TRUE.)
RETURN
C----------------------- End of Subroutine XSETUN ----------------------
END
*DECK IXSAV
INTEGER FUNCTION IXSAV (IPAR, IVALUE, ISET)
C***BEGIN PROLOGUE IXSAV
C***SUBSIDIARY
C***PURPOSE Save and recall error message control parameters.
C***CATEGORY R3C
C***TYPE ALL (IXSAV-A)
C***AUTHOR Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C IXSAV saves and recalls one of two error message parameters:
C LUNIT, the logical unit number to which messages are printed, and
C MESFLG, the message print flag.
C This is a modification of the SLATEC library routine J4SAVE.
C
C Saved local variables..
C LUNIT = Logical unit number for messages. The default is obtained
C by a call to IUMACH (may be machine-dependent).
C MESFLG = Print control flag..
C 1 means print all messages (the default).
C 0 means no printing.
C
C On input..
C IPAR = Parameter indicator (1 for LUNIT, 2 for MESFLG).
C IVALUE = The value to be set for the parameter, if ISET = .TRUE.
C ISET = Logical flag to indicate whether to read or write.
C If ISET = .TRUE., the parameter will be given
C the value IVALUE. If ISET = .FALSE., the parameter
C will be unchanged, and IVALUE is a dummy argument.
C
C On return..
C IXSAV = The (old) value of the parameter.
C
C***SEE ALSO XERRWD, XERRWV
C***ROUTINES CALLED IUMACH
C***REVISION HISTORY (YYMMDD)
C 921118 DATE WRITTEN
C 930329 Modified prologue to SLATEC format. (FNF)
C 930915 Added IUMACH call to get default output unit. (ACH)
C 930922 Minor cosmetic changes. (FNF)
C 010425 Type declaration for IUMACH added. (ACH)
C***END PROLOGUE IXSAV
C
C Subroutines called by IXSAV.. None
C Function routine called by IXSAV.. IUMACH
C-----------------------------------------------------------------------
C**End
LOGICAL ISET
INTEGER IPAR, IVALUE
C-----------------------------------------------------------------------
INTEGER IUMACH, LUNIT, MESFLG
C-----------------------------------------------------------------------
C The following Fortran-77 declaration is to cause the values of the
C listed (local) variables to be saved between calls to this routine.
C-----------------------------------------------------------------------
SAVE LUNIT, MESFLG
DATA LUNIT/-1/, MESFLG/1/
C
C***FIRST EXECUTABLE STATEMENT IXSAV
IF (IPAR .EQ. 1) THEN
IF (LUNIT .EQ. -1) LUNIT = IUMACH()
IXSAV = LUNIT
IF (ISET) LUNIT = IVALUE
ENDIF
C
IF (IPAR .EQ. 2) THEN
IXSAV = MESFLG
IF (ISET) MESFLG = IVALUE
ENDIF
C
RETURN
C----------------------- End of Function IXSAV -------------------------
END
*DECK IUMACH
INTEGER FUNCTION IUMACH()
C***BEGIN PROLOGUE IUMACH
C***PURPOSE Provide standard output unit number.
C***CATEGORY R1
C***TYPE INTEGER (IUMACH-I)
C***KEYWORDS MACHINE CONSTANTS
C***AUTHOR Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C *Usage:
C INTEGER LOUT, IUMACH
C LOUT = IUMACH()
C
C *Function Return Values:
C LOUT : the standard logical unit for Fortran output.
C
C***REFERENCES (NONE)
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 930915 DATE WRITTEN
C 930922 Made user-callable, and other cosmetic changes. (FNF)
C***END PROLOGUE IUMACH
C
C*Internal Notes:
C The built-in value of 6 is standard on a wide range of Fortran
C systems. This may be machine-dependent.
C**End
C***FIRST EXECUTABLE STATEMENT IUMACH
IUMACH = 6
C
RETURN
C----------------------- End of Function IUMACH ------------------------
END