diff --git a/ext/f2c_lapack/Makefile.in b/ext/f2c_lapack/Makefile.in index 7dd08bca1..bb6764363 100755 --- a/ext/f2c_lapack/Makefile.in +++ b/ext/f2c_lapack/Makefile.in @@ -64,6 +64,7 @@ dgetri.o \ dgetrs.o \ dlabad.o \ dlabrd.o \ +dlacon.o \ dlacpy.o \ dlamch.o \ dlange.o \ @@ -88,6 +89,7 @@ dlasrt.o \ dlassq.o \ dlasv2.o \ dlaswp.o \ +dlatrs.o \ dorg2r.o \ dorgbr.o \ dorgl2.o \ @@ -99,6 +101,7 @@ dorml2.o \ dormlq.o \ dormqr.o \ drscl.o \ +dtrcon.o \ dtrtri.o \ dtrti2.o \ dtrtrs.o \ diff --git a/ext/f2c_lapack/dlacon.c b/ext/f2c_lapack/dlacon.c new file mode 100644 index 000000000..2bc7054ab --- /dev/null +++ b/ext/f2c_lapack/dlacon.c @@ -0,0 +1,257 @@ +/* dlacon.f -- translated by f2c (version 20031025). + You must link the resulting object file with libf2c: + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip +*/ + +#include "f2c.h" + +/* Table of constant values */ + +static integer c__1 = 1; +static doublereal c_b11 = 1.; + +/* Subroutine */ int dlacon_(integer *n, doublereal *v, doublereal *x, + integer *isgn, doublereal *est, integer *kase) +{ + /* System generated locals */ + integer i__1; + doublereal d__1; + + /* Builtin functions */ + double d_sign(doublereal *, doublereal *); + integer i_dnnt(doublereal *); + + /* Local variables */ + static integer i__, j, iter; + static doublereal temp; + static integer jump; + extern doublereal dasum_(integer *, doublereal *, integer *); + static integer jlast; + extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, + doublereal *, integer *); + extern integer idamax_(integer *, doublereal *, integer *); + static doublereal altsgn, estold; + + +/* -- LAPACK auxiliary routine (version 3.0) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ +/* Courant Institute, Argonne National Lab, and Rice University */ +/* February 29, 1992 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DLACON estimates the 1-norm of a square, real matrix A. */ +/* Reverse communication is used for evaluating matrix-vector products. */ + +/* Arguments */ +/* ========= */ + +/* N (input) INTEGER */ +/* The order of the matrix. N >= 1. */ + +/* V (workspace) DOUBLE PRECISION array, dimension (N) */ +/* On the final return, V = A*W, where EST = norm(V)/norm(W) */ +/* (W is not returned). */ + +/* X (input/output) DOUBLE PRECISION array, dimension (N) */ +/* On an intermediate return, X should be overwritten by */ +/* A * X, if KASE=1, */ +/* A' * X, if KASE=2, */ +/* and DLACON must be re-called with all the other parameters */ +/* unchanged. */ + +/* ISGN (workspace) INTEGER array, dimension (N) */ + +/* EST (output) DOUBLE PRECISION */ +/* An estimate (a lower bound) for norm(A). */ + +/* KASE (input/output) INTEGER */ +/* On the initial call to DLACON, KASE should be 0. */ +/* On an intermediate return, KASE will be 1 or 2, indicating */ +/* whether X should be overwritten by A * X or A' * X. */ +/* On the final return from DLACON, KASE will again be 0. */ + +/* Further Details */ +/* ======= ======= */ + +/* Contributed by Nick Higham, University of Manchester. */ +/* Originally named SONEST, dated March 16, 1988. */ + +/* Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of */ +/* a real or complex matrix, with applications to condition estimation", */ +/* ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Save statement .. */ +/* .. */ +/* .. Executable Statements .. */ + + /* Parameter adjustments */ + --isgn; + --x; + --v; + + /* Function Body */ + if (*kase == 0) { + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + x[i__] = 1. / (doublereal) (*n); +/* L10: */ + } + *kase = 1; + jump = 1; + return 0; + } + + switch (jump) { + case 1: goto L20; + case 2: goto L40; + case 3: goto L70; + case 4: goto L110; + case 5: goto L140; + } + +/* ................ ENTRY (JUMP = 1) */ +/* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */ + +L20: + if (*n == 1) { + v[1] = x[1]; + *est = abs(v[1]); +/* ... QUIT */ + goto L150; + } + *est = dasum_(n, &x[1], &c__1); + + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + x[i__] = d_sign(&c_b11, &x[i__]); + isgn[i__] = i_dnnt(&x[i__]); +/* L30: */ + } + *kase = 2; + jump = 2; + return 0; + +/* ................ ENTRY (JUMP = 2) */ +/* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */ + +L40: + j = idamax_(n, &x[1], &c__1); + iter = 2; + +/* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */ + +L50: + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + x[i__] = 0.; +/* L60: */ + } + x[j] = 1.; + *kase = 1; + jump = 3; + return 0; + +/* ................ ENTRY (JUMP = 3) */ +/* X HAS BEEN OVERWRITTEN BY A*X. */ + +L70: + dcopy_(n, &x[1], &c__1, &v[1], &c__1); + estold = *est; + *est = dasum_(n, &v[1], &c__1); + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + d__1 = d_sign(&c_b11, &x[i__]); + if (i_dnnt(&d__1) != isgn[i__]) { + goto L90; + } +/* L80: */ + } +/* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */ + goto L120; + +L90: +/* TEST FOR CYCLING. */ + if (*est <= estold) { + goto L120; + } + + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + x[i__] = d_sign(&c_b11, &x[i__]); + isgn[i__] = i_dnnt(&x[i__]); +/* L100: */ + } + *kase = 2; + jump = 4; + return 0; + +/* ................ ENTRY (JUMP = 4) */ +/* X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */ + +L110: + jlast = j; + j = idamax_(n, &x[1], &c__1); + if (x[jlast] != (d__1 = x[j], abs(d__1)) && iter < 5) { + ++iter; + goto L50; + } + +/* ITERATION COMPLETE. FINAL STAGE. */ + +L120: + altsgn = 1.; + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + x[i__] = altsgn * ((doublereal) (i__ - 1) / (doublereal) (*n - 1) + + 1.); + altsgn = -altsgn; +/* L130: */ + } + *kase = 1; + jump = 5; + return 0; + +/* ................ ENTRY (JUMP = 5) */ +/* X HAS BEEN OVERWRITTEN BY A*X. */ + +L140: + temp = dasum_(n, &x[1], &c__1) / (doublereal) (*n * 3) * 2.; + if (temp > *est) { + dcopy_(n, &x[1], &c__1, &v[1], &c__1); + *est = temp; + } + +L150: + *kase = 0; + return 0; + +/* End of DLACON */ + +} /* dlacon_ */ + diff --git a/ext/f2c_lapack/dlatrs.c b/ext/f2c_lapack/dlatrs.c new file mode 100644 index 000000000..a08a0f09f --- /dev/null +++ b/ext/f2c_lapack/dlatrs.c @@ -0,0 +1,820 @@ +/* dlatrs.f -- translated by f2c (version 20031025). + You must link the resulting object file with libf2c: + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip +*/ + +#include "f2c.h" + +/* Table of constant values */ + +static integer c__1 = 1; +static doublereal c_b36 = .5; + +/* Subroutine */ int dlatrs_(char *uplo, char *trans, char *diag, char * + normin, integer *n, doublereal *a, integer *lda, doublereal *x, + doublereal *scale, doublereal *cnorm, integer *info, ftnlen uplo_len, + ftnlen trans_len, ftnlen diag_len, ftnlen normin_len) +{ + /* System generated locals */ + integer a_dim1, a_offset, i__1, i__2, i__3; + doublereal d__1, d__2, d__3; + + /* Local variables */ + static integer i__, j; + static doublereal xj, rec, tjj; + static integer jinc; + extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, + integer *); + static doublereal xbnd; + static integer imax; + static doublereal tmax, tjjs, xmax, grow, sumj; + extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, + integer *); + extern logical lsame_(char *, char *, ftnlen, ftnlen); + static doublereal tscal, uscal; + extern doublereal dasum_(integer *, doublereal *, integer *); + static integer jlast; + extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, + integer *, doublereal *, integer *); + static logical upper; + extern /* Subroutine */ int dtrsv_(char *, char *, char *, integer *, + doublereal *, integer *, doublereal *, integer *, ftnlen, ftnlen, + ftnlen); + extern doublereal dlamch_(char *, ftnlen); + extern integer idamax_(integer *, doublereal *, integer *); + extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen); + static doublereal bignum; + static logical notran; + static integer jfirst; + static doublereal smlnum; + static logical nounit; + + +/* -- LAPACK auxiliary routine (version 3.0) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ +/* Courant Institute, Argonne National Lab, and Rice University */ +/* June 30, 1992 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DLATRS solves one of the triangular systems */ + +/* A *x = s*b or A'*x = s*b */ + +/* with scaling to prevent overflow. Here A is an upper or lower */ +/* triangular matrix, A' denotes the transpose of A, x and b are */ +/* n-element vectors, and s is a scaling factor, usually less than */ +/* or equal to 1, chosen so that the components of x will be less than */ +/* the overflow threshold. If the unscaled problem will not cause */ +/* overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A */ +/* is singular (A(j,j) = 0 for some j), then s is set to 0 and a */ +/* non-trivial solution to A*x = 0 is returned. */ + +/* Arguments */ +/* ========= */ + +/* UPLO (input) CHARACTER*1 */ +/* Specifies whether the matrix A is upper or lower triangular. */ +/* = 'U': Upper triangular */ +/* = 'L': Lower triangular */ + +/* TRANS (input) CHARACTER*1 */ +/* Specifies the operation applied to A. */ +/* = 'N': Solve A * x = s*b (No transpose) */ +/* = 'T': Solve A'* x = s*b (Transpose) */ +/* = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose) */ + +/* DIAG (input) CHARACTER*1 */ +/* Specifies whether or not the matrix A is unit triangular. */ +/* = 'N': Non-unit triangular */ +/* = 'U': Unit triangular */ + +/* NORMIN (input) CHARACTER*1 */ +/* Specifies whether CNORM has been set or not. */ +/* = 'Y': CNORM contains the column norms on entry */ +/* = 'N': CNORM is not set on entry. On exit, the norms will */ +/* be computed and stored in CNORM. */ + +/* N (input) INTEGER */ +/* The order of the matrix A. N >= 0. */ + +/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */ +/* The triangular matrix A. If UPLO = 'U', the leading n by n */ +/* upper triangular part of the array A contains the upper */ +/* triangular matrix, and the strictly lower triangular part of */ +/* A is not referenced. If UPLO = 'L', the leading n by n lower */ +/* triangular part of the array A contains the lower triangular */ +/* matrix, and the strictly upper triangular part of A is not */ +/* referenced. If DIAG = 'U', the diagonal elements of A are */ +/* also not referenced and are assumed to be 1. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max (1,N). */ + +/* X (input/output) DOUBLE PRECISION array, dimension (N) */ +/* On entry, the right hand side b of the triangular system. */ +/* On exit, X is overwritten by the solution vector x. */ + +/* SCALE (output) DOUBLE PRECISION */ +/* The scaling factor s for the triangular system */ +/* A * x = s*b or A'* x = s*b. */ +/* If SCALE = 0, the matrix A is singular or badly scaled, and */ +/* the vector x is an exact or approximate solution to A*x = 0. */ + +/* CNORM (input or output) DOUBLE PRECISION array, dimension (N) */ + +/* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) */ +/* contains the norm of the off-diagonal part of the j-th column */ +/* of A. If TRANS = 'N', CNORM(j) must be greater than or equal */ +/* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) */ +/* must be greater than or equal to the 1-norm. */ + +/* If NORMIN = 'N', CNORM is an output argument and CNORM(j) */ +/* returns the 1-norm of the offdiagonal part of the j-th column */ +/* of A. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -k, the k-th argument had an illegal value */ + +/* Further Details */ +/* ======= ======= */ + +/* A rough bound on x is computed; if that is less than overflow, DTRSV */ +/* is called, otherwise, specific code is used which checks for possible */ +/* overflow or divide-by-zero at every operation. */ + +/* A columnwise scheme is used for solving A*x = b. The basic algorithm */ +/* if A is lower triangular is */ + +/* x[1:n] := b[1:n] */ +/* for j = 1, ..., n */ +/* x(j) := x(j) / A(j,j) */ +/* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] */ +/* end */ + +/* Define bounds on the components of x after j iterations of the loop: */ +/* M(j) = bound on x[1:j] */ +/* G(j) = bound on x[j+1:n] */ +/* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. */ + +/* Then for iteration j+1 we have */ +/* M(j+1) <= G(j) / | A(j+1,j+1) | */ +/* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | */ +/* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) */ + +/* where CNORM(j+1) is greater than or equal to the infinity-norm of */ +/* column j+1 of A, not counting the diagonal. Hence */ + +/* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) */ +/* 1<=i<=j */ +/* and */ + +/* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) */ +/* 1<=i< j */ + +/* Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the */ +/* reciprocal of the largest M(j), j=1,..,n, is larger than */ +/* max(underflow, 1/overflow). */ + +/* The bound on x(j) is also used to determine when a step in the */ +/* columnwise method can be performed without fear of overflow. If */ +/* the computed bound is greater than a large constant, x is scaled to */ +/* prevent overflow, but if the bound overflows, x is set to 0, x(j) to */ +/* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. */ + +/* Similarly, a row-wise scheme is used to solve A'*x = b. The basic */ +/* algorithm for A upper triangular is */ + +/* for j = 1, ..., n */ +/* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) */ +/* end */ + +/* We simultaneously compute two bounds */ +/* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j */ +/* M(j) = bound on x(i), 1<=i<=j */ + +/* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we */ +/* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. */ +/* Then the bound on x(j) is */ + +/* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | */ + +/* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) */ +/* 1<=i<=j */ + +/* and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater */ +/* than max(underflow, 1/overflow). */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Executable Statements .. */ + + /* Parameter adjustments */ + a_dim1 = *lda; + a_offset = 1 + a_dim1; + a -= a_offset; + --x; + --cnorm; + + /* Function Body */ + *info = 0; + upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1); + notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1); + nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1); + +/* Test the input parameters. */ + + if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) { + *info = -1; + } else if (! notran && ! lsame_(trans, "T", (ftnlen)1, (ftnlen)1) && ! + lsame_(trans, "C", (ftnlen)1, (ftnlen)1)) { + *info = -2; + } else if (! nounit && ! lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) { + *info = -3; + } else if (! lsame_(normin, "Y", (ftnlen)1, (ftnlen)1) && ! lsame_(normin, + "N", (ftnlen)1, (ftnlen)1)) { + *info = -4; + } else if (*n < 0) { + *info = -5; + } else if (*lda < max(1,*n)) { + *info = -7; + } + if (*info != 0) { + i__1 = -(*info); + xerbla_("DLATRS", &i__1, (ftnlen)6); + return 0; + } + +/* Quick return if possible */ + + if (*n == 0) { + return 0; + } + +/* Determine machine dependent parameters to control overflow. */ + + smlnum = dlamch_("Safe minimum", (ftnlen)12) / dlamch_("Precision", ( + ftnlen)9); + bignum = 1. / smlnum; + *scale = 1.; + + if (lsame_(normin, "N", (ftnlen)1, (ftnlen)1)) { + +/* Compute the 1-norm of each column, not including the diagonal. */ + + if (upper) { + +/* A is upper triangular. */ + + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = j - 1; + cnorm[j] = dasum_(&i__2, &a[j * a_dim1 + 1], &c__1); +/* L10: */ + } + } else { + +/* A is lower triangular. */ + + i__1 = *n - 1; + for (j = 1; j <= i__1; ++j) { + i__2 = *n - j; + cnorm[j] = dasum_(&i__2, &a[j + 1 + j * a_dim1], &c__1); +/* L20: */ + } + cnorm[*n] = 0.; + } + } + +/* Scale the column norms by TSCAL if the maximum element in CNORM is */ +/* greater than BIGNUM. */ + + imax = idamax_(n, &cnorm[1], &c__1); + tmax = cnorm[imax]; + if (tmax <= bignum) { + tscal = 1.; + } else { + tscal = 1. / (smlnum * tmax); + dscal_(n, &tscal, &cnorm[1], &c__1); + } + +/* Compute a bound on the computed solution vector to see if the */ +/* Level 2 BLAS routine DTRSV can be used. */ + + j = idamax_(n, &x[1], &c__1); + xmax = (d__1 = x[j], abs(d__1)); + xbnd = xmax; + if (notran) { + +/* Compute the growth in A * x = b. */ + + if (upper) { + jfirst = *n; + jlast = 1; + jinc = -1; + } else { + jfirst = 1; + jlast = *n; + jinc = 1; + } + + if (tscal != 1.) { + grow = 0.; + goto L50; + } + + if (nounit) { + +/* A is non-unit triangular. */ + +/* Compute GROW = 1/G(j) and XBND = 1/M(j). */ +/* Initially, G(0) = max{x(i), i=1,...,n}. */ + + grow = 1. / max(xbnd,smlnum); + xbnd = grow; + i__1 = jlast; + i__2 = jinc; + for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { + +/* Exit the loop if the growth factor is too small. */ + + if (grow <= smlnum) { + goto L50; + } + +/* M(j) = G(j-1) / abs(A(j,j)) */ + + tjj = (d__1 = a[j + j * a_dim1], abs(d__1)); +/* Computing MIN */ + d__1 = xbnd, d__2 = min(1.,tjj) * grow; + xbnd = min(d__1,d__2); + if (tjj + cnorm[j] >= smlnum) { + +/* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) */ + + grow *= tjj / (tjj + cnorm[j]); + } else { + +/* G(j) could overflow, set GROW to 0. */ + + grow = 0.; + } +/* L30: */ + } + grow = xbnd; + } else { + +/* A is unit triangular. */ + +/* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */ + +/* Computing MIN */ + d__1 = 1., d__2 = 1. / max(xbnd,smlnum); + grow = min(d__1,d__2); + i__2 = jlast; + i__1 = jinc; + for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { + +/* Exit the loop if the growth factor is too small. */ + + if (grow <= smlnum) { + goto L50; + } + +/* G(j) = G(j-1)*( 1 + CNORM(j) ) */ + + grow *= 1. / (cnorm[j] + 1.); +/* L40: */ + } + } +L50: + + ; + } else { + +/* Compute the growth in A' * x = b. */ + + if (upper) { + jfirst = 1; + jlast = *n; + jinc = 1; + } else { + jfirst = *n; + jlast = 1; + jinc = -1; + } + + if (tscal != 1.) { + grow = 0.; + goto L80; + } + + if (nounit) { + +/* A is non-unit triangular. */ + +/* Compute GROW = 1/G(j) and XBND = 1/M(j). */ +/* Initially, M(0) = max{x(i), i=1,...,n}. */ + + grow = 1. / max(xbnd,smlnum); + xbnd = grow; + i__1 = jlast; + i__2 = jinc; + for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { + +/* Exit the loop if the growth factor is too small. */ + + if (grow <= smlnum) { + goto L80; + } + +/* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) */ + + xj = cnorm[j] + 1.; +/* Computing MIN */ + d__1 = grow, d__2 = xbnd / xj; + grow = min(d__1,d__2); + +/* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) */ + + tjj = (d__1 = a[j + j * a_dim1], abs(d__1)); + if (xj > tjj) { + xbnd *= tjj / xj; + } +/* L60: */ + } + grow = min(grow,xbnd); + } else { + +/* A is unit triangular. */ + +/* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. */ + +/* Computing MIN */ + d__1 = 1., d__2 = 1. / max(xbnd,smlnum); + grow = min(d__1,d__2); + i__2 = jlast; + i__1 = jinc; + for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { + +/* Exit the loop if the growth factor is too small. */ + + if (grow <= smlnum) { + goto L80; + } + +/* G(j) = ( 1 + CNORM(j) )*G(j-1) */ + + xj = cnorm[j] + 1.; + grow /= xj; +/* L70: */ + } + } +L80: + ; + } + + if (grow * tscal > smlnum) { + +/* Use the Level 2 BLAS solve if the reciprocal of the bound on */ +/* elements of X is not too small. */ + + dtrsv_(uplo, trans, diag, n, &a[a_offset], lda, &x[1], &c__1, (ftnlen) + 1, (ftnlen)1, (ftnlen)1); + } else { + +/* Use a Level 1 BLAS solve, scaling intermediate results. */ + + if (xmax > bignum) { + +/* Scale X so that its components are less than or equal to */ +/* BIGNUM in absolute value. */ + + *scale = bignum / xmax; + dscal_(n, scale, &x[1], &c__1); + xmax = bignum; + } + + if (notran) { + +/* Solve A * x = b */ + + i__1 = jlast; + i__2 = jinc; + for (j = jfirst; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { + +/* Compute x(j) = b(j) / A(j,j), scaling x if necessary. */ + + xj = (d__1 = x[j], abs(d__1)); + if (nounit) { + tjjs = a[j + j * a_dim1] * tscal; + } else { + tjjs = tscal; + if (tscal == 1.) { + goto L100; + } + } + tjj = abs(tjjs); + if (tjj > smlnum) { + +/* abs(A(j,j)) > SMLNUM: */ + + if (tjj < 1.) { + if (xj > tjj * bignum) { + +/* Scale x by 1/b(j). */ + + rec = 1. / xj; + dscal_(n, &rec, &x[1], &c__1); + *scale *= rec; + xmax *= rec; + } + } + x[j] /= tjjs; + xj = (d__1 = x[j], abs(d__1)); + } else if (tjj > 0.) { + +/* 0 < abs(A(j,j)) <= SMLNUM: */ + + if (xj > tjj * bignum) { + +/* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM */ +/* to avoid overflow when dividing by A(j,j). */ + + rec = tjj * bignum / xj; + if (cnorm[j] > 1.) { + +/* Scale by 1/CNORM(j) to avoid overflow when */ +/* multiplying x(j) times column j. */ + + rec /= cnorm[j]; + } + dscal_(n, &rec, &x[1], &c__1); + *scale *= rec; + xmax *= rec; + } + x[j] /= tjjs; + xj = (d__1 = x[j], abs(d__1)); + } else { + +/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */ +/* scale = 0, and compute a solution to A*x = 0. */ + + i__3 = *n; + for (i__ = 1; i__ <= i__3; ++i__) { + x[i__] = 0.; +/* L90: */ + } + x[j] = 1.; + xj = 1.; + *scale = 0.; + xmax = 0.; + } +L100: + +/* Scale x if necessary to avoid overflow when adding a */ +/* multiple of column j of A. */ + + if (xj > 1.) { + rec = 1. / xj; + if (cnorm[j] > (bignum - xmax) * rec) { + +/* Scale x by 1/(2*abs(x(j))). */ + + rec *= .5; + dscal_(n, &rec, &x[1], &c__1); + *scale *= rec; + } + } else if (xj * cnorm[j] > bignum - xmax) { + +/* Scale x by 1/2. */ + + dscal_(n, &c_b36, &x[1], &c__1); + *scale *= .5; + } + + if (upper) { + if (j > 1) { + +/* Compute the update */ +/* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) */ + + i__3 = j - 1; + d__1 = -x[j] * tscal; + daxpy_(&i__3, &d__1, &a[j * a_dim1 + 1], &c__1, &x[1], + &c__1); + i__3 = j - 1; + i__ = idamax_(&i__3, &x[1], &c__1); + xmax = (d__1 = x[i__], abs(d__1)); + } + } else { + if (j < *n) { + +/* Compute the update */ +/* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) */ + + i__3 = *n - j; + d__1 = -x[j] * tscal; + daxpy_(&i__3, &d__1, &a[j + 1 + j * a_dim1], &c__1, & + x[j + 1], &c__1); + i__3 = *n - j; + i__ = j + idamax_(&i__3, &x[j + 1], &c__1); + xmax = (d__1 = x[i__], abs(d__1)); + } + } +/* L110: */ + } + + } else { + +/* Solve A' * x = b */ + + i__2 = jlast; + i__1 = jinc; + for (j = jfirst; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { + +/* Compute x(j) = b(j) - sum A(k,j)*x(k). */ +/* k<>j */ + + xj = (d__1 = x[j], abs(d__1)); + uscal = tscal; + rec = 1. / max(xmax,1.); + if (cnorm[j] > (bignum - xj) * rec) { + +/* If x(j) could overflow, scale x by 1/(2*XMAX). */ + + rec *= .5; + if (nounit) { + tjjs = a[j + j * a_dim1] * tscal; + } else { + tjjs = tscal; + } + tjj = abs(tjjs); + if (tjj > 1.) { + +/* Divide by A(j,j) when scaling x if A(j,j) > 1. */ + +/* Computing MIN */ + d__1 = 1., d__2 = rec * tjj; + rec = min(d__1,d__2); + uscal /= tjjs; + } + if (rec < 1.) { + dscal_(n, &rec, &x[1], &c__1); + *scale *= rec; + xmax *= rec; + } + } + + sumj = 0.; + if (uscal == 1.) { + +/* If the scaling needed for A in the dot product is 1, */ +/* call DDOT to perform the dot product. */ + + if (upper) { + i__3 = j - 1; + sumj = ddot_(&i__3, &a[j * a_dim1 + 1], &c__1, &x[1], + &c__1); + } else if (j < *n) { + i__3 = *n - j; + sumj = ddot_(&i__3, &a[j + 1 + j * a_dim1], &c__1, &x[ + j + 1], &c__1); + } + } else { + +/* Otherwise, use in-line code for the dot product. */ + + if (upper) { + i__3 = j - 1; + for (i__ = 1; i__ <= i__3; ++i__) { + sumj += a[i__ + j * a_dim1] * uscal * x[i__]; +/* L120: */ + } + } else if (j < *n) { + i__3 = *n; + for (i__ = j + 1; i__ <= i__3; ++i__) { + sumj += a[i__ + j * a_dim1] * uscal * x[i__]; +/* L130: */ + } + } + } + + if (uscal == tscal) { + +/* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) */ +/* was not used to scale the dotproduct. */ + + x[j] -= sumj; + xj = (d__1 = x[j], abs(d__1)); + if (nounit) { + tjjs = a[j + j * a_dim1] * tscal; + } else { + tjjs = tscal; + if (tscal == 1.) { + goto L150; + } + } + +/* Compute x(j) = x(j) / A(j,j), scaling if necessary. */ + + tjj = abs(tjjs); + if (tjj > smlnum) { + +/* abs(A(j,j)) > SMLNUM: */ + + if (tjj < 1.) { + if (xj > tjj * bignum) { + +/* Scale X by 1/abs(x(j)). */ + + rec = 1. / xj; + dscal_(n, &rec, &x[1], &c__1); + *scale *= rec; + xmax *= rec; + } + } + x[j] /= tjjs; + } else if (tjj > 0.) { + +/* 0 < abs(A(j,j)) <= SMLNUM: */ + + if (xj > tjj * bignum) { + +/* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. */ + + rec = tjj * bignum / xj; + dscal_(n, &rec, &x[1], &c__1); + *scale *= rec; + xmax *= rec; + } + x[j] /= tjjs; + } else { + +/* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and */ +/* scale = 0, and compute a solution to A'*x = 0. */ + + i__3 = *n; + for (i__ = 1; i__ <= i__3; ++i__) { + x[i__] = 0.; +/* L140: */ + } + x[j] = 1.; + *scale = 0.; + xmax = 0.; + } +L150: + ; + } else { + +/* Compute x(j) := x(j) / A(j,j) - sumj if the dot */ +/* product has already been divided by 1/A(j,j). */ + + x[j] = x[j] / tjjs - sumj; + } +/* Computing MAX */ + d__2 = xmax, d__3 = (d__1 = x[j], abs(d__1)); + xmax = max(d__2,d__3); +/* L160: */ + } + } + *scale /= tscal; + } + +/* Scale the column norms by 1/TSCAL for return. */ + + if (tscal != 1.) { + d__1 = 1. / tscal; + dscal_(n, &d__1, &cnorm[1], &c__1); + } + + return 0; + +/* End of DLATRS */ + +} /* dlatrs_ */ + diff --git a/ext/f2c_lapack/dtrcon.c b/ext/f2c_lapack/dtrcon.c new file mode 100644 index 000000000..6323a364e --- /dev/null +++ b/ext/f2c_lapack/dtrcon.c @@ -0,0 +1,242 @@ +/* dtrcon.f -- translated by f2c (version 20031025). + You must link the resulting object file with libf2c: + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip +*/ + +#include "f2c.h" + +/* Table of constant values */ + +static integer c__1 = 1; + +/* Subroutine */ int dtrcon_(char *norm, char *uplo, char *diag, integer *n, + doublereal *a, integer *lda, doublereal *rcond, doublereal *work, + integer *iwork, integer *info, ftnlen norm_len, ftnlen uplo_len, + ftnlen diag_len) +{ + /* System generated locals */ + integer a_dim1, a_offset, i__1; + doublereal d__1; + + /* Local variables */ + static integer ix, kase, kase1; + static doublereal scale; + extern logical lsame_(char *, char *, ftnlen, ftnlen); + extern /* Subroutine */ int drscl_(integer *, doublereal *, doublereal *, + integer *); + static doublereal anorm; + static logical upper; + static doublereal xnorm; + extern doublereal dlamch_(char *, ftnlen); + extern /* Subroutine */ int dlacon_(integer *, doublereal *, doublereal *, + integer *, doublereal *, integer *); + extern integer idamax_(integer *, doublereal *, integer *); + extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen); + extern doublereal dlantr_(char *, char *, char *, integer *, integer *, + doublereal *, integer *, doublereal *, ftnlen, ftnlen, ftnlen); + static doublereal ainvnm; + extern /* Subroutine */ int dlatrs_(char *, char *, char *, char *, + integer *, doublereal *, integer *, doublereal *, doublereal *, + doublereal *, integer *, ftnlen, ftnlen, ftnlen, ftnlen); + static logical onenrm; + static char normin[1]; + static doublereal smlnum; + static logical nounit; + + +/* -- LAPACK routine (version 3.0) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ +/* Courant Institute, Argonne National Lab, and Rice University */ +/* March 31, 1993 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DTRCON estimates the reciprocal of the condition number of a */ +/* triangular matrix A, in either the 1-norm or the infinity-norm. */ + +/* The norm of A is computed and an estimate is obtained for */ +/* norm(inv(A)), then the reciprocal of the condition number is */ +/* computed as */ +/* RCOND = 1 / ( norm(A) * norm(inv(A)) ). */ + +/* Arguments */ +/* ========= */ + +/* NORM (input) CHARACTER*1 */ +/* Specifies whether the 1-norm condition number or the */ +/* infinity-norm condition number is required: */ +/* = '1' or 'O': 1-norm; */ +/* = 'I': Infinity-norm. */ + +/* UPLO (input) CHARACTER*1 */ +/* = 'U': A is upper triangular; */ +/* = 'L': A is lower triangular. */ + +/* DIAG (input) CHARACTER*1 */ +/* = 'N': A is non-unit triangular; */ +/* = 'U': A is unit triangular. */ + +/* N (input) INTEGER */ +/* The order of the matrix A. N >= 0. */ + +/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */ +/* The triangular matrix A. If UPLO = 'U', the leading N-by-N */ +/* upper triangular part of the array A contains the upper */ +/* triangular matrix, and the strictly lower triangular part of */ +/* A is not referenced. If UPLO = 'L', the leading N-by-N lower */ +/* triangular part of the array A contains the lower triangular */ +/* matrix, and the strictly upper triangular part of A is not */ +/* referenced. If DIAG = 'U', the diagonal elements of A are */ +/* also not referenced and are assumed to be 1. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,N). */ + +/* RCOND (output) DOUBLE PRECISION */ +/* The reciprocal of the condition number of the matrix A, */ +/* computed as RCOND = 1/(norm(A) * norm(inv(A))). */ + +/* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */ + +/* IWORK (workspace) INTEGER array, dimension (N) */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Executable Statements .. */ + +/* Test the input parameters. */ + + /* Parameter adjustments */ + a_dim1 = *lda; + a_offset = 1 + a_dim1; + a -= a_offset; + --work; + --iwork; + + /* Function Body */ + *info = 0; + upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1); + onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O", (ftnlen)1, ( + ftnlen)1); + nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1); + + if (! onenrm && ! lsame_(norm, "I", (ftnlen)1, (ftnlen)1)) { + *info = -1; + } else if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) { + *info = -2; + } else if (! nounit && ! lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) { + *info = -3; + } else if (*n < 0) { + *info = -4; + } else if (*lda < max(1,*n)) { + *info = -6; + } + if (*info != 0) { + i__1 = -(*info); + xerbla_("DTRCON", &i__1, (ftnlen)6); + return 0; + } + +/* Quick return if possible */ + + if (*n == 0) { + *rcond = 1.; + return 0; + } + + *rcond = 0.; + smlnum = dlamch_("Safe minimum", (ftnlen)12) * (doublereal) max(1,*n); + +/* Compute the norm of the triangular matrix A. */ + + anorm = dlantr_(norm, uplo, diag, n, n, &a[a_offset], lda, &work[1], ( + ftnlen)1, (ftnlen)1, (ftnlen)1); + +/* Continue only if ANORM > 0. */ + + if (anorm > 0.) { + +/* Estimate the norm of the inverse of A. */ + + ainvnm = 0.; + *(unsigned char *)normin = 'N'; + if (onenrm) { + kase1 = 1; + } else { + kase1 = 2; + } + kase = 0; +L10: + dlacon_(n, &work[*n + 1], &work[1], &iwork[1], &ainvnm, &kase); + if (kase != 0) { + if (kase == kase1) { + +/* Multiply by inv(A). */ + + dlatrs_(uplo, "No transpose", diag, normin, n, &a[a_offset], + lda, &work[1], &scale, &work[(*n << 1) + 1], info, ( + ftnlen)1, (ftnlen)12, (ftnlen)1, (ftnlen)1); + } else { + +/* Multiply by inv(A'). */ + + dlatrs_(uplo, "Transpose", diag, normin, n, &a[a_offset], lda, + &work[1], &scale, &work[(*n << 1) + 1], info, ( + ftnlen)1, (ftnlen)9, (ftnlen)1, (ftnlen)1); + } + *(unsigned char *)normin = 'Y'; + +/* Multiply by 1/SCALE if doing so will not cause overflow. */ + + if (scale != 1.) { + ix = idamax_(n, &work[1], &c__1); + xnorm = (d__1 = work[ix], abs(d__1)); + if (scale < xnorm * smlnum || scale == 0.) { + goto L20; + } + drscl_(n, &scale, &work[1], &c__1); + } + goto L10; + } + +/* Compute the estimate of the reciprocal condition number. */ + + if (ainvnm != 0.) { + *rcond = 1. / anorm / ainvnm; + } + } + +L20: + return 0; + +/* End of DTRCON */ + +} /* dtrcon_ */ + diff --git a/ext/lapack/Makefile.in b/ext/lapack/Makefile.in index a0db0ac4c..c936ba49f 100755 --- a/ext/lapack/Makefile.in +++ b/ext/lapack/Makefile.in @@ -35,6 +35,7 @@ dgetri.o \ dgetrs.o \ dlabad.o \ dlabrd.o \ +dlacon.o \ dlacpy.o \ dlamch.o \ dlange.o \ @@ -45,6 +46,7 @@ dlarfb.o \ dlarfg.o \ dlarft.o \ dlartg.o \ +dlatrs.o \ dlas2.o \ dlascl.o \ dlaset.o \ @@ -68,6 +70,7 @@ dorml2.o \ dormlq.o \ dormqr.o \ drscl.o \ +dtrcon.o \ dtrtrs.o \ dgerfs.o \ dgecon.o \ diff --git a/ext/lapack/dlacon.f b/ext/lapack/dlacon.f new file mode 100644 index 000000000..4efa0e816 --- /dev/null +++ b/ext/lapack/dlacon.f @@ -0,0 +1,204 @@ + SUBROUTINE DLACON( N, V, X, ISGN, EST, KASE ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* February 29, 1992 +* +* .. Scalar Arguments .. + INTEGER KASE, N + DOUBLE PRECISION EST +* .. +* .. Array Arguments .. + INTEGER ISGN( * ) + DOUBLE PRECISION V( * ), X( * ) +* .. +* +* Purpose +* ======= +* +* DLACON estimates the 1-norm of a square, real matrix A. +* Reverse communication is used for evaluating matrix-vector products. +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix. N >= 1. +* +* V (workspace) DOUBLE PRECISION array, dimension (N) +* On the final return, V = A*W, where EST = norm(V)/norm(W) +* (W is not returned). +* +* X (input/output) DOUBLE PRECISION array, dimension (N) +* On an intermediate return, X should be overwritten by +* A * X, if KASE=1, +* A' * X, if KASE=2, +* and DLACON must be re-called with all the other parameters +* unchanged. +* +* ISGN (workspace) INTEGER array, dimension (N) +* +* EST (output) DOUBLE PRECISION +* An estimate (a lower bound) for norm(A). +* +* KASE (input/output) INTEGER +* On the initial call to DLACON, KASE should be 0. +* On an intermediate return, KASE will be 1 or 2, indicating +* whether X should be overwritten by A * X or A' * X. +* On the final return from DLACON, KASE will again be 0. +* +* Further Details +* ======= ======= +* +* Contributed by Nick Higham, University of Manchester. +* Originally named SONEST, dated March 16, 1988. +* +* Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of +* a real or complex matrix, with applications to condition estimation", +* ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. +* +* ===================================================================== +* +* .. Parameters .. + INTEGER ITMAX + PARAMETER ( ITMAX = 5 ) + DOUBLE PRECISION ZERO, ONE, TWO + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, ITER, J, JLAST, JUMP + DOUBLE PRECISION ALTSGN, ESTOLD, TEMP +* .. +* .. External Functions .. + INTEGER IDAMAX + DOUBLE PRECISION DASUM + EXTERNAL IDAMAX, DASUM +* .. +* .. External Subroutines .. + EXTERNAL DCOPY +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, NINT, SIGN +* .. +* .. Save statement .. + SAVE +* .. +* .. Executable Statements .. +* + IF( KASE.EQ.0 ) THEN + DO 10 I = 1, N + X( I ) = ONE / DBLE( N ) + 10 CONTINUE + KASE = 1 + JUMP = 1 + RETURN + END IF +* + GO TO ( 20, 40, 70, 110, 140 )JUMP +* +* ................ ENTRY (JUMP = 1) +* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. +* + 20 CONTINUE + IF( N.EQ.1 ) THEN + V( 1 ) = X( 1 ) + EST = ABS( V( 1 ) ) +* ... QUIT + GO TO 150 + END IF + EST = DASUM( N, X, 1 ) +* + DO 30 I = 1, N + X( I ) = SIGN( ONE, X( I ) ) + ISGN( I ) = NINT( X( I ) ) + 30 CONTINUE + KASE = 2 + JUMP = 2 + RETURN +* +* ................ ENTRY (JUMP = 2) +* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. +* + 40 CONTINUE + J = IDAMAX( N, X, 1 ) + ITER = 2 +* +* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. +* + 50 CONTINUE + DO 60 I = 1, N + X( I ) = ZERO + 60 CONTINUE + X( J ) = ONE + KASE = 1 + JUMP = 3 + RETURN +* +* ................ ENTRY (JUMP = 3) +* X HAS BEEN OVERWRITTEN BY A*X. +* + 70 CONTINUE + CALL DCOPY( N, X, 1, V, 1 ) + ESTOLD = EST + EST = DASUM( N, V, 1 ) + DO 80 I = 1, N + IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) ) + $ GO TO 90 + 80 CONTINUE +* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. + GO TO 120 +* + 90 CONTINUE +* TEST FOR CYCLING. + IF( EST.LE.ESTOLD ) + $ GO TO 120 +* + DO 100 I = 1, N + X( I ) = SIGN( ONE, X( I ) ) + ISGN( I ) = NINT( X( I ) ) + 100 CONTINUE + KASE = 2 + JUMP = 4 + RETURN +* +* ................ ENTRY (JUMP = 4) +* X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. +* + 110 CONTINUE + JLAST = J + J = IDAMAX( N, X, 1 ) + IF( ( X( JLAST ).NE.ABS( X( J ) ) ) .AND. ( ITER.LT.ITMAX ) ) THEN + ITER = ITER + 1 + GO TO 50 + END IF +* +* ITERATION COMPLETE. FINAL STAGE. +* + 120 CONTINUE + ALTSGN = ONE + DO 130 I = 1, N + X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) + ALTSGN = -ALTSGN + 130 CONTINUE + KASE = 1 + JUMP = 5 + RETURN +* +* ................ ENTRY (JUMP = 5) +* X HAS BEEN OVERWRITTEN BY A*X. +* + 140 CONTINUE + TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) ) + IF( TEMP.GT.EST ) THEN + CALL DCOPY( N, X, 1, V, 1 ) + EST = TEMP + END IF +* + 150 CONTINUE + KASE = 0 + RETURN +* +* End of DLACON +* + END diff --git a/ext/lapack/dlatrs.f b/ext/lapack/dlatrs.f new file mode 100644 index 000000000..591c966d2 --- /dev/null +++ b/ext/lapack/dlatrs.f @@ -0,0 +1,702 @@ + SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, + $ CNORM, INFO ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* June 30, 1992 +* +* .. Scalar Arguments .. + CHARACTER DIAG, NORMIN, TRANS, UPLO + INTEGER INFO, LDA, N + DOUBLE PRECISION SCALE +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * ) +* .. +* +* Purpose +* ======= +* +* DLATRS solves one of the triangular systems +* +* A *x = s*b or A'*x = s*b +* +* with scaling to prevent overflow. Here A is an upper or lower +* triangular matrix, A' denotes the transpose of A, x and b are +* n-element vectors, and s is a scaling factor, usually less than +* or equal to 1, chosen so that the components of x will be less than +* the overflow threshold. If the unscaled problem will not cause +* overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A +* is singular (A(j,j) = 0 for some j), then s is set to 0 and a +* non-trivial solution to A*x = 0 is returned. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the matrix A is upper or lower triangular. +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* TRANS (input) CHARACTER*1 +* Specifies the operation applied to A. +* = 'N': Solve A * x = s*b (No transpose) +* = 'T': Solve A'* x = s*b (Transpose) +* = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose) +* +* DIAG (input) CHARACTER*1 +* Specifies whether or not the matrix A is unit triangular. +* = 'N': Non-unit triangular +* = 'U': Unit triangular +* +* NORMIN (input) CHARACTER*1 +* Specifies whether CNORM has been set or not. +* = 'Y': CNORM contains the column norms on entry +* = 'N': CNORM is not set on entry. On exit, the norms will +* be computed and stored in CNORM. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* The triangular matrix A. If UPLO = 'U', the leading n by n +* upper triangular part of the array A contains the upper +* triangular matrix, and the strictly lower triangular part of +* A is not referenced. If UPLO = 'L', the leading n by n lower +* triangular part of the array A contains the lower triangular +* matrix, and the strictly upper triangular part of A is not +* referenced. If DIAG = 'U', the diagonal elements of A are +* also not referenced and are assumed to be 1. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max (1,N). +* +* X (input/output) DOUBLE PRECISION array, dimension (N) +* On entry, the right hand side b of the triangular system. +* On exit, X is overwritten by the solution vector x. +* +* SCALE (output) DOUBLE PRECISION +* The scaling factor s for the triangular system +* A * x = s*b or A'* x = s*b. +* If SCALE = 0, the matrix A is singular or badly scaled, and +* the vector x is an exact or approximate solution to A*x = 0. +* +* CNORM (input or output) DOUBLE PRECISION array, dimension (N) +* +* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) +* contains the norm of the off-diagonal part of the j-th column +* of A. If TRANS = 'N', CNORM(j) must be greater than or equal +* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) +* must be greater than or equal to the 1-norm. +* +* If NORMIN = 'N', CNORM is an output argument and CNORM(j) +* returns the 1-norm of the offdiagonal part of the j-th column +* of A. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* +* Further Details +* ======= ======= +* +* A rough bound on x is computed; if that is less than overflow, DTRSV +* is called, otherwise, specific code is used which checks for possible +* overflow or divide-by-zero at every operation. +* +* A columnwise scheme is used for solving A*x = b. The basic algorithm +* if A is lower triangular is +* +* x[1:n] := b[1:n] +* for j = 1, ..., n +* x(j) := x(j) / A(j,j) +* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] +* end +* +* Define bounds on the components of x after j iterations of the loop: +* M(j) = bound on x[1:j] +* G(j) = bound on x[j+1:n] +* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. +* +* Then for iteration j+1 we have +* M(j+1) <= G(j) / | A(j+1,j+1) | +* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | +* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) +* +* where CNORM(j+1) is greater than or equal to the infinity-norm of +* column j+1 of A, not counting the diagonal. Hence +* +* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) +* 1<=i<=j +* and +* +* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) +* 1<=i< j +* +* Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the +* reciprocal of the largest M(j), j=1,..,n, is larger than +* max(underflow, 1/overflow). +* +* The bound on x(j) is also used to determine when a step in the +* columnwise method can be performed without fear of overflow. If +* the computed bound is greater than a large constant, x is scaled to +* prevent overflow, but if the bound overflows, x is set to 0, x(j) to +* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. +* +* Similarly, a row-wise scheme is used to solve A'*x = b. The basic +* algorithm for A upper triangular is +* +* for j = 1, ..., n +* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) +* end +* +* We simultaneously compute two bounds +* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j +* M(j) = bound on x(i), 1<=i<=j +* +* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we +* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. +* Then the bound on x(j) is +* +* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | +* +* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) +* 1<=i<=j +* +* and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater +* than max(underflow, 1/overflow). +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, HALF, ONE + PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL NOTRAN, NOUNIT, UPPER + INTEGER I, IMAX, J, JFIRST, JINC, JLAST + DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS, + $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER IDAMAX + DOUBLE PRECISION DASUM, DDOT, DLAMCH + EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH +* .. +* .. External Subroutines .. + EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + NOTRAN = LSAME( TRANS, 'N' ) + NOUNIT = LSAME( DIAG, 'N' ) +* +* Test the input parameters. +* + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. + $ LSAME( TRANS, 'C' ) ) THEN + INFO = -2 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -3 + ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. + $ LSAME( NORMIN, 'N' ) ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLATRS', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine machine dependent parameters to control overflow. +* + SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) + BIGNUM = ONE / SMLNUM + SCALE = ONE +* + IF( LSAME( NORMIN, 'N' ) ) THEN +* +* Compute the 1-norm of each column, not including the diagonal. +* + IF( UPPER ) THEN +* +* A is upper triangular. +* + DO 10 J = 1, N + CNORM( J ) = DASUM( J-1, A( 1, J ), 1 ) + 10 CONTINUE + ELSE +* +* A is lower triangular. +* + DO 20 J = 1, N - 1 + CNORM( J ) = DASUM( N-J, A( J+1, J ), 1 ) + 20 CONTINUE + CNORM( N ) = ZERO + END IF + END IF +* +* Scale the column norms by TSCAL if the maximum element in CNORM is +* greater than BIGNUM. +* + IMAX = IDAMAX( N, CNORM, 1 ) + TMAX = CNORM( IMAX ) + IF( TMAX.LE.BIGNUM ) THEN + TSCAL = ONE + ELSE + TSCAL = ONE / ( SMLNUM*TMAX ) + CALL DSCAL( N, TSCAL, CNORM, 1 ) + END IF +* +* Compute a bound on the computed solution vector to see if the +* Level 2 BLAS routine DTRSV can be used. +* + J = IDAMAX( N, X, 1 ) + XMAX = ABS( X( J ) ) + XBND = XMAX + IF( NOTRAN ) THEN +* +* Compute the growth in A * x = b. +* + IF( UPPER ) THEN + JFIRST = N + JLAST = 1 + JINC = -1 + ELSE + JFIRST = 1 + JLAST = N + JINC = 1 + END IF +* + IF( TSCAL.NE.ONE ) THEN + GROW = ZERO + GO TO 50 + END IF +* + IF( NOUNIT ) THEN +* +* A is non-unit triangular. +* +* Compute GROW = 1/G(j) and XBND = 1/M(j). +* Initially, G(0) = max{x(i), i=1,...,n}. +* + GROW = ONE / MAX( XBND, SMLNUM ) + XBND = GROW + DO 30 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 50 +* +* M(j) = G(j-1) / abs(A(j,j)) +* + TJJ = ABS( A( J, J ) ) + XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) + IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN +* +* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) +* + GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) + ELSE +* +* G(j) could overflow, set GROW to 0. +* + GROW = ZERO + END IF + 30 CONTINUE + GROW = XBND + ELSE +* +* A is unit triangular. +* +* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. +* + GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) + DO 40 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 50 +* +* G(j) = G(j-1)*( 1 + CNORM(j) ) +* + GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) + 40 CONTINUE + END IF + 50 CONTINUE +* + ELSE +* +* Compute the growth in A' * x = b. +* + IF( UPPER ) THEN + JFIRST = 1 + JLAST = N + JINC = 1 + ELSE + JFIRST = N + JLAST = 1 + JINC = -1 + END IF +* + IF( TSCAL.NE.ONE ) THEN + GROW = ZERO + GO TO 80 + END IF +* + IF( NOUNIT ) THEN +* +* A is non-unit triangular. +* +* Compute GROW = 1/G(j) and XBND = 1/M(j). +* Initially, M(0) = max{x(i), i=1,...,n}. +* + GROW = ONE / MAX( XBND, SMLNUM ) + XBND = GROW + DO 60 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 80 +* +* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) +* + XJ = ONE + CNORM( J ) + GROW = MIN( GROW, XBND / XJ ) +* +* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) +* + TJJ = ABS( A( J, J ) ) + IF( XJ.GT.TJJ ) + $ XBND = XBND*( TJJ / XJ ) + 60 CONTINUE + GROW = MIN( GROW, XBND ) + ELSE +* +* A is unit triangular. +* +* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. +* + GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) + DO 70 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 80 +* +* G(j) = ( 1 + CNORM(j) )*G(j-1) +* + XJ = ONE + CNORM( J ) + GROW = GROW / XJ + 70 CONTINUE + END IF + 80 CONTINUE + END IF +* + IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN +* +* Use the Level 2 BLAS solve if the reciprocal of the bound on +* elements of X is not too small. +* + CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 ) + ELSE +* +* Use a Level 1 BLAS solve, scaling intermediate results. +* + IF( XMAX.GT.BIGNUM ) THEN +* +* Scale X so that its components are less than or equal to +* BIGNUM in absolute value. +* + SCALE = BIGNUM / XMAX + CALL DSCAL( N, SCALE, X, 1 ) + XMAX = BIGNUM + END IF +* + IF( NOTRAN ) THEN +* +* Solve A * x = b +* + DO 110 J = JFIRST, JLAST, JINC +* +* Compute x(j) = b(j) / A(j,j), scaling x if necessary. +* + XJ = ABS( X( J ) ) + IF( NOUNIT ) THEN + TJJS = A( J, J )*TSCAL + ELSE + TJJS = TSCAL + IF( TSCAL.EQ.ONE ) + $ GO TO 100 + END IF + TJJ = ABS( TJJS ) + IF( TJJ.GT.SMLNUM ) THEN +* +* abs(A(j,j)) > SMLNUM: +* + IF( TJJ.LT.ONE ) THEN + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by 1/b(j). +* + REC = ONE / XJ + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF + X( J ) = X( J ) / TJJS + XJ = ABS( X( J ) ) + ELSE IF( TJJ.GT.ZERO ) THEN +* +* 0 < abs(A(j,j)) <= SMLNUM: +* + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM +* to avoid overflow when dividing by A(j,j). +* + REC = ( TJJ*BIGNUM ) / XJ + IF( CNORM( J ).GT.ONE ) THEN +* +* Scale by 1/CNORM(j) to avoid overflow when +* multiplying x(j) times column j. +* + REC = REC / CNORM( J ) + END IF + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + X( J ) = X( J ) / TJJS + XJ = ABS( X( J ) ) + ELSE +* +* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and +* scale = 0, and compute a solution to A*x = 0. +* + DO 90 I = 1, N + X( I ) = ZERO + 90 CONTINUE + X( J ) = ONE + XJ = ONE + SCALE = ZERO + XMAX = ZERO + END IF + 100 CONTINUE +* +* Scale x if necessary to avoid overflow when adding a +* multiple of column j of A. +* + IF( XJ.GT.ONE ) THEN + REC = ONE / XJ + IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN +* +* Scale x by 1/(2*abs(x(j))). +* + REC = REC*HALF + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + END IF + ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN +* +* Scale x by 1/2. +* + CALL DSCAL( N, HALF, X, 1 ) + SCALE = SCALE*HALF + END IF +* + IF( UPPER ) THEN + IF( J.GT.1 ) THEN +* +* Compute the update +* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) +* + CALL DAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X, + $ 1 ) + I = IDAMAX( J-1, X, 1 ) + XMAX = ABS( X( I ) ) + END IF + ELSE + IF( J.LT.N ) THEN +* +* Compute the update +* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) +* + CALL DAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1, + $ X( J+1 ), 1 ) + I = J + IDAMAX( N-J, X( J+1 ), 1 ) + XMAX = ABS( X( I ) ) + END IF + END IF + 110 CONTINUE +* + ELSE +* +* Solve A' * x = b +* + DO 160 J = JFIRST, JLAST, JINC +* +* Compute x(j) = b(j) - sum A(k,j)*x(k). +* k<>j +* + XJ = ABS( X( J ) ) + USCAL = TSCAL + REC = ONE / MAX( XMAX, ONE ) + IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN +* +* If x(j) could overflow, scale x by 1/(2*XMAX). +* + REC = REC*HALF + IF( NOUNIT ) THEN + TJJS = A( J, J )*TSCAL + ELSE + TJJS = TSCAL + END IF + TJJ = ABS( TJJS ) + IF( TJJ.GT.ONE ) THEN +* +* Divide by A(j,j) when scaling x if A(j,j) > 1. +* + REC = MIN( ONE, REC*TJJ ) + USCAL = USCAL / TJJS + END IF + IF( REC.LT.ONE ) THEN + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF +* + SUMJ = ZERO + IF( USCAL.EQ.ONE ) THEN +* +* If the scaling needed for A in the dot product is 1, +* call DDOT to perform the dot product. +* + IF( UPPER ) THEN + SUMJ = DDOT( J-1, A( 1, J ), 1, X, 1 ) + ELSE IF( J.LT.N ) THEN + SUMJ = DDOT( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) + END IF + ELSE +* +* Otherwise, use in-line code for the dot product. +* + IF( UPPER ) THEN + DO 120 I = 1, J - 1 + SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) + 120 CONTINUE + ELSE IF( J.LT.N ) THEN + DO 130 I = J + 1, N + SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) + 130 CONTINUE + END IF + END IF +* + IF( USCAL.EQ.TSCAL ) THEN +* +* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) +* was not used to scale the dotproduct. +* + X( J ) = X( J ) - SUMJ + XJ = ABS( X( J ) ) + IF( NOUNIT ) THEN + TJJS = A( J, J )*TSCAL + ELSE + TJJS = TSCAL + IF( TSCAL.EQ.ONE ) + $ GO TO 150 + END IF +* +* Compute x(j) = x(j) / A(j,j), scaling if necessary. +* + TJJ = ABS( TJJS ) + IF( TJJ.GT.SMLNUM ) THEN +* +* abs(A(j,j)) > SMLNUM: +* + IF( TJJ.LT.ONE ) THEN + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale X by 1/abs(x(j)). +* + REC = ONE / XJ + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF + X( J ) = X( J ) / TJJS + ELSE IF( TJJ.GT.ZERO ) THEN +* +* 0 < abs(A(j,j)) <= SMLNUM: +* + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. +* + REC = ( TJJ*BIGNUM ) / XJ + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + X( J ) = X( J ) / TJJS + ELSE +* +* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and +* scale = 0, and compute a solution to A'*x = 0. +* + DO 140 I = 1, N + X( I ) = ZERO + 140 CONTINUE + X( J ) = ONE + SCALE = ZERO + XMAX = ZERO + END IF + 150 CONTINUE + ELSE +* +* Compute x(j) := x(j) / A(j,j) - sumj if the dot +* product has already been divided by 1/A(j,j). +* + X( J ) = X( J ) / TJJS - SUMJ + END IF + XMAX = MAX( XMAX, ABS( X( J ) ) ) + 160 CONTINUE + END IF + SCALE = SCALE / TSCAL + END IF +* +* Scale the column norms by 1/TSCAL for return. +* + IF( TSCAL.NE.ONE ) THEN + CALL DSCAL( N, ONE / TSCAL, CNORM, 1 ) + END IF +* + RETURN +* +* End of DLATRS +* + END diff --git a/ext/lapack/dtrcon.f b/ext/lapack/dtrcon.f new file mode 100644 index 000000000..8da58c760 --- /dev/null +++ b/ext/lapack/dtrcon.f @@ -0,0 +1,193 @@ + SUBROUTINE DTRCON( NORM, UPLO, DIAG, N, A, LDA, RCOND, WORK, + $ IWORK, INFO ) +* +* -- LAPACK routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* March 31, 1993 +* +* .. Scalar Arguments .. + CHARACTER DIAG, NORM, UPLO + INTEGER INFO, LDA, N + DOUBLE PRECISION RCOND +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DTRCON estimates the reciprocal of the condition number of a +* triangular matrix A, in either the 1-norm or the infinity-norm. +* +* The norm of A is computed and an estimate is obtained for +* norm(inv(A)), then the reciprocal of the condition number is +* computed as +* RCOND = 1 / ( norm(A) * norm(inv(A)) ). +* +* Arguments +* ========= +* +* NORM (input) CHARACTER*1 +* Specifies whether the 1-norm condition number or the +* infinity-norm condition number is required: +* = '1' or 'O': 1-norm; +* = 'I': Infinity-norm. +* +* UPLO (input) CHARACTER*1 +* = 'U': A is upper triangular; +* = 'L': A is lower triangular. +* +* DIAG (input) CHARACTER*1 +* = 'N': A is non-unit triangular; +* = 'U': A is unit triangular. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* The triangular matrix A. If UPLO = 'U', the leading N-by-N +* upper triangular part of the array A contains the upper +* triangular matrix, and the strictly lower triangular part of +* A is not referenced. If UPLO = 'L', the leading N-by-N lower +* triangular part of the array A contains the lower triangular +* matrix, and the strictly upper triangular part of A is not +* referenced. If DIAG = 'U', the diagonal elements of A are +* also not referenced and are assumed to be 1. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* RCOND (output) DOUBLE PRECISION +* The reciprocal of the condition number of the matrix A, +* computed as RCOND = 1/(norm(A) * norm(inv(A))). +* +* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) +* +* IWORK (workspace) INTEGER array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL NOUNIT, ONENRM, UPPER + CHARACTER NORMIN + INTEGER IX, KASE, KASE1 + DOUBLE PRECISION AINVNM, ANORM, SCALE, SMLNUM, XNORM +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER IDAMAX + DOUBLE PRECISION DLAMCH, DLANTR + EXTERNAL LSAME, IDAMAX, DLAMCH, DLANTR +* .. +* .. External Subroutines .. + EXTERNAL DLACON, DLATRS, DRSCL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) + NOUNIT = LSAME( DIAG, 'N' ) +* + IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN + INFO = -1 + ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -2 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DTRCON', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) THEN + RCOND = ONE + RETURN + END IF +* + RCOND = ZERO + SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( 1, N ) ) +* +* Compute the norm of the triangular matrix A. +* + ANORM = DLANTR( NORM, UPLO, DIAG, N, N, A, LDA, WORK ) +* +* Continue only if ANORM > 0. +* + IF( ANORM.GT.ZERO ) THEN +* +* Estimate the norm of the inverse of A. +* + AINVNM = ZERO + NORMIN = 'N' + IF( ONENRM ) THEN + KASE1 = 1 + ELSE + KASE1 = 2 + END IF + KASE = 0 + 10 CONTINUE + CALL DLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE ) + IF( KASE.NE.0 ) THEN + IF( KASE.EQ.KASE1 ) THEN +* +* Multiply by inv(A). +* + CALL DLATRS( UPLO, 'No transpose', DIAG, NORMIN, N, A, + $ LDA, WORK, SCALE, WORK( 2*N+1 ), INFO ) + ELSE +* +* Multiply by inv(A'). +* + CALL DLATRS( UPLO, 'Transpose', DIAG, NORMIN, N, A, LDA, + $ WORK, SCALE, WORK( 2*N+1 ), INFO ) + END IF + NORMIN = 'Y' +* +* Multiply by 1/SCALE if doing so will not cause overflow. +* + IF( SCALE.NE.ONE ) THEN + IX = IDAMAX( N, WORK, 1 ) + XNORM = ABS( WORK( IX ) ) + IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO ) + $ GO TO 20 + CALL DRSCL( N, SCALE, WORK, 1 ) + END IF + GO TO 10 + END IF +* +* Compute the estimate of the reciprocal condition number. +* + IF( AINVNM.NE.ZERO ) + $ RCOND = ( ONE / ANORM ) / AINVNM + END IF +* + 20 CONTINUE + RETURN +* +* End of DTRCON +* + END