/* Subroutine */ int dsvdc_(doublereal *x, integer *ldx, integer *n, integer *
p, doublereal *s, doublereal *e, doublereal *u, integer *ldu,
doublereal *v, integer *ldv, doublereal *work, integer *job, integer *
info)
{
/* System generated locals */
integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2,
i__3;
doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7;
/* Builtin functions */
double d_sign(doublereal *, doublereal *), sqrt(doublereal);
/* Local variables */
static doublereal b, c__, f, g;
static integer i__, j, k, l, m;
static doublereal t, t1, el;
static integer kk;
static doublereal cs;
static integer ll, mm, ls;
static doublereal sl;
static integer lu;
static doublereal sm, sn;
static integer lm1, mm1, lp1, mp1, nct, ncu, lls, nrt;
static doublereal emm1, smm1;
static integer kase;
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
integer *);
static integer jobu, iter;
extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
doublereal *, integer *, doublereal *, doublereal *);
static doublereal test;
extern doublereal dnrm2_(integer *, doublereal *, integer *);
static integer nctp1, nrtp1;
extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
integer *);
static doublereal scale, shift;
extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
doublereal *, integer *), drotg_(doublereal *, doublereal *,
doublereal *, doublereal *);
static integer maxit;
extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *);
static logical wantu, wantv;
static doublereal ztest;
/* dsvdc is a subroutine to reduce a double precision nxp matrix x */
/* by orthogonal transformations u and v to diagonal form. the */
/* diagonal elements s(i) are the singular values of x. the */
/* columns of u are the corresponding left singular vectors, */
/* and the columns of v the right singular vectors. */
/* on entry */
/* x double precision(ldx,p), where ldx.ge.n. */
/* x contains the matrix whose singular value */
/* decomposition is to be computed. x is */
/* destroyed by dsvdc. */
/* ldx integer. */
/* ldx is the leading dimension of the array x. */
/* n integer. */
/* n is the number of rows of the matrix x. */
/* p integer. */
/* p is the number of columns of the matrix x. */
/* ldu integer. */
/* ldu is the leading dimension of the array u. */
/* (see below). */
/* ldv integer. */
/* ldv is the leading dimension of the array v. */
/* (see below). */
/* work double precision(n). */
/* work is a scratch array. */
/* job integer. */
/* job controls the computation of the singular */
/* vectors. it has the decimal expansion ab */
/* with the following meaning */
/* a.eq.0 do not compute the left singular */
/* vectors. */
/* a.eq.1 return the n left singular vectors */
/* in u. */
/* a.ge.2 return the first min(n,p) singular */
/* vectors in u. */
/* b.eq.0 do not compute the right singular */
/* vectors. */
/* b.eq.1 return the right singular vectors */
/* in v. */
/* on return */
/* s double precision(mm), where mm=min(n+1,p). */
//.........这里部分代码省略.........
/* Subroutine */ int dsytd2_(char* uplo, integer* n, doublereal* a, integer *
lda, doublereal* d__, doublereal* e, doublereal* tau, integer* info) {
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3;
/* Local variables */
integer i__;
extern doublereal ddot_(integer*, doublereal*, integer*, doublereal*,
integer*);
doublereal taui;
extern /* Subroutine */ int dsyr2_(char*, integer*, doublereal*,
doublereal*, integer*, doublereal*, integer*, doublereal*,
integer*);
doublereal alpha;
extern logical lsame_(char*, char*);
extern /* Subroutine */ int daxpy_(integer*, doublereal*, doublereal*,
integer*, doublereal*, integer*);
logical upper;
extern /* Subroutine */ int dsymv_(char*, integer*, doublereal*,
doublereal*, integer*, doublereal*, integer*, doublereal*,
doublereal*, integer*), dlarfg_(integer*, doublereal*,
doublereal*, integer*, doublereal*), xerbla_(char*, integer *
);
/* -- LAPACK routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal */
/* form T by an orthogonal similarity transformation: Q' * A * Q = T. */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER*1 */
/* Specifies whether the upper or lower triangular part of the */
/* symmetric matrix A is stored: */
/* = 'U': Upper triangular */
/* = 'L': Lower triangular */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
/* On entry, the symmetric matrix A. If UPLO = 'U', the leading */
/* n-by-n upper triangular part of A contains the upper */
/* triangular part of the matrix A, and the strictly lower */
/* triangular part of A is not referenced. If UPLO = 'L', the */
/* leading n-by-n lower triangular part of A contains the lower */
/* triangular part of the matrix A, and the strictly upper */
/* triangular part of A is not referenced. */
/* On exit, if UPLO = 'U', the diagonal and first superdiagonal */
/* of A are overwritten by the corresponding elements of the */
/* tridiagonal matrix T, and the elements above the first */
/* superdiagonal, with the array TAU, represent the orthogonal */
/* matrix Q as a product of elementary reflectors; if UPLO */
/* = 'L', the diagonal and first subdiagonal of A are over- */
/* written by the corresponding elements of the tridiagonal */
/* matrix T, and the elements below the first subdiagonal, with */
/* the array TAU, represent the orthogonal matrix Q as a product */
/* of elementary reflectors. See Further Details. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* D (output) DOUBLE PRECISION array, dimension (N) */
/* The diagonal elements of the tridiagonal matrix T: */
/* D(i) = A(i,i). */
/* E (output) DOUBLE PRECISION array, dimension (N-1) */
/* The off-diagonal elements of the tridiagonal matrix T: */
/* E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. */
/* TAU (output) DOUBLE PRECISION array, dimension (N-1) */
/* The scalar factors of the elementary reflectors (see Further */
/* Details). */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value. */
/* Further Details */
/* =============== */
/* If UPLO = 'U', the matrix Q is represented as a product of elementary */
/* reflectors */
/* Q = H(n-1) . . . H(2) H(1). */
/* Each H(i) has the form */
//.........这里部分代码省略.........
开发者ID:353，项目名称:viewercv，代码行数:101，代码来源:dsytd2.c
示例12: squares
//.........这里部分代码省略.........
The dimension of the array WORK. LWORK >= max(1,M+N+P).
For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
where NB is an upper bound for the optimal blocksizes for
DGEQRF, SGERQF, DORMQR and SORMRQ.
If LWORK = -1, then a workspace query is assumed; the routine
only calculates the optimal size of the WORK array, returns
this value as the first entry of the WORK array, and no error
message related to LWORK is issued by XERBLA.
INFO (output) INTEGER
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal value.
=====================================================================
Test the input parameters
Parameter adjustments */
/* Table of constant values */
static integer c__1 = 1;
static integer c_n1 = -1;
static doublereal c_b29 = -1.;
static doublereal c_b31 = 1.;
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
/* Local variables */
static integer lopt;
extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
doublereal *, doublereal *, integer *, doublereal *, integer *,
doublereal *, doublereal *, integer *), dcopy_(integer *,
doublereal *, integer *, doublereal *, integer *), daxpy_(integer
*, doublereal *, doublereal *, integer *, doublereal *, integer *)
, dtrmv_(char *, char *, char *, integer *, doublereal *, integer
*, doublereal *, integer *), dtrsv_(char *
, char *, char *, integer *, doublereal *, integer *, doublereal *
, integer *);
static integer nb, mn, nr;
extern /* Subroutine */ int dggrqf_(integer *, integer *, integer *,
doublereal *, integer *, doublereal *, doublereal *, integer *,
doublereal *, doublereal *, integer *, integer *), xerbla_(char *,
integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
static integer nb1, nb2, nb3, nb4;
extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *,
integer *, doublereal *, integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *, integer *),
dormrq_(char *, char *, integer *, integer *, integer *,
doublereal *, integer *, doublereal *, doublereal *, integer *,
doublereal *, integer *, integer *);
static integer lwkopt;
static logical lquery;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
b_dim1 = *ldb;
b_offset = 1 + b_dim1 * 1;
b -= b_offset;
--c__;
void syrkBFS8( double *C, double *A, int n, int x, int r, double alpha ) {
int nhalf = n/2;
int xNew = x/4;
int xNewer = x/8;
int rrank = getRelativeRank(x,xNew);
int rrank2 = getRelativeRank(xNew,xNewer);
int nOldTri = getSizeTri(r-1,x);
int nOldSq = getSizeSq(r-1,x);
double *C11 = C;
double *C21 = C + nOldTri;
double *C22 = C21 + nOldSq;
double *A11 = A;
double *A21 = A+nOldSq;
double *A12 = A21+nOldSq;
double *A22 = A12+nOldSq;
// first do the 4-way re-arrangement.
int nCSize;
if( rrank == 0 || rrank == 3 )
nCSize = 4*nOldTri;
else
nCSize = 4*nOldSq;
double *C21c = (double*) malloc( nOldSq*sizeof(double) );
//int Csizes[] = {nOldTri,nOldSq,0,nOldTri};
int Csizes2[] = {nOldTri,nOldSq,nOldSq,nOldTri};
double *nC = (double*) malloc( 4*Csizes2[rrank]*sizeof(double) );
//startTimer(TIMER_COMM_SYRK);
//reduceBy( 4, x, C1, nC, Csizes );
//stopTimer(TIMER_COMM_SYRK);
double *A1[] = {A11,A21,A22,A22};
double *A2[] = {A12,A11,A12,A21};
int Asizes[] = {nOldSq,nOldSq,nOldSq,nOldSq};
double *nA1 = (double*) malloc( 4*nOldSq*sizeof(double) );
double *nA2 = (double*) malloc( 4*nOldSq*sizeof(double) );
startTimer(TIMER_COMM_SYRK);
reduceBy( 4, x, A1, nA1, Asizes );
reduceBy( 4, x, A2, nA2, Asizes );
stopTimer(TIMER_COMM_SYRK);
if( rrank == 1 || rrank == 2 ) { // these two do the calls to mult
mult( nC, nA1, nA2, nhalf, xNew, r-1, 0. );
} else { // these two will do the recursive syrk calls. First, we need to split them up further
double *nCcopy = (double*) malloc( 4*nOldTri*sizeof(double) );
double *nnC = (double*) malloc( 8*nOldTri*sizeof(double) );
double *nnA = (double*) malloc( 8*nOldSq*sizeof(double) );
double *nnA1[] = {nA1,nA2};
int nAsizes[] = {4*nOldSq,4*nOldSq};
startTimer(TIMER_COMM_SYRK);
reduceBy( 2, xNew, nnA1, nnA, nAsizes );
stopTimer(TIMER_COMM_SYRK);
double *nnC1[] = {nC,nCcopy};
int nCsizes2[] = {4*nOldTri,4*nOldTri};
startTimer(TIMER_COMM_SYRK);
//reduceBy( 2, xNew, nnC1, nnC, nCsizes );
stopTimer(TIMER_COMM_SYRK);
syrk( nnC, nnA, nhalf, xNewer, r-1, 0. );
startTimer(TIMER_COMM_SYRK);
expandBy( 2, xNew, nnC1, nnC, nCsizes2 );
stopTimer(TIMER_COMM_SYRK);
// final additions
int ione = 1;
double done = 1.;
int s = 4*nOldTri;
daxpy_( &s, &done, nCcopy, &ione, nC, &ione );
free(nCcopy);
free(nnC);
free(nnA);
}
free(nA1);
free(nA2);
// recollect the answers, final additions
double *expC11, *expC21, *expC22;
if( alpha == 0. ) {
expC11 = C11;
expC21 = C21;
expC22 = C22;
} else {
expC11 = (double*) malloc( nOldTri*sizeof(double) );
expC21 = (double*) malloc( nOldSq*sizeof(double) );
expC22 = (double*) malloc( nOldTri*sizeof(double) );
}
double *C1[] = {expC11, expC21, C21c, expC22};
startTimer(TIMER_COMM_SYRK);
expandBy( 4, x, C1, nC, Csizes2 );
stopTimer(TIMER_COMM_SYRK);
free(nC);
int ione = 1;
double done = 1.;
if( alpha != 0 ) { // only correct for alpha=1
daxpy_( &nOldTri, &done, expC11, &ione, C11, &ione );
daxpy_( &nOldSq, &done, expC21, &ione, C21, &ione );
daxpy_( &nOldTri, &done, expC22, &ione, C22, &ione );
free(expC11);
free(expC21);
free(expC22);
}
//.........这里部分代码省略.........
开发者ID:benjamingr，项目名称:CAPS，代码行数:101，代码来源:syrk.cpp
示例14: ddot_
/* Subroutine */ int dspgst_(integer *itype, char *uplo, integer *n,
doublereal *ap, doublereal *bp, integer *info)
{
/* System generated locals */
integer i__1, i__2;
doublereal d__1;
/* Local variables */
integer j, k, j1, k1, jj, kk;
doublereal ct, ajj;
integer j1j1;
doublereal akk;
integer k1k1;
doublereal bjj, bkk;
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
integer *);
extern /* Subroutine */ int dspr2_(char *, integer *, doublereal *,
doublereal *, integer *, doublereal *, integer *, doublereal *), dscal_(integer *, doublereal *, doublereal *, integer *);
extern logical lsame_(char *, char *);
extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *), dspmv_(char *, integer *,
doublereal *, doublereal *, doublereal *, integer *, doublereal *,
doublereal *, integer *);
logical upper;
extern /* Subroutine */ int dtpmv_(char *, char *, char *, integer *,
doublereal *, doublereal *, integer *),
dtpsv_(char *, char *, char *, integer *, doublereal *,
doublereal *, integer *), xerbla_(char *,
integer *);
/* -- LAPACK routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DSPGST reduces a real symmetric-definite generalized eigenproblem */
/* to standard form, using packed storage. */
/* If ITYPE = 1, the problem is A*x = lambda*B*x, */
/* and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) */
/* If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or */
/* B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. */
/* B must have been previously factorized as U**T*U or L*L**T by DPPTRF. */
/* Arguments */
/* ========= */
/* ITYPE (input) INTEGER */
/* = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T); */
/* = 2 or 3: compute U*A*U**T or L**T*A*L. */
/* UPLO (input) CHARACTER*1 */
/* = 'U': Upper triangle of A is stored and B is factored as */
/* U**T*U; */
/* = 'L': Lower triangle of A is stored and B is factored as */
/* L*L**T. */
/* N (input) INTEGER */
/* The order of the matrices A and B. 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, if INFO = 0, the transformed matrix, stored in the */
/* same format as A. */
/* BP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
/* The triangular factor from the Cholesky factorization of B, */
/* stored in the same format as A, as returned by DPPTRF. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Executable Statements .. */
//.........这里部分代码省略.........
/* Subroutine */ int dtzrqf_(integer *m, integer *n, doublereal *a, integer *
lda, doublereal *tau, integer *info)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer i__, k, m1;
extern /* Subroutine */ int dger_(integer *, integer *, doublereal *,
doublereal *, integer *, doublereal *, integer *, doublereal *,
integer *), dgemv_(char *, integer *, integer *, doublereal *,
doublereal *, integer *, doublereal *, integer *, doublereal *,
doublereal *, integer *), dcopy_(integer *, doublereal *,
integer *, doublereal *, integer *), daxpy_(integer *, doublereal
*, doublereal *, integer *, doublereal *, integer *), dlarfg_(
integer *, doublereal *, doublereal *, integer *, doublereal *),
xerbla_(char *, integer *);
/* -- LAPACK routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* This routine is deprecated and has been replaced by routine DTZRZF. */
/* DTZRQF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A */
/* to upper triangular form by means of orthogonal transformations. */
/* The upper trapezoidal matrix A is factored as */
/* A = ( R 0 ) * Z, */
/* where Z is an N-by-N orthogonal matrix and R is an M-by-M upper */
/* triangular matrix. */
/* 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 >= M. */
/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
/* On entry, the leading M-by-N upper trapezoidal part of the */
/* array A must contain the matrix to be factorized. */
/* On exit, the leading M-by-M upper triangular part of A */
/* contains the upper triangular matrix R, and elements M+1 to */
/* N of the first M rows of A, with the array TAU, represent the */
/* orthogonal matrix Z as a product of M elementary reflectors. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,M). */
/* TAU (output) DOUBLE PRECISION array, dimension (M) */
/* The scalar factors of the elementary reflectors. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* Further Details */
/* =============== */
/* The factorization is obtained by Householder's method. The kth */
/* transformation matrix, Z( k ), which is used to introduce zeros into */
/* the ( m - k + 1 )th row of A, is given in the form */
/* Z( k ) = ( I 0 ), */
/* ( 0 T( k ) ) */
/* where */
/* T( k ) = I - tau*u( k )*u( k )', u( k ) = ( 1 ), */
/* ( 0 ) */
/* ( z( k ) ) */
/* tau is a scalar and z( k ) is an ( n - m ) element vector. */
/* tau and z( k ) are chosen to annihilate the elements of the kth row */
/* of X. */
/* The scalar tau is returned in the kth element of TAU and the vector */
/* u( k ) in the kth row of A, such that the elements of z( k ) are */
/* in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in */
/* the upper triangular part of A. */
/* Z is given by */
/* Z = Z( 1 ) * Z( 2 ) * ... * Z( m ). */
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论