//.........这里部分代码省略.........
/* generator. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Save statement .. */
/* .. */
/* .. Data statements .. */
/* .. */
/* .. Executable Statements .. */
/* Set some constants for use in the subroutine. */
if (first) {
first = FALSE_;
eps = slamch_("Precision");
badc2 = .1f / eps;
badc1 = sqrt(badc2);
small = slamch_("Safe minimum");
large = 1.f / small;
/* If it looks like we're on a Cray, take the square root of */
/* SMALL and LARGE to avoid overflow and underflow problems. */
slabad_(&small, &large);
small = small / eps * .25f;
large = 1.f / small;
}
s_copy(c2, path + 1, (ftnlen)2, (ftnlen)2);
/* Set some parameters we don't plan to change. */
*(unsigned char *)dist = 'S';
*mode = 3;
/* xQR, xLQ, xQL, xRQ: Set parameters to generate a general */
/* M x N matrix. */
if (lsamen_(&c__2, c2, "QR") || lsamen_(&c__2, c2,
"LQ") || lsamen_(&c__2, c2, "QL") || lsamen_(&c__2, c2, "RQ")) {
/* Set TYPE, the type of matrix to be generated. */
*(unsigned char *)type__ = 'N';
/* Set the lower and upper bandwidths. */
if (*imat == 1) {
*kl = 0;
*ku = 0;
} else if (*imat == 2) {
*kl = 0;
/* Computing MAX */
i__1 = *n - 1;
*ku = max(i__1,0);
} else if (*imat == 3) {
/* Subroutine */ int stbt03_(char *uplo, char *trans, char *diag, integer *n,
integer *kd, integer *nrhs, real *ab, integer *ldab, real *scale,
real *cnorm, real *tscal, real *x, integer *ldx, real *b, integer *
ldb, real *work, real *resid)
{
/* System generated locals */
integer ab_dim1, ab_offset, b_dim1, b_offset, x_dim1, x_offset, i__1;
real r__1, r__2, r__3;
/* Local variables */
static integer j;
extern logical lsame_(char *, char *);
extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
static real xscal;
extern /* Subroutine */ int stbmv_(char *, char *, char *, integer *,
integer *, real *, integer *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *);
static real tnorm, xnorm;
extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *,
real *, integer *), slabad_(real *, real *);
static integer ix;
extern doublereal slamch_(char *);
static real bignum;
extern integer isamax_(integer *, real *, integer *);
static real smlnum, eps, err;
#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
#define x_ref(a_1,a_2) x[(a_2)*x_dim1 + a_1]
#define ab_ref(a_1,a_2) ab[(a_2)*ab_dim1 + a_1]
/* -- LAPACK test routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
February 29, 1992
Purpose
=======
STBT03 computes the residual for the solution to a scaled triangular
system of equations A*x = s*b or A'*x = s*b when A is a
triangular band matrix. Here A' is the transpose of A, s is a scalar,
and x and b are N by NRHS matrices. The test ratio is the maximum
over the number of right hand sides of
norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
where op(A) denotes A or A' and EPS is the machine epsilon.
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': A *x = b (No transpose)
= 'T': A'*x = b (Transpose)
= 'C': A'*x = 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
N (input) INTEGER
The order of the matrix A. N >= 0.
KD (input) INTEGER
The number of superdiagonals or subdiagonals of the
triangular band matrix A. KD >= 0.
NRHS (input) INTEGER
The number of right hand sides, i.e., the number of columns
of the matrices X and B. NRHS >= 0.
AB (input) REAL array, dimension (LDAB,N)
The upper or lower triangular band matrix A, stored in the
first kd+1 rows of the array. The j-th column of A is stored
in the j-th column of the array AB as follows:
if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
LDAB (input) INTEGER
The leading dimension of the array AB. LDAB >= KD+1.
SCALE (input) REAL
The scaling factor s used in solving the triangular system.
CNORM (input) REAL array, dimension (N)
The 1-norms of the columns of A, not counting the diagonal.
TSCAL (input) REAL
The scaling factor used in computing the 1-norms in CNORM.
CNORM actually contains the column norms of TSCAL*A.
X (input) REAL array, dimension (LDX,NRHS)
The computed solution vectors for the system of linear
//.........这里部分代码省略.........
开发者ID:zangel，项目名称:uquad，代码行数:101，代码来源:stbt03.c
示例8: sqrt
/* Subroutine */ int snaitr_(integer *ido, char *bmat, integer *n, integer *k,
integer *np, integer *nb, real *resid, real *rnorm, real *v, integer
*ldv, real *h__, integer *ldh, integer *ipntr, real *workd, integer *
info, ftnlen bmat_len)
{
/* Initialized data */
static logical first = TRUE_;
/* System generated locals */
integer h_dim1, h_offset, v_dim1, v_offset, i__1, i__2;
real r__1, r__2;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static integer i__, j;
static real t0, t1, t2, t3, t4, t5;
static integer jj, ipj, irj, ivj;
static real ulp, tst1;
static integer ierr, iter;
static real unfl, ovfl;
extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
static integer itry;
static real temp1;
static logical orth1, orth2, step3, step4;
extern doublereal snrm2_(integer *, real *, integer *);
static real betaj;
extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
static integer infol;
extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *,
real *, integer *, real *, integer *, real *, real *, integer *,
ftnlen);
static real xtemp[2];
extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
integer *);
static real wnorm;
extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *,
real *, integer *), ivout_(integer *, integer *, integer *,
integer *, char *, ftnlen), smout_(integer *, integer *, integer *
, real *, integer *, integer *, char *, ftnlen), svout_(integer *,
integer *, real *, integer *, char *, ftnlen), sgetv0_(integer *,
char *, integer *, logical *, integer *, integer *, real *,
integer *, real *, real *, integer *, real *, integer *, ftnlen);
static real rnorm1;
extern /* Subroutine */ int slabad_(real *, real *);
extern doublereal slamch_(char *, ftnlen);
extern /* Subroutine */ int second_(real *), slascl_(char *, integer *,
integer *, real *, real *, integer *, integer *, real *, integer *
, integer *, ftnlen);
static logical rstart;
static integer msglvl;
static real smlnum;
extern doublereal slanhs_(char *, integer *, real *, integer *, real *,
ftnlen);
/* %----------------------------------------------------% */
/* | Include files for debugging and timing information | */
/* %----------------------------------------------------% */
/* \SCCS Information: @(#) */
/* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
/* %---------------------------------% */
/* | See debug.doc for documentation | */
/* %---------------------------------% */
/* %------------------% */
/* | Scalar Arguments | */
/* %------------------% */
/* %--------------------------------% */
/* | See stat.doc for documentation | */
/* %--------------------------------% */
/* \SCCS Information: @(#) */
/* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
/* %-----------------% */
/* | Array Arguments | */
/* %-----------------% */
/* %------------% */
/* | Parameters | */
/* %------------% */
/* %---------------% */
/* | Local Scalars | */
/* %---------------% */
/* %-----------------------% */
/* | Local Array Arguments | */
//.........这里部分代码省略.........
/* Subroutine */ int clahqr_(logical *wantt, logical *wantz, integer *n,
integer *ilo, integer *ihi, complex *h__, integer *ldh, complex *w,
integer *iloz, integer *ihiz, complex *z__, integer *ldz, integer *
info)
{
/* System generated locals */
integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4;
real r__1, r__2, r__3, r__4, r__5, r__6;
complex q__1, q__2, q__3, q__4, q__5, q__6, q__7;
/* Builtin functions */
double r_imag(complex *);
void r_cnjg(complex *, complex *);
double c_abs(complex *);
void c_sqrt(complex *, complex *), pow_ci(complex *, complex *, integer *)
;
/* Local variables */
integer i__, j, k, l, m;
real s;
complex t, u, v[2], x, y;
integer i1, i2;
complex t1;
real t2;
complex v2;
real aa, ab, ba, bb, h10;
complex h11;
real h21;
complex h22, sc;
integer nh, nz;
real sx;
integer jhi;
complex h11s;
integer jlo, its;
real ulp;
complex sum;
real tst;
complex temp;
extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
integer *), ccopy_(integer *, complex *, integer *, complex *,
integer *);
real rtemp;
extern /* Subroutine */ int slabad_(real *, real *), clarfg_(integer *,
complex *, complex *, integer *, complex *);
extern /* Complex */ VOID cladiv_(complex *, complex *, complex *);
extern doublereal slamch_(char *);
real safmin, safmax, smlnum;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CLAHQR is an auxiliary routine called by CHSEQR to update the */
/* eigenvalues and Schur decomposition already computed by CHSEQR, by */
/* dealing with the Hessenberg submatrix in rows and columns ILO to */
/* IHI. */
/* Arguments */
/* ========= */
/* WANTT (input) LOGICAL */
/* = .TRUE. : the full Schur form T is required; */
/* = .FALSE.: only eigenvalues are required. */
/* WANTZ (input) LOGICAL */
/* = .TRUE. : the matrix of Schur vectors Z is required; */
/* = .FALSE.: Schur vectors are not required. */
/* N (input) INTEGER */
/* The order of the matrix H. N >= 0. */
/* ILO (input) INTEGER */
/* IHI (input) INTEGER */
/* It is assumed that H is already upper triangular in rows and */
/* columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). */
/* CLAHQR works primarily with the Hessenberg submatrix in rows */
/* and columns ILO to IHI, but applies transformations to all of */
/* H if WANTT is .TRUE.. */
/* 1 <= ILO <= max(1,IHI); IHI <= N. */
/* H (input/output) COMPLEX array, dimension (LDH,N) */
/* On entry, the upper Hessenberg matrix H. */
/* On exit, if INFO is zero and if WANTT is .TRUE., then H */
/* is upper triangular in rows and columns ILO:IHI. If INFO */
/* is zero and if WANTT is .FALSE., then the contents of H */
/* are unspecified on exit. The output state of H in case */
/* INF is positive is below under the description of INFO. */
/* LDH (input) INTEGER */
/* The leading dimension of the array H. LDH >= max(1,N). */
//.........这里部分代码省略.........
/* Subroutine */ int ctrevc_(char *side, char *howmny, logical *select,
integer *n, complex *t, integer *ldt, complex *vl, integer *ldvl,
complex *vr, integer *ldvr, integer *mm, integer *m, complex *work,
real *rwork, integer *info)
{
/* System generated locals */
integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1,
i__2, i__3, i__4, i__5;
real r__1, r__2, r__3;
complex q__1, q__2;
/* Builtin functions */
double r_imag(complex *);
void r_cnjg(complex *, complex *);
/* Local variables */
integer i__, j, k, ii, ki, is;
real ulp;
logical allv;
real unfl, ovfl, smin;
logical over;
real scale;
extern logical lsame_(char *, char *);
extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
, complex *, integer *, complex *, integer *, complex *, complex *
, integer *);
real remax;
extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
complex *, integer *);
logical leftv, bothv, somev;
extern /* Subroutine */ int slabad_(real *, real *);
extern integer icamax_(integer *, complex *, integer *);
extern doublereal slamch_(char *);
extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer
*), xerbla_(char *, integer *), clatrs_(char *, char *,
char *, char *, integer *, complex *, integer *, complex *, real *
, real *, integer *);
extern doublereal scasum_(integer *, complex *, integer *);
logical rightv;
real smlnum;
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* CTREVC computes some or all of the right and/or left eigenvectors of */
/* a complex upper triangular matrix T. */
/* Matrices of this type are produced by the Schur factorization of */
/* a complex general matrix: A = Q*T*Q**H, as computed by CHSEQR. */
/* The right eigenvector x and the left eigenvector y of T corresponding */
/* to an eigenvalue w are defined by: */
/* T*x = w*x, (y**H)*T = w*(y**H) */
/* where y**H denotes the conjugate transpose of the vector y. */
/* The eigenvalues are not input to this routine, but are read directly */
/* from the diagonal of T. */
/* This routine returns the matrices X and/or Y of right and left */
/* eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an */
/* input matrix. If Q is the unitary factor that reduces a matrix A to */
/* Schur form T, then Q*X and Q*Y are the matrices of right and left */
/* eigenvectors of A. */
/* Arguments */
/* ========= */
/* SIDE (input) CHARACTER*1 */
/* = 'R': compute right eigenvectors only; */
/* = 'L': compute left eigenvectors only; */
/* = 'B': compute both right and left eigenvectors. */
/* HOWMNY (input) CHARACTER*1 */
/* = 'A': compute all right and/or left eigenvectors; */
/* = 'B': compute all right and/or left eigenvectors, */
/* backtransformed using the matrices supplied in */
/* VR and/or VL; */
/* = 'S': compute selected right and/or left eigenvectors, */
/* as indicated by the logical array SELECT. */
/* SELECT (input) LOGICAL array, dimension (N) */
/* If HOWMNY = 'S', SELECT specifies the eigenvectors to be */
/* computed. */
/* The eigenvector corresponding to the j-th eigenvalue is */
/* computed if SELECT(j) = .TRUE.. */
/* Not referenced if HOWMNY = 'A' or 'B'. */
/* N (input) INTEGER */
/* The order of the matrix T. N >= 0. */
//.........这里部分代码省略.........
/* Subroutine */ int stbt06_(real *rcond, real *rcondc, char *uplo, char *
diag, integer *n, integer *kd, real *ab, integer *ldab, real *work,
real *rat)
{
/* System generated locals */
integer ab_dim1, ab_offset;
real r__1, r__2;
/* Local variables */
static real rmin, rmax, anorm;
extern /* Subroutine */ int slabad_(real *, real *);
extern doublereal slamch_(char *);
static real bignum;
extern doublereal slantb_(char *, char *, char *, integer *, integer *,
real *, integer *, real *);
static real smlnum, eps;
/* -- LAPACK test routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
February 29, 1992
Purpose
=======
STBT06 computes a test ratio comparing RCOND (the reciprocal
condition number of a triangular matrix A) and RCONDC, the estimate
computed by STBCON. Information about the triangular matrix A is
used if one estimate is zero and the other is non-zero to decide if
underflow in the estimate is justified.
Arguments
=========
RCOND (input) REAL
The estimate of the reciprocal condition number obtained by
forming the explicit inverse of the matrix A and computing
RCOND = 1/( norm(A) * norm(inv(A)) ).
RCONDC (input) REAL
The estimate of the reciprocal condition number computed by
STBCON.
UPLO (input) CHARACTER
Specifies whether the matrix A is upper or lower triangular.
= 'U': Upper triangular
= 'L': Lower triangular
DIAG (input) CHARACTER
Specifies whether or not the matrix A is unit triangular.
= 'N': Non-unit triangular
= 'U': Unit triangular
N (input) INTEGER
The order of the matrix A. N >= 0.
KD (input) INTEGER
The number of superdiagonals or subdiagonals of the
triangular band matrix A. KD >= 0.
AB (input) REAL array, dimension (LDAB,N)
The upper or lower triangular band matrix A, stored in the
first kd+1 rows of the array. The j-th column of A is stored
in the j-th column of the array AB as follows:
if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
LDAB (input) INTEGER
The leading dimension of the array AB. LDAB >= KD+1.
WORK (workspace) REAL array, dimension (N)
RAT (output) REAL
The test ratio. If both RCOND and RCONDC are nonzero,
RAT = MAX( RCOND, RCONDC )/MIN( RCOND, RCONDC ) - 1.
If RAT = 0, the two estimates are exactly the same.
=====================================================================
Parameter adjustments */
ab_dim1 = *ldab;
ab_offset = 1 + ab_dim1 * 1;
ab -= ab_offset;
--work;
/* Function Body */
eps = slamch_("Epsilon");
rmax = dmax(*rcond,*rcondc);
rmin = dmin(*rcond,*rcondc);
/* Do the easy cases first. */
if (rmin < 0.f) {
/* Invalid value for RCOND or RCONDC, return 1/EPS. */
*rat = 1.f / eps;
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论