*DECK SGEFA SUBROUTINE SGEFA (A, LDA, N, IPVT, INFO) C***BEGIN PROLOGUE SGEFA C***PURPOSE Factor a matrix using Gaussian elimination. C***CATEGORY D2A1 C***TYPE SINGLE 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 SGEFA factors a real matrix by Gaussian elimination. C C SGEFA is usually called by SGECO, but it can be called C directly with a saving in time if RCOND is not needed. C (Time for SGECO) = (1 + 9/N)*(Time for SGEFA) . C C On Entry C C A REAL(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 SGESL or SGEDI will divide by zero C if called. Use RCOND in SGECO 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 ISAMAX, SAXPY, SSCAL 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 SGEFA INTEGER LDA,N,IPVT(*),INFO REAL A(LDA,*) C REAL T INTEGER ISAMAX,J,K,KP1,L,NM1 C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C C***FIRST EXECUTABLE STATEMENT SGEFA 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 = ISAMAX(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.0E0) 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.0E0/A(K,K) CALL SSCAL(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 SAXPY(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.0E0) INFO = N RETURN END *DECK SGESL SUBROUTINE SGESL (A, LDA, N, IPVT, B, JOB) C***BEGIN PROLOGUE SGESL C***PURPOSE Solve the real system A*X=B or TRANS(A)*X=B using the C factors of SGECO or SGEFA. C***CATEGORY D2A1 C***TYPE SINGLE 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 SGESL solves the real system C A * X = B or TRANS(A) * X = B C using the factors computed by SGECO or SGEFA. C C On Entry C C A REAL(LDA, N) C the output from SGECO or SGEFA. 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 SGECO or SGEFA. C C B REAL(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 SGECO has set RCOND .GT. 0.0 C or SGEFA has set INFO .EQ. 0 . C C To compute INVERSE(A) * C where C is a matrix C with P columns C CALL SGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND is too small) GO TO ... C DO 10 J = 1, P C CALL SGESL(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 SAXPY, SDOT 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 SGESL INTEGER LDA,N,IPVT(*),JOB REAL A(LDA,*),B(*) C REAL SDOT,T INTEGER K,KB,L,NM1 C***FIRST EXECUTABLE STATEMENT SGESL 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 SAXPY(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 SAXPY(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 = SDOT(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) + SDOT(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 SGBFA SUBROUTINE SGBFA (ABD, LDA, N, ML, MU, IPVT, INFO) C***BEGIN PROLOGUE SGBFA C***PURPOSE Factor a band matrix using Gaussian elimination. C***CATEGORY D2A2 C***TYPE SINGLE 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 SGBFA factors a real band matrix by elimination. C C SGBFA is usually called by SBGCO, but it can be called C directly with a saving in time if RCOND is not needed. C C On Entry C C ABD REAL(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 SGBSL will divide by zero if C called. Use RCOND in SBGCO 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 ISAMAX, SAXPY, SSCAL 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 SGBFA INTEGER LDA,N,ML,MU,IPVT(*),INFO REAL ABD(LDA,*) C REAL T INTEGER I,ISAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1 C C***FIRST EXECUTABLE STATEMENT SGBFA 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.0E0 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.0E0 40 CONTINUE 50 CONTINUE C C FIND L = PIVOT INDEX C LM = MIN(ML,N-K) L = ISAMAX(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.0E0) 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.0E0/ABD(M,K) CALL SSCAL(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 SAXPY(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.0E0) INFO = N RETURN END *DECK SGBSL SUBROUTINE SGBSL (ABD, LDA, N, ML, MU, IPVT, B, JOB) C***BEGIN PROLOGUE SGBSL C***PURPOSE Solve the real band system A*X=B or TRANS(A)*X=B using C the factors computed by SGBCO or SGBFA. C***CATEGORY D2A2 C***TYPE SINGLE 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 SGBSL solves the real band system C A * X = B or TRANS(A) * X = B C using the factors computed by SBGCO or SGBFA. C C On Entry C C ABD REAL(LDA, N) C the output from SBGCO or SGBFA. 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 SBGCO or SGBFA. C C B REAL(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 SBGCO has set RCOND .GT. 0.0 C or SGBFA has set INFO .EQ. 0 . C C To compute INVERSE(A) * C where C is a matrix C with P columns C CALL SBGCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) C If (RCOND is too small) GO TO ... C DO 10 J = 1, P C CALL SGBSL(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 SAXPY, SDOT 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 SGBSL INTEGER LDA,N,ML,MU,IPVT(*),JOB REAL ABD(LDA,*),B(*) C REAL SDOT,T INTEGER K,KB,L,LA,LB,LM,M,NM1 C***FIRST EXECUTABLE STATEMENT SGBSL 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 SAXPY(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 SAXPY(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 = SDOT(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) + SDOT(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 SAXPY SUBROUTINE SAXPY (N, SA, SX, INCX, SY, INCY) C***BEGIN PROLOGUE SAXPY C***PURPOSE Compute a constant times a vector plus a vector. C***CATEGORY D1A7 C***TYPE SINGLE 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 SA single precision scalar multiplier C SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C SY single precision result (unchanged if N .LE. 0) C C Overwrite single precision SY with single precision SA*SX +SY. C For I = 0 to N-1, replace SY(LY+I*INCY) with SA*SX(LX+I*INCX) + C SY(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 SAXPY REAL SX(*), SY(*), SA C***FIRST EXECUTABLE STATEMENT SAXPY IF (N.LE.0 .OR. SA.EQ.0.0E0) 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 SY(IY) = SY(IY) + SA*SX(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 SY(I) = SY(I) + SA*SX(I) 30 CONTINUE IF (N .LT. 4) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I+1) = SY(I+1) + SA*SX(I+1) SY(I+2) = SY(I+2) + SA*SX(I+2) SY(I+3) = SY(I+3) + SA*SX(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 SY(I) = SA*SX(I) + SY(I) 70 CONTINUE RETURN END *DECK SCOPY SUBROUTINE SCOPY (N, SX, INCX, SY, INCY) C***BEGIN PROLOGUE SCOPY C***PURPOSE Copy a vector. C***CATEGORY D1A5 C***TYPE SINGLE 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 SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C SY copy of vector SX (unchanged if N .LE. 0) C C Copy single precision SX to single precision SY. C For I = 0 to N-1, copy SX(LX+I*INCX) to SY(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 SCOPY REAL SX(*), SY(*) C***FIRST EXECUTABLE STATEMENT SCOPY 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 SY(IY) = SX(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 SY(I) = SX(I) 30 CONTINUE IF (N .LT. 7) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 SY(I) = SX(I) SY(I+1) = SX(I+1) SY(I+2) = SX(I+2) SY(I+3) = SX(I+3) SY(I+4) = SX(I+4) SY(I+5) = SX(I+5) SY(I+6) = SX(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 SY(I) = SX(I) 70 CONTINUE RETURN END *DECK SDOT REAL FUNCTION SDOT (N, SX, INCX, SY, INCY) C***BEGIN PROLOGUE SDOT C***PURPOSE Compute the inner product of two vectors. C***CATEGORY D1A4 C***TYPE SINGLE 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 SX single precision vector with N elements C INCX storage spacing between elements of SX C SY single precision vector with N elements C INCY storage spacing between elements of SY C C --Output-- C SDOT single precision dot product (zero if N .LE. 0) C C Returns the dot product of single precision SX and SY. C SDOT = sum for I = 0 to N-1 of SX(LX+I*INCX) * SY(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 SDOT REAL SX(*), SY(*) C***FIRST EXECUTABLE STATEMENT SDOT SDOT = 0.0E0 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 SDOT = SDOT + SX(IX)*SY(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 SDOT = SDOT + SX(I)*SY(I) 30 CONTINUE IF (N .LT. 5) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SDOT = SDOT + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2) + 1 SX(I+3)*SY(I+3) + SX(I+4)*SY(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 SDOT = SDOT + SX(I)*SY(I) 70 CONTINUE RETURN END *DECK SNRM2 REAL FUNCTION SNRM2 (N, SX, INCX) C***BEGIN PROLOGUE SNRM2 C***PURPOSE Compute the Euclidean length (L2 norm) of a vector. C***CATEGORY D1A3B C***TYPE SINGLE 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 SX single precision vector with N elements C INCX storage spacing between elements of SX C C --Output-- C SNRM2 single precision result (zero if N .LE. 0) C C Euclidean norm of the N-vector stored in SX 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 SNRM2 INTEGER NEXT REAL SX(*), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE SAVE CUTLO, CUTHI, ZERO, ONE DATA ZERO, ONE /0.0E0, 1.0E0/ C DATA CUTLO, CUTHI /4.441E-16, 1.304E19/ C***FIRST EXECUTABLE STATEMENT SNRM2 IF (N .GT. 0) GO TO 10 SNRM2 = 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(SX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF (SX(I) .EQ. ZERO) GO TO 200 IF (ABS(SX(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 / SX(I)) / SX(I) 105 XMAX = ABS(SX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF (ABS(SX(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(SX(I)) .LE. XMAX) GO TO 115 SUM = ONE + SUM * (XMAX / SX(I))**2 XMAX = ABS(SX(I)) GO TO 200 C 115 SUM = SUM + (SX(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(SX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + SX(J)**2 SNRM2 = 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 SNRM2 = XMAX * SQRT(SUM) 300 CONTINUE RETURN END *DECK SSCAL SUBROUTINE SSCAL (N, SA, SX, INCX) C***BEGIN PROLOGUE SSCAL C***PURPOSE Multiply a vector by a constant. C***CATEGORY D1A6 C***TYPE SINGLE 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 SA single precision scale factor C SX single precision vector with N elements C INCX storage spacing between elements of SX C C --Output-- C SX single precision result (unchanged if N .LE. 0) C C Replace single precision SX by single precision SA*SX. C For I = 0 to N-1, replace SX(IX+I*INCX) with SA * SX(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 SSCAL REAL SA, SX(*) INTEGER I, INCX, IX, M, MP1, N C***FIRST EXECUTABLE STATEMENT SSCAL 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 SX(IX) = SA*SX(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 SX(I) = SA*SX(I) 30 CONTINUE IF (N .LT. 5) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SX(I) = SA*SX(I) SX(I+1) = SA*SX(I+1) SX(I+2) = SA*SX(I+2) SX(I+3) = SA*SX(I+3) SX(I+4) = SA*SX(I+4) 50 CONTINUE RETURN END *DECK ISAMAX INTEGER FUNCTION ISAMAX (N, SX, INCX) C***BEGIN PROLOGUE ISAMAX C***PURPOSE Find the smallest index of that component of a vector C having the maximum magnitude. C***CATEGORY D1A2 C***TYPE SINGLE 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 SX single precision vector with N elements C INCX storage spacing between elements of SX C C --Output-- C ISAMAX smallest index (zero if N .LE. 0) C C Find smallest index of maximum magnitude of single precision SX. C ISAMAX = first I, I = 1 to N, to maximize ABS(SX(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 861211 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 920618 Slight restructuring of code. (RWC, WRB) C***END PROLOGUE ISAMAX REAL SX(*), SMAX, XMAG INTEGER I, INCX, IX, N C***FIRST EXECUTABLE STATEMENT ISAMAX ISAMAX = 0 IF (N .LE. 0) RETURN ISAMAX = 1 IF (N .EQ. 1) RETURN C 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 SMAX = ABS(SX(IX)) IX = IX + INCX DO 10 I = 2,N XMAG = ABS(SX(IX)) IF (XMAG .GT. SMAX) THEN ISAMAX = I SMAX = XMAG ENDIF IX = IX + INCX 10 CONTINUE RETURN C C Code for increments equal to 1. C 20 SMAX = ABS(SX(1)) DO 30 I = 2,N XMAG = ABS(SX(I)) IF (XMAG .GT. SMAX) THEN ISAMAX = I SMAX = XMAG ENDIF 30 CONTINUE RETURN END *DECK XERRWV SUBROUTINE XERRWV (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2) C***BEGIN PROLOGUE XERRWV C***SUBSIDIARY C***PURPOSE Write error message with values. C***CATEGORY R3C C***TYPE SINGLE PRECISION (XERRWV-S, XERRWD-D) C***AUTHOR Hindmarsh, Alan C., (LLNL) C***DESCRIPTION C C Subroutines XERRWV, 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 single precision and are printed C in E21.13 format. C C***ROUTINES CALLED IXSAV C***REVISION HISTORY (YYMMDD) C 791129 DATE WRITTEN C 890413 Heavily revised, with Common eliminated. (ACH, PNB) 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 XERRWV 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 XERRWV.. None C Function routine called by XERRWV.. IXSAV C----------------------------------------------------------------------- C**End C C Declare arguments. C REAL 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 XERRWV 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 =',E21.13) IF (NR .EQ. 2) WRITE (LUNIT, 50) R1,R2 50 FORMAT(6X,'In above, R1 =',E21.13,3X,'R2 =',E21.13) C C Abort the run if LEVEL = 2. C 100 IF (LEVEL .NE. 2) RETURN STOP C----------------------- End of Subroutine XERRWV ---------------------- 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