/*! search the index of element having the largest absolute value
in 0-based numbering system */
inline void idamax(long& i, long& j, const _dgematrix& mat)
{ VERBOSE_REPORT;
long index( idamax_(mat.m*mat.n, mat.array, 1) -1 );
i =index%mat.m;
j =index/mat.m;
mat.destroy();
}
/*! search the index of element having the largest absolute value
in 0-based numbering system */
inline void idamax(long& i, long& j, const dssmatrix& mat)
{
#ifdef CPPL_VERBOSE
std::cerr << "# [MARK] idamax(long&, long&, const dssmatrix&)"
<< std::endl;
#endif//CPPL_VERBOSE
long index( idamax_(mat.VOL, mat.Array, 1) -1 );
i =mat.Indx[index];
j =mat.Jndx[index];
}
/*! return the index of element having the largest absolute value
in 0-based numbering system */
inline long idamax(const _drovector& vec)
{
#ifdef CPPL_VERBOSE
std::cerr << "# [MARK] idamax(const _drovector&)"
<< std::endl;
#endif//CPPL_VERBOSE
long i( idamax_(vec.L, vec.Array, 1) -1 );
vec.destroy();
return i;
}
/*! search the index of element having the largest absolute value
in 0-based numbering system */
inline void idamax(long& i, long& j, const dsymatrix& mat)
{VERBOSE_REPORT;
i=j=0;
double val(0.);
for(long J=0; J<mat.n; J++){
long I( J +idamax_(mat.m-J, mat.darray[J]+J, 1) -1 );
if( val<fabs(mat.darray[J][I]) ){
val =fabs(mat.darray[J][I]);
i=I;
j=J;
}
}
}
/* Subroutine */ int zptcon_(integer *n, doublereal *d__, doublecomplex *e,
doublereal *anorm, doublereal *rcond, doublereal *rwork, integer *
info)
{
/* -- LAPACK routine (version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
September 30, 1994
Purpose
=======
ZPTCON computes the reciprocal of the condition number (in the
1-norm) of a complex Hermitian positive definite tridiagonal matrix
using the factorization A = L*D*L**H or A = U**H*D*U computed by
ZPTTRF.
Norm(inv(A)) is computed by a direct method, and the reciprocal of
the condition number is computed as
RCOND = 1 / (ANORM * norm(inv(A))).
Arguments
=========
N (input) INTEGER
The order of the matrix A. N >= 0.
D (input) DOUBLE PRECISION array, dimension (N)
The n diagonal elements of the diagonal matrix D from the
factorization of A, as computed by ZPTTRF.
E (input) COMPLEX*16 array, dimension (N-1)
The (n-1) off-diagonal elements of the unit bidiagonal factor
U or L from the factorization of A, as computed by ZPTTRF.
ANORM (input) DOUBLE PRECISION
The 1-norm of the original matrix A.
RCOND (output) DOUBLE PRECISION
The reciprocal of the condition number of the matrix A,
computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
1-norm of inv(A) computed in this routine.
RWORK (workspace) DOUBLE PRECISION array, dimension (N)
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
Further Details
===============
The method used is described in Nicholas J. Higham, "Efficient
Algorithms for Computing the Condition Number of a Tridiagonal
Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
=====================================================================
Test the input arguments.
Parameter adjustments */
/* Table of constant values */
static integer c__1 = 1;
/* System generated locals */
integer i__1;
doublereal d__1;
/* Builtin functions */
double z_abs(doublecomplex *);
/* Local variables */
static integer i__, ix;
extern integer idamax_(integer *, doublereal *, integer *);
extern /* Subroutine */ int xerbla_(char *, integer *);
static doublereal ainvnm;
--rwork;
--e;
--d__;
/* Function Body */
*info = 0;
if (*n < 0) {
*info = -1;
} else if (*anorm < 0.) {
*info = -4;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("ZPTCON", &i__1);
return 0;
}
/* Quick return if possible */
*rcond = 0.;
if (*n == 0) {
*rcond = 1.;
//.........这里部分代码省略.........
/*< SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) >*/
/* Subroutine */ int dsptrf_(char *uplo, integer *n, doublereal *ap, integer *ipiv,
integer *info, ftnlen uplo_len)
{
/* System generated locals */
integer i__1;
doublereal d__1, d__2, d__3;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
doublereal c__;
integer j, k;
doublereal s, t, r1, r2;
integer kc, kk, kp, kx, knc, kpc=0, npp, imax=0, jmax;
extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
doublereal *, integer *, doublereal *, doublereal *), dspr_(char *
, integer *, doublereal *, doublereal *, integer *, doublereal *,
ftnlen);
doublereal alpha;
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
integer *);
extern logical lsame_(const char *, const char *, ftnlen, ftnlen);
extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
doublereal *, integer *);
integer kstep;
logical upper;
extern /* Subroutine */ int dlaev2_(doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *);
doublereal absakk;
extern integer idamax_(integer *, doublereal *, integer *);
extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
doublereal colmax, rowmax;
(void)uplo_len;
/* -- LAPACK routine (version 2.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* March 31, 1993 */
/* .. Scalar Arguments .. */
/*< CHARACTER UPLO >*/
/*< INTEGER INFO, N >*/
/* .. */
/* .. Array Arguments .. */
/*< INTEGER IPIV( * ) >*/
/*< DOUBLE PRECISION AP( * ) >*/
/* .. */
/* Purpose */
/* ======= */
/* DSPTRF computes the factorization of a real symmetric matrix A stored */
/* in packed format using the Bunch-Kaufman diagonal pivoting method: */
/* A = U*D*U**T or A = L*D*L**T */
/* where U (or L) is a product of permutation and unit upper (lower) */
/* triangular matrices, and D is symmetric and block diagonal with */
/* 1-by-1 and 2-by-2 diagonal blocks. */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER*1 */
/* = 'U': Upper triangle of A is stored; */
/* = 'L': Lower triangle of A is stored. */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
/* On entry, the upper or lower triangle of the symmetric matrix */
/* A, packed columnwise in a linear array. The j-th column of A */
/* is stored in the array AP as follows: */
/* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
/* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
/* On exit, the block diagonal matrix D and the multipliers used */
/* to obtain the factor U or L, stored as a packed triangular */
/* matrix overwriting A (see below for further details). */
/* IPIV (output) INTEGER array, dimension (N) */
/* Details of the interchanges and the block structure of D. */
/* If IPIV(k) > 0, then rows and columns k and IPIV(k) were */
/* interchanged and D(k,k) is a 1-by-1 diagonal block. */
/* If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and */
/* columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) */
/* is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = */
/* IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were */
/* interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* > 0: if INFO = i, D(i,i) is exactly zero. The factorization */
/* has been completed, but the block diagonal matrix D is */
/* exactly singular, and division by zero will occur if it */
//.........这里部分代码省略.........
/* Subroutine */ int dgbtrf_(integer *m, integer *n, integer *kl, integer *ku,
doublereal *ab, integer *ldab, integer *ipiv, integer *info)
{
/* System generated locals */
integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5, i__6;
doublereal d__1;
/* Local variables */
integer i__, j, i2, i3, j2, j3, k2, jb, nb, ii, jj, jm, ip, jp, km, ju,
kv, nw;
extern /* Subroutine */ int dger_(integer *, integer *, doublereal *,
doublereal *, integer *, doublereal *, integer *, doublereal *,
integer *);
doublereal temp;
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
integer *), dgemm_(char *, char *, integer *, integer *, integer *
, doublereal *, doublereal *, integer *, doublereal *, integer *,
doublereal *, doublereal *, integer *), dcopy_(
integer *, doublereal *, integer *, doublereal *, integer *),
dswap_(integer *, doublereal *, integer *, doublereal *, integer *
);
doublereal work13[4160] /* was [65][64] */, work31[4160] /*
was [65][64] */;
extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
integer *, integer *, doublereal *, doublereal *, integer *,
doublereal *, integer *), dgbtf2_(
integer *, integer *, integer *, integer *, doublereal *, integer
*, integer *, integer *);
extern integer idamax_(integer *, doublereal *, integer *);
extern /* Subroutine */ int xerbla_(char *, integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *);
extern /* Subroutine */ int dlaswp_(integer *, doublereal *, integer *,
integer *, integer *, integer *, integer *);
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGBTRF computes an LU factorization of a real m-by-n band matrix A */
/* using partial pivoting with row interchanges. */
/* This is the blocked version of the algorithm, calling Level 3 BLAS. */
/* Arguments */
/* ========= */
/* M (input) INTEGER */
/* The number of rows of the matrix A. M >= 0. */
/* N (input) INTEGER */
/* The number of columns of the matrix A. N >= 0. */
/* KL (input) INTEGER */
/* The number of subdiagonals within the band of A. KL >= 0. */
/* KU (input) INTEGER */
/* The number of superdiagonals within the band of A. KU >= 0. */
/* AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) */
/* On entry, the matrix A in band storage, in rows KL+1 to */
/* 2*KL+KU+1; rows 1 to KL of the array need not be set. */
/* The j-th column of A is stored in the j-th column of the */
/* array AB as follows: */
/* AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) */
/* On exit, details of the factorization: U is stored as an */
/* upper triangular band matrix with KL+KU superdiagonals in */
/* rows 1 to KL+KU+1, and the multipliers used during the */
/* factorization are stored in rows KL+KU+2 to 2*KL+KU+1. */
/* See below for further details. */
/* LDAB (input) INTEGER */
/* The leading dimension of the array AB. LDAB >= 2*KL+KU+1. */
/* IPIV (output) INTEGER array, dimension (min(M,N)) */
/* The pivot indices; for 1 <= i <= min(M,N), row i of the */
/* matrix was interchanged with row IPIV(i). */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* > 0: if INFO = +i, U(i,i) is exactly zero. The factorization */
/* has been completed, but the factor U is exactly */
/* singular, and division by zero will occur if it is used */
/* to solve a system of equations. */
/* Further Details */
/* =============== */
/* The band storage scheme is illustrated by the following example, when */
//.........这里部分代码省略.........
int
dlacon_(int *n, double *v, double *x, int *isgn, double *est, int *kase)
{
/*
Purpose
=======
DLACON estimates the 1-norm of a square matrix A.
Reverse communication is used for evaluating matrix-vector products.
Arguments
=========
N (input) INT
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) INT array, dimension (N)
EST (output) DOUBLE PRECISION
An estimate (a lower bound) for norm(A).
KASE (input/output) INT
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 CONEST, 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.
=====================================================================
*/
/* Table of constant values */
int c__1 = 1;
double zero = 0.0;
double one = 1.0;
/* Local variables */
static int iter;
static int jump, jlast;
static double altsgn, estold;
static int i, j;
double temp;
extern int idamax_(int *, double *, int *);
extern double dasum_(int *, double *, int *);
extern int dcopy_(int *, double *, int *, double *, int *);
#define d_sign(a, b) (b >= 0 ? fabs(a) : -fabs(a)) /* Copy sign */
#define i_dnnt(a) \
( a>=0 ? floor(a+.5) : -floor(.5-a) ) /* Round to nearest integer */
if ( *kase == 0 ) {
for (i = 0; i < *n; ++i) {
x[i] = 1. / (double) (*n);
}
*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[0] = x[0];
*est = fabs(v[0]);
/* ... QUIT */
goto L150;
}
*est = dasum_(n, x, &c__1);
for (i = 0; i < *n; ++i) {
x[i] = d_sign(one, x[i]);
//.........这里部分代码省略.........
/* Subroutine */ int zlaqps_(integer *m, integer *n, integer *offset, integer
*nb, integer *kb, doublecomplex *a, integer *lda, integer *jpvt,
doublecomplex *tau, doublereal *vn1, doublereal *vn2, doublecomplex *
auxv, doublecomplex *f, integer *ldf)
{
/* System generated locals */
integer a_dim1, a_offset, f_dim1, f_offset, i__1, i__2, i__3;
doublereal d__1, d__2;
doublecomplex z__1;
/* Builtin functions */
double sqrt(doublereal);
void d_cnjg(doublecomplex *, doublecomplex *);
double z_abs(doublecomplex *);
integer i_dnnt(doublereal *);
/* Local variables */
integer j, k, rk;
doublecomplex akk;
integer pvt;
doublereal temp, temp2, tol3z;
integer itemp;
extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *,
integer *, doublecomplex *, doublecomplex *, integer *,
doublecomplex *, integer *, doublecomplex *, doublecomplex *,
integer *), zgemv_(char *, integer *, integer *,
doublecomplex *, doublecomplex *, integer *, doublecomplex *,
integer *, doublecomplex *, doublecomplex *, integer *),
zswap_(integer *, doublecomplex *, integer *, doublecomplex *,
integer *);
extern doublereal dznrm2_(integer *, doublecomplex *, integer *), dlamch_(
char *);
extern integer idamax_(integer *, doublereal *, integer *);
integer lsticc;
extern /* Subroutine */ int zlarfp_(integer *, doublecomplex *,
doublecomplex *, integer *, doublecomplex *);
integer lastrk;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZLAQPS computes a step of QR factorization with column pivoting */
/* of a complex M-by-N matrix A by using Blas-3. It tries to factorize */
/* NB columns from A starting from the row OFFSET+1, and updates all */
/* of the matrix with Blas-3 xGEMM. */
/* In some cases, due to catastrophic cancellations, it cannot */
/* factorize NB columns. Hence, the actual number of factorized */
/* columns is returned in KB. */
/* Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. */
/* Arguments */
/* ========= */
/* M (input) INTEGER */
/* The number of rows of the matrix A. M >= 0. */
/* N (input) INTEGER */
/* The number of columns of the matrix A. N >= 0 */
/* OFFSET (input) INTEGER */
/* The number of rows of A that have been factorized in */
/* previous steps. */
/* NB (input) INTEGER */
/* The number of columns to factorize. */
/* KB (output) INTEGER */
/* The number of columns actually factorized. */
/* A (input/output) COMPLEX*16 array, dimension (LDA,N) */
/* On entry, the M-by-N matrix A. */
/* On exit, block A(OFFSET+1:M,1:KB) is the triangular */
/* factor obtained and block A(1:OFFSET,1:N) has been */
/* accordingly pivoted, but no factorized. */
/* The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has */
/* been updated. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,M). */
/* JPVT (input/output) INTEGER array, dimension (N) */
/* JPVT(I) = K <==> Column K of the full matrix A has been */
/* permuted into position I in AP. */
/* TAU (output) COMPLEX*16 array, dimension (KB) */
/* The scalar factors of the elementary reflectors. */
/* VN1 (input/output) DOUBLE PRECISION array, dimension (N) */
//.........这里部分代码省略.........
/* Subroutine */extern "C" int dlaqp2_(integer *m, integer *n, integer *offset,
doublereal *a, integer *lda, integer *jpvt, doublereal *tau,
doublereal *vn1, doublereal *vn2, doublereal *work)
{
/* -- 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, 1999
Purpose
=======
DLAQP2 computes a QR factorization with column pivoting of
the block A(OFFSET+1:M,1:N).
The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A. M >= 0.
N (input) INTEGER
The number of columns of the matrix A. N >= 0.
OFFSET (input) INTEGER
The number of rows of the matrix A that must be pivoted
but no factorized. OFFSET >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the M-by-N matrix A.
On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
the triangular factor obtained; the elements in block
A(OFFSET+1:M,1:N) below the diagonal, together with the
array TAU, represent the orthogonal matrix Q as a product of
elementary reflectors. Block A(1:OFFSET,1:N) has been
accordingly pivoted, but no factorized.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
JPVT (input/output) INTEGER array, dimension (N)
On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
to the front of A*P (a leading column); if JPVT(i) = 0,
the i-th column of A is a free column.
On exit, if JPVT(i) = k, then the i-th column of A*P
was the k-th column of A.
TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
The scalar factors of the elementary reflectors.
VN1 (input/output) DOUBLE PRECISION array, dimension (N)
The vector with the partial column norms.
VN2 (input/output) DOUBLE PRECISION array, dimension (N)
The vector with the exact column norms.
WORK (workspace) DOUBLE PRECISION array, dimension (N)
Further Details
===============
Based on contributions by
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
X. Sun, Computer Science Dept., Duke University, USA
=====================================================================
Parameter adjustments */
/* Table of constant values */
static integer c__1 = 1;
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
doublereal d__1, d__2;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
static doublereal temp;
extern doublereal dnrm2_(integer *, doublereal *, integer *);
static doublereal temp2;
static integer i__, j;
extern /* Subroutine */ int dlarf_(char *, integer *, integer *,
doublereal *, integer *, doublereal *, doublereal *, integer *,
doublereal *);
static integer offpi, itemp;
extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
doublereal *, integer *);
static integer mn;
extern /* Subroutine */ int dlarfg_(integer *, doublereal *, doublereal *,
integer *, doublereal *);
extern integer idamax_(integer *, doublereal *, integer *);
static doublereal aii;
static integer pvt;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论