--- rpl/lapack/lapack/dgejsv.f 2012/08/22 09:48:13 1.9
+++ rpl/lapack/lapack/dgejsv.f 2020/05/21 21:45:56 1.20
@@ -2,18 +2,18 @@
*
* =========== DOCUMENTATION ===========
*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download DGEJSV + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
+*> Download DGEJSV + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
*> [TXT]
-*> \endhtmlonly
+*> \endhtmlonly
*
* Definition:
* ===========
@@ -21,7 +21,7 @@
* SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
* M, N, A, LDA, SVA, U, LDU, V, LDV,
* WORK, LWORK, IWORK, INFO )
-*
+*
* .. Scalar Arguments ..
* IMPLICIT NONE
* INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
@@ -32,7 +32,7 @@
* INTEGER IWORK( * )
* CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
* ..
-*
+*
*
*> \par Purpose:
* =============
@@ -51,6 +51,8 @@
*> the right singular vectors of [A], respectively. The matrices [U] and [V]
*> are computed and stored in the arrays U and V, respectively. The diagonal
*> of [SIGMA] is computed and stored in the array SVA.
+*> DGEJSV can sometimes compute tiny singular values and their singular vectors much
+*> more accurately than other SVD routines, see below under Further Details.
*> \endverbatim
*
* Arguments:
@@ -80,12 +82,12 @@
*> desirable, then this option is advisable. The input matrix A
*> is preprocessed with QR factorization with FULL (row and
*> column) pivoting.
-*> = 'G' Computation as with 'F' with an additional estimate of the
+*> = 'G': Computation as with 'F' with an additional estimate of the
*> condition number of B, where A=D*B. If A has heavily weighted
*> rows, then using this condition number gives too pessimistic
*> error bound.
*> = 'A': Small singular values are the noise and the matrix is treated
-*> as numerically rank defficient. The error in the computed
+*> as numerically rank deficient. The error in the computed
*> singular values is bounded by f(m,n)*epsilon*||A||.
*> The computed SVD A = U * S * V^t restores A up to
*> f(m,n)*epsilon*||A||.
@@ -131,7 +133,7 @@
*> specified range. If A .NE. 0 is scaled so that the largest singular
*> value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
*> the licence to kill columns of A whose norm in c*A is less than
-*> DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
+*> DSQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
*> where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
*> = 'N': Do not kill small columns of c*A. This option assumes that
*> BLAS and QR factorizations and triangular solvers are
@@ -228,14 +230,14 @@
*> If JOBU = 'F', then U contains on exit the M-by-M matrix of
*> the left singular vectors, including an ONB
*> of the orthogonal complement of the Range(A).
-*> If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
+*> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
*> then U is used as workspace if the procedure
*> replaces A with A^t. In that case, [V] is computed
*> in U as left singular vectors of A^t and then
*> copied back to the V array. This 'W' option is just
*> a reminder to the caller that in this case U is
*> reserved as workspace of length N*N.
-*> If JOBU = 'N' U is not referenced.
+*> If JOBU = 'N' U is not referenced, unless JOBT='T'.
*> \endverbatim
*>
*> \param[in] LDU
@@ -250,14 +252,14 @@
*> V is DOUBLE PRECISION array, dimension ( LDV, N )
*> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
*> the right singular vectors;
-*> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
+*> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
*> then V is used as workspace if the pprocedure
*> replaces A with A^t. In that case, [U] is computed
*> in V as right singular vectors of A^t and then
*> copied back to the U array. This 'W' option is just
*> a reminder to the caller that in this case V is
*> reserved as workspace of length N*N.
-*> If JOBV = 'N' V is not referenced.
+*> If JOBV = 'N' V is not referenced, unless JOBT='T'.
*> \endverbatim
*>
*> \param[in] LDV
@@ -269,14 +271,14 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is DOUBLE PRECISION array, dimension at least LWORK.
-*> On exit, if N.GT.0 .AND. M.GT.0 (else not referenced),
+*> WORK is DOUBLE PRECISION array, dimension (LWORK)
+*> On exit, if N > 0 .AND. M > 0 (else not referenced),
*> WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
*> that SCALE*SVA(1:N) are the computed singular values
*> of A. (See the description of SVA().)
*> WORK(2) = See the description of WORK(1).
*> WORK(3) = SCONDA is an estimate for the condition number of
-*> column equilibrated A. (If JOBA .EQ. 'E' or 'G')
+*> column equilibrated A. (If JOBA = 'E' or 'G')
*> SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
*> It is computed using DPOCON. It holds
*> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
@@ -295,7 +297,7 @@
*> triangular factor in the first QR factorization.
*> WORK(5) = an estimate of the scaled condition number of the
*> triangular factor in the second QR factorization.
-*> The following two parameters are computed if JOBT .EQ. 'T'.
+*> The following two parameters are computed if JOBT = 'T'.
*> They are provided for a developer/implementer who is familiar
*> with the details of the method.
*>
@@ -311,47 +313,47 @@
*> Length of WORK to confirm proper allocation of work space.
*> LWORK depends on the job:
*>
-*> If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
-*> -> .. no scaled condition estimate required (JOBE.EQ.'N'):
+*> If only SIGMA is needed (JOBU = 'N', JOBV = 'N') and
+*> -> .. no scaled condition estimate required (JOBE = 'N'):
*> LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
*> ->> For optimal performance (blocked code) the optimal value
*> is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
*> block size for DGEQP3 and DGEQRF.
-*> In general, optimal LWORK is computed as
-*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7).
+*> In general, optimal LWORK is computed as
+*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7).
*> -> .. an estimate of the scaled condition number of A is
*> required (JOBA='E', 'G'). In this case, LWORK is the maximum
*> of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7).
-*> ->> For optimal performance (blocked code) the optimal value
+*> ->> For optimal performance (blocked code) the optimal value
*> is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7).
*> In general, the optimal length LWORK is computed as
-*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
+*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
*> N+N*N+LWORK(DPOCON),7).
*>
-*> If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
+*> If SIGMA and the right singular vectors are needed (JOBV = 'V'),
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
*> -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
-*> where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ,
+*> where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF,
*> DORMLQ. In general, the optimal length LWORK is computed as
-*> LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON),
-*> N+LWORK(DGELQ), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
+*> LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON),
+*> N+LWORK(DGELQF), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
*>
*> If SIGMA and the left singular vectors are needed
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
*> -> For optimal performance:
-*> if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
-*> if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
+*> if JOBU = 'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
+*> if JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
*> where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
*> In general, the optimal length LWORK is computed as
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
-*> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
-*> Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or
-*> M*NB (for JOBU.EQ.'F').
-*>
-*> If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and
-*> -> if JOBV.EQ.'V'
-*> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
-*> -> if JOBV.EQ.'J' the minimal requirement is
+*> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
+*> Here LWORK(DORMQR) equals N*NB (for JOBU = 'U') or
+*> M*NB (for JOBU = 'F').
+*>
+*> If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
+*> -> if JOBV = 'V'
+*> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
+*> -> if JOBV = 'J' the minimal requirement is
*> LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
*> -> For optimal performance, LWORK should be additionally
*> larger than N+M*NB, where NB is the optimal block size
@@ -360,14 +362,14 @@
*>
*> \param[out] IWORK
*> \verbatim
-*> IWORK is INTEGER array, dimension M+3*N.
+*> IWORK is INTEGER array, dimension (M+3*N).
*> On exit,
*> IWORK(1) = the numerical rank determined after the initial
*> QR factorization with pivoting. See the descriptions
*> of JOBA and JOBR.
*> IWORK(2) = the number of the computed nonzero singular values
*> IWORK(3) = if nonzero, a warning message:
-*> If IWORK(3).EQ.1 then some of the column norms of A
+*> If IWORK(3) = 1 then some of the column norms of A
*> were denormalized floats. The requested high accuracy
*> is not warranted by the data.
*> \endverbatim
@@ -375,23 +377,23 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value.
-*> = 0 : successfull exit;
-*> > 0 : DGEJSV did not converge in the maximal allowed number
-*> of sweeps. The computed values may be inaccurate.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value.
+*> = 0: successful exit;
+*> > 0: DGEJSV did not converge in the maximal allowed number
+*> of sweeps. The computed values may be inaccurate.
*> \endverbatim
*
* Authors:
* ========
*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
*
-*> \date November 2011
+*> \date June 2016
*
-*> \ingroup doubleGEcomputational
+*> \ingroup doubleGEsing
*
*> \par Further Details:
* =====================
@@ -426,8 +428,8 @@
*> The rank revealing QR factorization (in this code: DGEQP3) should be
*> implemented as in [3]. We have a new version of DGEQP3 under development
*> that is more robust than the current one in LAPACK, with a cleaner cut in
-*> rank defficient cases. It will be available in the SIGMA library [4].
-*> If M is much larger than N, it is obvious that the inital QRF with
+*> rank deficient cases. It will be available in the SIGMA library [4].
+*> If M is much larger than N, it is obvious that the initial QRF with
*> column pivoting can be preprocessed by the QRF without pivoting. That
*> well known trick is not used in DGEJSV because in some cases heavy row
*> weighting can be treated with complete pivoting. The overhead in cases
@@ -474,10 +476,10 @@
$ M, N, A, LDA, SVA, U, LDU, V, LDV,
$ WORK, LWORK, IWORK, INFO )
*
-* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK computational routine (version 3.7.1) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2011
+* June 2016
*
* .. Scalar Arguments ..
IMPLICIT NONE
@@ -506,8 +508,7 @@
$ NOSCAL, ROWPIV, RSVEC, TRANSP
* ..
* .. Intrinsic Functions ..
- INTRINSIC DABS, DLOG, DMAX1, DMIN1, DBLE,
- $ MAX0, MIN0, IDNINT, DSIGN, DSQRT
+ INTRINSIC DABS, DLOG, MAX, MIN, DBLE, IDNINT, DSIGN, DSQRT
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DNRM2
@@ -561,19 +562,19 @@
ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
INFO = - 13
ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
- INFO = - 14
+ INFO = - 15
ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND.
- & (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR.
+ & (LWORK .LT. MAX(7,4*N+1,2*M+N))) .OR.
& (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND.
- & (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR.
- & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1)))
+ & (LWORK .LT. MAX(7,4*N+N*N,2*M+N))) .OR.
+ & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1)))
& .OR.
- & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1)))
+ & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1)))
& .OR.
- & (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND.
- & (LWORK.LT.MAX0(2*M+N,6*N+2*N*N)))
+ & (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND.
+ & (LWORK.LT.MAX(2*M+N,6*N+2*N*N)))
& .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND.
- & LWORK.LT.MAX0(2*M+N,4*N+N*N,2*N+N*N+6)))
+ & LWORK.LT.MAX(2*M+N,4*N+N*N,2*N+N*N+6)))
& THEN
INFO = - 17
ELSE
@@ -589,7 +590,11 @@
*
* Quick return for void matrix (Y3K safe)
* #:)
- IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
+ IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
+ IWORK(1:3) = 0
+ WORK(1:7) = 0
+ RETURN
+ ENDIF
*
* Determine whether the matrix U should be M x N or M x M
*
@@ -644,8 +649,8 @@
AAPP = ZERO
AAQQ = BIG
DO 4781 p = 1, N
- AAPP = DMAX1( AAPP, SVA(p) )
- IF ( SVA(p) .NE. ZERO ) AAQQ = DMIN1( AAQQ, SVA(p) )
+ AAPP = MAX( AAPP, SVA(p) )
+ IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) )
4781 CONTINUE
*
* Quick return for zero M x N matrix
@@ -715,6 +720,7 @@
IWORK(1) = 0
IWORK(2) = 0
END IF
+ IWORK(3) = 0
IF ( ERREST ) WORK(3) = ONE
IF ( LSVEC .AND. RSVEC ) THEN
WORK(4) = ONE
@@ -749,14 +755,14 @@
* in one pass through the vector
WORK(M+N+p) = XSC * SCALEM
WORK(N+p) = XSC * (SCALEM*DSQRT(TEMP1))
- AATMAX = DMAX1( AATMAX, WORK(N+p) )
- IF (WORK(N+p) .NE. ZERO) AATMIN = DMIN1(AATMIN,WORK(N+p))
+ AATMAX = MAX( AATMAX, WORK(N+p) )
+ IF (WORK(N+p) .NE. ZERO) AATMIN = MIN(AATMIN,WORK(N+p))
1950 CONTINUE
ELSE
DO 1904 p = 1, M
WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) )
- AATMAX = DMAX1( AATMAX, WORK(M+N+p) )
- AATMIN = DMIN1( AATMIN, WORK(M+N+p) )
+ AATMAX = MAX( AATMAX, WORK(M+N+p) )
+ AATMIN = MIN( AATMIN, WORK(M+N+p) )
1904 CONTINUE
END IF
*
@@ -947,7 +953,7 @@
IF ( L2ABER ) THEN
* Standard absolute error bound suffices. All sigma_i with
* sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
-* agressive enforcement of lower numerical rank by introducing a
+* aggressive enforcement of lower numerical rank by introducing a
* backward error of the order of N*EPSLN*||A||.
TEMP1 = DSQRT(DBLE(N))*EPSLN
DO 3001 p = 2, N
@@ -959,9 +965,9 @@
3001 CONTINUE
3002 CONTINUE
ELSE IF ( L2RANK ) THEN
-* .. similarly as above, only slightly more gentle (less agressive).
+* .. similarly as above, only slightly more gentle (less aggressive).
* Sudden drop on the diagonal of R1 is used as the criterion for
-* close-to-rank-defficient.
+* close-to-rank-deficient.
TEMP1 = DSQRT(SFMIN)
DO 3401 p = 2, N
IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR.
@@ -994,7 +1000,7 @@
MAXPRJ = ONE
DO 3051 p = 2, N
TEMP1 = DABS(A(p,p)) / SVA(IWORK(p))
- MAXPRJ = DMIN1( MAXPRJ, TEMP1 )
+ MAXPRJ = MIN( MAXPRJ, TEMP1 )
3051 CONTINUE
IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
END IF
@@ -1052,7 +1058,7 @@
* Singular Values only
*
* .. transpose A(1:NR,1:N)
- DO 1946 p = 1, MIN0( N-1, NR )
+ DO 1946 p = 1, MIN( N-1, NR )
CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
1946 CONTINUE
*
@@ -1288,7 +1294,7 @@
CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
$ WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
CONDR1 = ONE / DSQRT(TEMP1)
-* .. here need a second oppinion on the condition number
+* .. here need a second opinion on the condition number
* .. then assume worst case scenario
* R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
* more conservative <=> CONDR1 .LT. DSQRT(DBLE(N))
@@ -1308,7 +1314,7 @@
XSC = DSQRT(SMALL)/EPSLN
DO 3959 p = 2, NR
DO 3958 q = 1, p - 1
- TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
+ TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
IF ( DABS(V(q,p)) .LE. TEMP1 )
$ V(q,p) = DSIGN( TEMP1, V(q,p) )
3958 CONTINUE
@@ -1329,7 +1335,7 @@
ELSE
*
* .. ill-conditioned case: second QRF with pivoting
-* Note that windowed pivoting would be equaly good
+* Note that windowed pivoting would be equally good
* numerically, and more run-time efficient. So, in
* an optimal implementation, the next call to DGEQP3
* should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
@@ -1347,7 +1353,7 @@
XSC = DSQRT(SMALL)
DO 3969 p = 2, NR
DO 3968 q = 1, p - 1
- TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
+ TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
IF ( DABS(V(q,p)) .LE. TEMP1 )
$ V(q,p) = DSIGN( TEMP1, V(q,p) )
3968 CONTINUE
@@ -1360,7 +1366,7 @@
XSC = DSQRT(SMALL)
DO 8970 p = 2, NR
DO 8971 q = 1, p - 1
- TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
+ TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
V(p,q) = - DSIGN( TEMP1, V(q,p) )
8971 CONTINUE
8970 CONTINUE
@@ -1382,7 +1388,7 @@
*
IF ( CONDR2 .GE. COND_OK ) THEN
* .. save the Householder vectors used for Q3
-* (this overwrittes the copy of R2, as it will not be
+* (this overwrites the copy of R2, as it will not be
* needed in this branch, but it does not overwritte the
* Huseholder vectors of Q2.).
CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
@@ -1632,7 +1638,7 @@
*
* This branch deploys a preconditioned Jacobi SVD with explicitly
* accumulated rotations. It is included as optional, mainly for
-* experimental purposes. It does perfom well, and can also be used.
+* experimental purposes. It does perform well, and can also be used.
* In this implementation, this branch will be automatically activated
* if the condition number sigma_max(A) / sigma_min(A) is predicted
* to be greater than the overflow threshold. This is because the
@@ -1671,7 +1677,7 @@
XSC = DSQRT(SMALL/EPSLN)
DO 9970 q = 2, NR
DO 9971 p = 1, q - 1
- TEMP1 = XSC * DMIN1(DABS(U(p,p)),DABS(U(q,q)))
+ TEMP1 = XSC * MIN(DABS(U(p,p)),DABS(U(q,q)))
U(p,q) = - DSIGN( TEMP1, U(q,p) )
9971 CONTINUE
9970 CONTINUE