Added in condition number routines for triangular matrices
This commit is contained in:
Harry Moffat
parent
7f485a0ba3
commit
24ac9bd133
8 changed files with 2424 additions and 0 deletions
|
|
@ -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 \
|
||||
|
|
|
|||
257
ext/f2c_lapack/dlacon.c
Normal file
257
ext/f2c_lapack/dlacon.c
Normal file
|
|
@ -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_ */
|
||||
|
||||
820
ext/f2c_lapack/dlatrs.c
Normal file
820
ext/f2c_lapack/dlatrs.c
Normal file
|
|
@ -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_ */
|
||||
|
||||
242
ext/f2c_lapack/dtrcon.c
Normal file
242
ext/f2c_lapack/dtrcon.c
Normal file
|
|
@ -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_ */
|
||||
|
||||
|
|
@ -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 \
|
||||
|
|
|
|||
204
ext/lapack/dlacon.f
Normal file
204
ext/lapack/dlacon.f
Normal file
|
|
@ -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
|
||||
702
ext/lapack/dlatrs.f
Normal file
702
ext/lapack/dlatrs.f
Normal file
|
|
@ -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
|
||||
193
ext/lapack/dtrcon.f
Normal file
193
ext/lapack/dtrcon.f
Normal file
|
|
@ -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
|
||||
Loading…
Add table
Reference in a new issue