本文整理汇总了C++中dgemm_函数的典型用法代码示例。如果您正苦于以下问题：C++ dgemm_函数的具体用法？C++ dgemm_怎么用？C++ dgemm_使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了dgemm_函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: dgemm_
/* Subroutine */ int dget03_(integer *n, doublereal *a, integer *lda,
doublereal *ainv, integer *ldainv, doublereal *work, integer *ldwork,
doublereal *rwork, doublereal *rcond, doublereal *resid)
{
/* System generated locals */
integer a_dim1, a_offset, ainv_dim1, ainv_offset, work_dim1, work_offset,
i__1;
/* Local variables */
integer i__;
doublereal eps;
extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
integer *, doublereal *, doublereal *, integer *, doublereal *,
integer *, doublereal *, doublereal *, integer *);
doublereal anorm;
extern doublereal dlamch_(char *), dlange_(char *, integer *,
integer *, doublereal *, integer *, doublereal *);
doublereal ainvnm;
/* -- LAPACK test routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGET03 computes the residual for a general matrix times its inverse: */
/* norm( I - AINV*A ) / ( N * norm(A) * norm(AINV) * EPS ), */
/* where EPS is the machine epsilon. */
/* Arguments */
/* ========== */
/* N (input) INTEGER */
/* The number of rows and columns of the matrix A. N >= 0. */
/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
/* The original N x N matrix A. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* AINV (input) DOUBLE PRECISION array, dimension (LDAINV,N) */
/* The inverse of the matrix A. */
/* LDAINV (input) INTEGER */
/* The leading dimension of the array AINV. LDAINV >= max(1,N). */
/* WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,N) */
/* LDWORK (input) INTEGER */
/* The leading dimension of the array WORK. LDWORK >= max(1,N). */
/* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */
/* RCOND (output) DOUBLE PRECISION */
/* The reciprocal of the condition number of A, computed as */
/* ( 1/norm(A) ) / norm(AINV). */
/* RESID (output) DOUBLE PRECISION */
/* norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS ) */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Quick exit if N = 0. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
ainv_dim1 = *ldainv;
ainv_offset = 1 + ainv_dim1;
ainv -= ainv_offset;
work_dim1 = *ldwork;
work_offset = 1 + work_dim1;
work -= work_offset;
--rwork;
/* Function Body */
if (*n <= 0) {
*rcond = 1.;
*resid = 0.;
//.........这里部分代码省略.........
示例2: dimension
//.........这里部分代码省略.........
The contents of A on exit are illustrated by the following examples:
m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
( v1 v2 v3 v4 v5 )
where d and e denote diagonal and off-diagonal elements of B, vi
denotes an element of the vector defining H(i), and ui an element of
the vector defining G(i).
=====================================================================
Test the input parameters
Parameter adjustments */
/* Table of constant values */
static integer c__1 = 1;
static integer c_n1 = -1;
static integer c__3 = 3;
static integer c__2 = 2;
static doublereal c_b21 = -1.;
static doublereal c_b22 = 1.;
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
/* Local variables */
static integer i__, j;
extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
integer *, doublereal *, doublereal *, integer *, doublereal *,
integer *, doublereal *, doublereal *, integer *);
static integer nbmin, iinfo, minmn;
extern /* Subroutine */ int dgebd2_(integer *, integer *, doublereal *,
integer *, doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *);
static integer nb;
extern /* Subroutine */ int dlabrd_(integer *, integer *, integer *,
doublereal *, integer *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, integer *, doublereal *, integer *);
static integer nx;
static doublereal ws;
extern /* Subroutine */ int xerbla_(char *, integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
static integer ldwrkx, ldwrky, lwkopt;
static logical lquery;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
a_dim1 = *lda;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
--d__;
--e;
--tauq;
--taup;
--work;
/* Function Body */
*info = 0;
/* Computing MAX */示例3: sqrt
void CheMPS2::TensorQ::AddTermsABLeft(TensorA * denA, TensorB * denB, TensorT * denT, double * workmem, double * workmem2){
for (int ikappa=0; ikappa<nKappa; ikappa++){
const int ID = denBK->directProd(sectorI1[ikappa],Idiff);
int dimLU = denBK->gCurrentDim(index, sectorN1[ikappa], sectorTwoS1[ikappa], sectorI1[ikappa]);
int dimLD = denBK->gCurrentDim(index, sectorN1[ikappa]+1, sectorTwoSD[ikappa], ID);
//case 1
const int IRD = denBK->directProd(ID, denBK->gIrrep(index));
for (int TwoSRD=sectorTwoSD[ikappa]-1; TwoSRD<=sectorTwoSD[ikappa]+1; TwoSRD+=2){
int dimRU = denBK->gCurrentDim(index+1, sectorN1[ikappa], sectorTwoS1[ikappa], sectorI1[ikappa]);
int dimRD = denBK->gCurrentDim(index+1, sectorN1[ikappa]+2, TwoSRD, IRD);
if ((dimRU>0) && (dimRD>0)){
int fase = ((((TwoSRD + sectorTwoS1[ikappa] + 2)/2)%2)!=0)?-1:1;
double factorB = fase * sqrt(3.0/(sectorTwoSD[ikappa]+1.0)) * (TwoSRD+1) * gsl_sf_coupling_6j(1,1,2,sectorTwoS1[ikappa],TwoSRD,sectorTwoSD[ikappa]);
double alpha;
double * mem;
if (TwoSRD == sectorTwoS1[ikappa]){
fase = ((((sectorTwoS1[ikappa]+1-sectorTwoSD[ikappa])/2)%2)!=0)?-1:1;
double factorA = fase * sqrt( 0.5 * (sectorTwoS1[ikappa]+1.0) / (sectorTwoSD[ikappa]+1.0) );
double * BlockA = denA->gStorage( sectorN1[ikappa], sectorTwoS1[ikappa], sectorI1[ikappa], sectorN1[ikappa]+2, TwoSRD, IRD );
double * BlockB = denB->gStorage( sectorN1[ikappa], sectorTwoS1[ikappa], sectorI1[ikappa], sectorN1[ikappa]+2, TwoSRD, IRD );
mem = workmem;
for (int cnt=0; cnt<dimRU*dimRD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
alpha = 1.0;
} else {
alpha = factorB;
mem = denB->gStorage( sectorN1[ikappa], sectorTwoS1[ikappa], sectorI1[ikappa], sectorN1[ikappa]+2, TwoSRD, IRD );
}
double * BlockTup = denT->gStorage(sectorN1[ikappa], sectorTwoS1[ikappa],sectorI1[ikappa],sectorN1[ikappa], sectorTwoS1[ikappa], sectorI1[ikappa]);
double * BlockTdo = denT->gStorage(sectorN1[ikappa]+1,sectorTwoSD[ikappa],ID, sectorN1[ikappa]+2,TwoSRD, IRD);
char notr = 'N';
double beta = 0.0; //set
dgemm_(¬r,¬r,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,mem,&dimRU,&beta,workmem2,&dimLU);
alpha = 1.0;
beta = 1.0; //add
char trans = 'T';
dgemm_(¬r,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
}
}
//case 2
const int IRU = denBK->directProd(sectorI1[ikappa],denBK->gIrrep(index));
for (int TwoSRU=sectorTwoS1[ikappa]-1; TwoSRU<=sectorTwoS1[ikappa]+1; TwoSRU+=2){
int dimRU = denBK->gCurrentDim(index+1, sectorN1[ikappa]+1, TwoSRU, IRU);
int dimRD = denBK->gCurrentDim(index+1, sectorN1[ikappa]+3, sectorTwoSD[ikappa], ID);
if ((dimRU>0) && (dimRD>0)){
int fase = ((((sectorTwoS1[ikappa] + sectorTwoSD[ikappa] + 1)/2)%2)!=0)?-1:1;
double factorB = fase * sqrt(3.0*(TwoSRU+1)) * gsl_sf_coupling_6j(1,1,2,TwoSRU,sectorTwoSD[ikappa],sectorTwoS1[ikappa]);
double alpha;
double * mem;
if (TwoSRU == sectorTwoSD[ikappa]){
double factorA = - sqrt(0.5);
double * BlockA = denA->gStorage( sectorN1[ikappa]+1, TwoSRU, IRU, sectorN1[ikappa]+3, sectorTwoSD[ikappa], ID );
double * BlockB = denB->gStorage( sectorN1[ikappa]+1, TwoSRU, IRU, sectorN1[ikappa]+3, sectorTwoSD[ikappa], ID );
mem = workmem;
for (int cnt=0; cnt<dimRU*dimRD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
alpha = 1.0;
} else {
alpha = factorB;
mem = denB->gStorage( sectorN1[ikappa]+1, TwoSRU, IRU, sectorN1[ikappa]+3, sectorTwoSD[ikappa], ID );
}
double * BlockTup = denT->gStorage(sectorN1[ikappa], sectorTwoS1[ikappa], sectorI1[ikappa], sectorN1[ikappa]+1, TwoSRU, IRU);
double * BlockTdo = denT->gStorage(sectorN1[ikappa]+1, sectorTwoSD[ikappa], ID, sectorN1[ikappa]+3, sectorTwoSD[ikappa], ID);
char notr = 'N';
double beta = 0.0; //set
dgemm_(¬r,¬r,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,mem,&dimRU,&beta,workmem2,&dimLU);
alpha = 1.0;
beta = 1.0; //add
char trans = 'T';
dgemm_(¬r,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
}
}
}
//.........这里部分代码省略.........
示例4: UPLO
/* Subroutine */ int dpbtrf_(char *uplo, integer *n, integer *kd, doublereal *
ab, integer *ldab, integer *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
Purpose
=======
DPBTRF computes the Cholesky factorization of a real symmetric
positive definite band matrix A.
The factorization has the form
A = U**T * U, if UPLO = 'U', or
A = L * L**T, if UPLO = 'L',
where U is an upper triangular matrix and L is lower triangular.
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.
KD (input) INTEGER
The number of superdiagonals of the matrix A if UPLO = 'U',
or the number of subdiagonals if UPLO = 'L'. KD >= 0.
AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
On entry, the upper or lower triangle of the symmetric 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).
On exit, if INFO = 0, the triangular factor U or L from the
Cholesky factorization A = U**T*U or A = L*L**T of the band
matrix A, in the same storage format as A.
LDAB (input) INTEGER
The leading dimension of the array AB. LDAB >= KD+1.
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, the leading minor of order i is not
positive definite, and the factorization could not be
completed.
Further Details
===============
The band storage scheme is illustrated by the following example, when
N = 6, KD = 2, and UPLO = 'U':
On entry: On exit:
* * a13 a24 a35 a46 * * u13 u24 u35 u46
* a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
Similarly, if UPLO = 'L' the format of A is as follows:
On entry: On exit:
a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66
a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 *
a31 a42 a53 a64 * * l31 l42 l53 l64 * *
Array elements marked * are not used by the routine.
Contributed by
Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
=====================================================================
Test the input parameters.
Parameter adjustments */
/* Table of constant values */
static integer c__1 = 1;
static integer c_n1 = -1;
static doublereal c_b18 = 1.;
static doublereal c_b21 = -1.;
static integer c__33 = 33;
/* System generated locals */
integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4;
/* Local variables */
static doublereal work[1056] /* was [33][32] */;
static integer i__, j;
extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
//.........这里部分代码省略.........
示例5: dlauum_
int dlauum_(char *uplo, int *n, double *a, int *
lda, int *info)
{
/* System generated locals */
int a_dim1, a_offset, i__1, i__2, i__3, i__4;
/* Local variables */
int i__, ib, nb;
extern int dgemm_(char *, char *, int *, int *,
int *, double *, double *, int *, double *,
int *, double *, double *, int *);
extern int lsame_(char *, char *);
extern int dtrmm_(char *, char *, char *, char *,
int *, int *, double *, double *, int *,
double *, int *);
int upper;
extern int dsyrk_(char *, char *, int *, int *,
double *, double *, int *, double *, double *,
int *), dlauu2_(char *, int *,
double *, int *, int *), xerbla_(char *,
int *);
extern int ilaenv_(int *, char *, char *, int *, int *,
int *, int *);
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLAUUM computes the product U * U' or L' * L, where the triangular */
/* factor U or L is stored in the upper or lower triangular part of */
/* the array A. */
/* If UPLO = 'U' or 'u' then the upper triangle of the result is stored, */
/* overwriting the factor U in A. */
/* If UPLO = 'L' or 'l' then the lower triangle of the result is stored, */
/* overwriting the factor L in A. */
/* This is the blocked form of the algorithm, calling Level 3 BLAS. */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER*1 */
/* Specifies whether the triangular factor stored in the array A */
/* is upper or lower triangular: */
/* = 'U': Upper triangular */
/* = 'L': Lower triangular */
/* N (input) INTEGER */
/* The order of the triangular factor U or L. N >= 0. */
/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
/* On entry, the triangular factor U or L. */
/* On exit, if UPLO = 'U', the upper triangle of A is */
/* overwritten with the upper triangle of the product U * U'; */
/* if UPLO = 'L', the lower triangle of A is overwritten with */
/* the lower triangle of the product L' * L. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= MAX(1,N). */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -k, the k-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;
/* Function Body */
*info = 0;
upper = lsame_(uplo, "U");
if (! upper && ! lsame_(uplo, "L")) {
*info = -1;
//.........这里部分代码省略.........
示例6: log
/* Subroutine */ int dlaed0_(integer *icompq, integer *qsiz, integer *n,
doublereal *d__, doublereal *e, doublereal *q, integer *ldq,
doublereal *qstore, integer *ldqs, doublereal *work, integer *iwork,
integer *info)
{
/* System generated locals */
integer q_dim1, q_offset, qstore_dim1, qstore_offset, i__1, i__2;
doublereal d__1;
/* Builtin functions */
double log(doublereal);
integer pow_ii(integer *, integer *);
/* Local variables */
static integer i__, j, k, iq, lgn, msd2, smm1, spm1, spm2;
static doublereal temp;
static integer curr;
extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
integer *, doublereal *, doublereal *, integer *, doublereal *,
integer *, doublereal *, doublereal *, integer *);
static integer iperm;
extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
doublereal *, integer *);
static integer indxq, iwrem;
extern /* Subroutine */ int dlaed1_(integer *, doublereal *, doublereal *,
integer *, integer *, doublereal *, integer *, doublereal *,
integer *, integer *);
static integer iqptr;
extern /* Subroutine */ int dlaed7_(integer *, integer *, integer *,
integer *, integer *, integer *, doublereal *, doublereal *,
integer *, integer *, doublereal *, integer *, doublereal *,
integer *, integer *, integer *, integer *, integer *, doublereal
*, doublereal *, integer *, integer *);
static integer tlvls;
extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
doublereal *, integer *, doublereal *, integer *);
static integer igivcl;
extern /* Subroutine */ int xerbla_(char *, integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
static integer igivnm, submat, curprb, subpbs, igivpt;
extern /* Subroutine */ int dsteqr_(char *, integer *, doublereal *,
doublereal *, doublereal *, integer *, doublereal *, integer *);
static integer curlvl, matsiz, iprmpt, smlsiz;
/* -- LAPACK routine (version 3.1) --
Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
November 2006
Purpose
=======
DLAED0 computes all eigenvalues and corresponding eigenvectors of a
symmetric tridiagonal matrix using the divide and conquer method.
Arguments
=========
ICOMPQ (input) INTEGER
= 0: Compute eigenvalues only.
= 1: Compute eigenvectors of original dense symmetric matrix
also. On entry, Q contains the orthogonal matrix used
to reduce the original matrix to tridiagonal form.
= 2: Compute eigenvalues and eigenvectors of tridiagonal
matrix.
QSIZ (input) INTEGER
The dimension of the orthogonal matrix used to reduce
the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
N (input) INTEGER
The dimension of the symmetric tridiagonal matrix. N >= 0.
D (input/output) DOUBLE PRECISION array, dimension (N)
On entry, the main diagonal of the tridiagonal matrix.
On exit, its eigenvalues.
E (input) DOUBLE PRECISION array, dimension (N-1)
The off-diagonal elements of the tridiagonal matrix.
On exit, E has been destroyed.
Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
On entry, Q must contain an N-by-N orthogonal matrix.
If ICOMPQ = 0 Q is not referenced.
If ICOMPQ = 1 On entry, Q is a subset of the columns of the
orthogonal matrix used to reduce the full
matrix to tridiagonal form corresponding to
the subset of the full matrix which is being
decomposed at this time.
If ICOMPQ = 2 On entry, Q will be the identity matrix.
On exit, Q contains the eigenvectors of the
tridiagonal matrix.
LDQ (input) INTEGER
The leading dimension of the array Q. If eigenvectors are
desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
//.........这里部分代码省略.........
示例7: dlamch_
//.........这里部分代码省略.........
result[1] = 0.;
result[2] = 0.;
result[3] = 0.;
result[4] = 0.;
return 0;
}
/* Copy the last k columns of the factorization to the array Q */
dlaset_("Full", m, m, &c_b4, &c_b4, &q[q_offset], lda);
if (*k > 0 && *m > *k) {
i__1 = *m - *k;
dlacpy_("Full", &i__1, k, &af[(*n - *k + 1) * af_dim1 + 1], lda, &q[(*
m - *k + 1) * q_dim1 + 1], lda);
}
if (*k > 1) {
i__1 = *k - 1;
i__2 = *k - 1;
dlacpy_("Upper", &i__1, &i__2, &af[*m - *k + 1 + (*n - *k + 2) *
af_dim1], lda, &q[*m - *k + 1 + (*m - *k + 2) * q_dim1], lda);
}
/* Generate the m-by-m matrix Q */
s_copy(srnamc_1.srnamt, "DORGQL", (ftnlen)32, (ftnlen)6);
dorgql_(m, m, k, &q[q_offset], lda, &tau[minmn - *k + 1], &work[1], lwork,
&info);
for (iside = 1; iside <= 2; ++iside) {
if (iside == 1) {
*(unsigned char *)side = 'L';
mc = *m;
nc = *n;
} else {
*(unsigned char *)side = 'R';
mc = *n;
nc = *m;
}
/* Generate MC by NC matrix C */
i__1 = nc;
for (j = 1; j <= i__1; ++j) {
dlarnv_(&c__2, iseed, &mc, &c__[j * c_dim1 + 1]);
/* L10: */
}
cnorm = dlange_("1", &mc, &nc, &c__[c_offset], lda, &rwork[1]);
if (cnorm == 0.) {
cnorm = 1.;
}
for (itrans = 1; itrans <= 2; ++itrans) {
if (itrans == 1) {
*(unsigned char *)trans = 'N';
} else {
*(unsigned char *)trans = 'T';
}
/* Copy C */
dlacpy_("Full", &mc, &nc, &c__[c_offset], lda, &cc[cc_offset],
lda);
/* Apply Q or Q' to C */
s_copy(srnamc_1.srnamt, "DORMQL", (ftnlen)32, (ftnlen)6);
if (*k > 0) {
dormql_(side, trans, &mc, &nc, k, &af[(*n - *k + 1) * af_dim1
+ 1], lda, &tau[minmn - *k + 1], &cc[cc_offset], lda,
&work[1], lwork, &info);
}
/* Form explicit product and subtract */
if (lsame_(side, "L")) {
dgemm_(trans, "No transpose", &mc, &nc, &mc, &c_b22, &q[
q_offset], lda, &c__[c_offset], lda, &c_b23, &cc[
cc_offset], lda);
} else {
dgemm_("No transpose", trans, &mc, &nc, &nc, &c_b22, &c__[
c_offset], lda, &q[q_offset], lda, &c_b23, &cc[
cc_offset], lda);
}
/* Compute error in the difference */
resid = dlange_("1", &mc, &nc, &cc[cc_offset], lda, &rwork[1]);
result[(iside - 1 << 1) + itrans] = resid / ((doublereal) max(1,*
m) * cnorm * eps);
/* L20: */
}
/* L30: */
}
return 0;
/* End of DQLT03 */
} /* dqlt03_ */示例8: main
//.........这里部分代码省略.........
#if NBLIS >= 6
{
bli_trsm( BLIS_LEFT,
&BLIS_ONE,
&a,
&c );
}
#endif
#endif
#ifdef NBLAS
#if NBLAS >= 1
for ( ii = 0; ii < 2000000000; ++ii )
{
f77_char transa = 'N';
f77_char transb = 'N';
f77_int mm = bli_obj_length( c );
f77_int kk = bli_obj_width_after_trans( a );
f77_int nn = bli_obj_width( c );
f77_int lda = bli_obj_col_stride( a );
f77_int ldb = bli_obj_col_stride( b );
f77_int ldc = bli_obj_col_stride( c );
double* alphap = bli_obj_buffer( alpha );
double* ap = bli_obj_buffer( a );
double* bp = bli_obj_buffer( b );
double* betap = bli_obj_buffer( beta );
double* cp = bli_obj_buffer( c );
dgemm_( &transa,
&transb,
&mm,
&nn,
&kk,
alphap,
ap, &lda,
bp, &ldb,
betap,
cp, &ldc );
}
#endif
#if NBLAS >= 2
{
f77_char side = 'L';
f77_char uplo = 'L';
f77_int mm = bli_obj_length( c );
f77_int nn = bli_obj_width( c );
f77_int lda = bli_obj_col_stride( a );
f77_int ldb = bli_obj_col_stride( b );
f77_int ldc = bli_obj_col_stride( c );
double* alphap = bli_obj_buffer( alpha );
double* ap = bli_obj_buffer( a );
double* bp = bli_obj_buffer( b );
double* betap = bli_obj_buffer( beta );
double* cp = bli_obj_buffer( c );
dsymm_( &side,
&uplo,
&mm,
&nn,
alphap,示例9: dgemm_
/* Subroutine */ int dgetrf_(integer* m, integer* n, doublereal* a, integer *
lda, integer* ipiv, integer* info) {
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
/* Local variables */
integer i__, j, jb, nb;
extern /* Subroutine */ int dgemm_(char*, char*, integer*, integer*,
integer*, doublereal*, doublereal*, integer*, doublereal*,
integer*, doublereal*, doublereal*, integer*);
integer iinfo;
extern /* Subroutine */ int dtrsm_(char*, char*, char*, char*,
integer*, integer*, doublereal*, doublereal*, integer*,
doublereal*, integer*), dgetf2_(
integer*, integer*, doublereal*, integer*, integer*, integer
*), 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.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGETRF computes an LU factorization of a general M-by-N matrix A */
/* using partial pivoting with row interchanges. */
/* The factorization has the form */
/* A = P * L * U */
/* where P is a permutation matrix, L is lower triangular with unit */
/* diagonal elements (lower trapezoidal if m > n), and U is upper */
/* triangular (upper trapezoidal if m < n). */
/* This is the right-looking Level 3 BLAS version of the algorithm. */
/* 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. */
/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
/* On entry, the M-by-N matrix to be factored. */
/* On exit, the factors L and U from the factorization */
/* A = P*L*U; the unit diagonal elements of L are not stored. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,M). */
/* 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. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--ipiv;
/* Function Body */
*info = 0;
if (*m < 0) {
*info = -1;
//.........这里部分代码省略.........
示例10: dgemm_
/* Subroutine */ int dgqrts_(integer *n, integer *m, integer *p, doublereal *
a, doublereal *af, doublereal *q, doublereal *r__, integer *lda,
doublereal *taua, doublereal *b, doublereal *bf, doublereal *z__,
doublereal *t, doublereal *bwk, integer *ldb, doublereal *taub,
doublereal *work, integer *lwork, doublereal *rwork, doublereal *
result)
{
/* System generated locals */
integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, bf_dim1,
bf_offset, bwk_dim1, bwk_offset, q_dim1, q_offset, r_dim1,
r_offset, t_dim1, t_offset, z_dim1, z_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
static integer info;
static doublereal unfl;
extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
integer *, doublereal *, doublereal *, integer *, doublereal *,
integer *, doublereal *, doublereal *, integer *);
static doublereal resid, anorm, bnorm;
extern /* Subroutine */ int dsyrk_(char *, char *, integer *, integer *,
doublereal *, doublereal *, integer *, doublereal *, doublereal *,
integer *);
extern doublereal dlamch_(char *), dlange_(char *, integer *,
integer *, doublereal *, integer *, doublereal *);
extern /* Subroutine */ int dggqrf_(integer *, integer *, integer *,
doublereal *, integer *, doublereal *, doublereal *, integer *,
doublereal *, doublereal *, integer *, integer *), dlacpy_(char *,
integer *, integer *, doublereal *, integer *, doublereal *,
integer *), dlaset_(char *, integer *, integer *,
doublereal *, doublereal *, doublereal *, integer *);
extern doublereal dlansy_(char *, char *, integer *, doublereal *,
integer *, doublereal *);
extern /* Subroutine */ int dorgqr_(integer *, integer *, integer *,
doublereal *, integer *, doublereal *, doublereal *, integer *,
integer *), dorgrq_(integer *, integer *, integer *, doublereal *,
integer *, doublereal *, doublereal *, integer *, integer *);
static doublereal ulp;
#define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
#define t_ref(a_1,a_2) t[(a_2)*t_dim1 + a_1]
#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
#define af_ref(a_1,a_2) af[(a_2)*af_dim1 + a_1]
#define bf_ref(a_1,a_2) bf[(a_2)*bf_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
September 30, 1994
Purpose
=======
DGQRTS tests DGGQRF, which computes the GQR factorization of an
N-by-M matrix A and a N-by-P matrix B: A = Q*R and B = Q*T*Z.
Arguments
=========
N (input) INTEGER
The number of rows of the matrices A and B. N >= 0.
M (input) INTEGER
The number of columns of the matrix A. M >= 0.
P (input) INTEGER
The number of columns of the matrix B. P >= 0.
A (input) DOUBLE PRECISION array, dimension (LDA,M)
The N-by-M matrix A.
AF (output) DOUBLE PRECISION array, dimension (LDA,N)
Details of the GQR factorization of A and B, as returned
by DGGQRF, see SGGQRF for further details.
Q (output) DOUBLE PRECISION array, dimension (LDA,N)
The M-by-M orthogonal matrix Q.
R (workspace) DOUBLE PRECISION array, dimension (LDA,MAX(M,N))
LDA (input) INTEGER
The leading dimension of the arrays A, AF, R and Q.
LDA >= max(M,N).
TAUA (output) DOUBLE PRECISION array, dimension (min(M,N))
The scalar factors of the elementary reflectors, as returned
by DGGQRF.
B (input) DOUBLE PRECISION array, dimension (LDB,P)
On entry, the N-by-P matrix A.
BF (output) DOUBLE PRECISION array, dimension (LDB,N)
Details of the GQR factorization of A and B, as returned
by DGGQRF, see SGGQRF for further details.
Z (output) DOUBLE PRECISION array, dimension (LDB,P)
The P-by-P orthogonal matrix Z.
//.........这里部分代码省略.........
示例11: glm_gibbs
extern "C" void glm_gibbs(double * rZ, double * rxo, double * rlam, int * rmodelprior, double * rpriorprob, double * rbeta1, double * rbeta2, int * rburnin, int * rniter, int * rscalemixture, double * ralpha, int * rno, int * rna, int * rp, double * B_mcmc, double * prob_mcmc, int * gamma_mcmc, double * phi_mcmc, double * lam_mcmc, double * B_rb, double * prob_rb, double * intercept_mcmc, double * xo_scale)
{
GetRNGstate();
//MCMC Variables//
int burnin=*rburnin;
int niter=*rniter;
//Dimensions//
int p=*rp;
int no=*rno;
int na=*rna;
//Phi Variables//
double phi=1.0;
//Yo Variables//
std::vector<double> Z(rZ, rZ+no);
std::vector<double> xo(rxo, rxo+no*p);
standardize_xo(xo,xo_scale,no,p);
std::vector<double> xoyo(p);
double yobar=0;
std::vector<double> xoxo(p*p);
dgemm_( &transT, &transN, &p, &p, &no, &unity, &*xo.begin(), &no, &*xo.begin(), &no, &inputscale0, &*xoxo.begin(), &p );
//Construct Xa//
std::vector<double> xa(p*(p+1)/2); //Triangular Packed Storage
std::vector<double> d(p);
chol_xa(xa,xoxo,d,p);
//Reserve Memory for Submatrices//
std::vector<double> xog; xog.reserve(no*p);
std::vector<double> xogyo; xogyo.reserve(p);
std::vector<double> xogxog_Lamg; xogxog_Lamg.reserve(p*p);
std::vector<double> xag; xag.reserve(na*p);
//Ya Variables//
std::vector<double> xaya(p);
//Beta Variables//
double intercept=0;
std::vector<double> Bols(p);
std::vector<double> B(p,0.0);
std::vector<double> Bg; Bg.reserve(p);
//Lambda Variables//
int scalemixture=*rscalemixture;
double alpha=*ralpha;
std::vector<double> lam(rlam,rlam+p);
std::vector<double> lamg; lamg.reserve(p); //vector instead of diagonal pxp matrix
//Gamma Variables//
std::vector<int> gamma(p,1);
int p_gamma=std::accumulate(gamma.begin(),gamma.end(),0);
bool gamma_diff=true;
int modelprior=*rmodelprior;
//Probability Variables//
std::vector<double> prob(p);
std::vector<double> odds(p);
std::vector<double> priorprob(rpriorprob,rpriorprob+p);
//Theta Variables//
double theta=0.5;
double beta1=*rbeta1;
double beta2=*rbeta2;
//Store Initial Values//
std::copy(B.begin(),B.end(),B_mcmc);
std::copy(prob.begin(),prob.end(),prob_mcmc);
std::copy(gamma.begin(),gamma.end(),gamma_mcmc);
std::copy(lam.begin(),lam.end(),lam_mcmc);
//Run Gibbs Sampler//
for (int t = 1; t < niter; ++t)
{
//Form Submatrices//
if(p_gamma) submatrices_uncollapsed(gamma_diff,B,xog,xag,lamg,Bg,gamma,lam,xo,xa,p_gamma,p,no,na);
//Draw xoyo//
draw_xoyo(Z,xoyo,yobar,xo,xog,Bg,phi,no,p,p_gamma,intercept);
//Draw xaya//
draw_uncollapsed_xaya(xaya,xa,xag,Bg,phi,na,p,p_gamma);
//Compute Probabilities//
if(modelprior==1)
{
bernoulli_probabilities(prob,odds,Bols,d,xoyo,xaya,priorprob,lam,phi);
}else if(modelprior==2)
{
betabinomial_probabilities(prob,odds,Bols,d,xoyo,xaya,theta,lam,phi);
}else
{
uniform_probabilities(prob,odds,Bols,d,xoyo,xaya,lam,phi);
}
//Draw Gamma//
//.........这里部分代码省略.........
示例12: COMPZ
//.........这里部分代码省略.........
Further Details
===============
Based on contributions by
Jeff Rutter, Computer Science Division, University of California
at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee.
=====================================================================
Test the input parameters.
Parameter adjustments */
/* Table of constant values */
static integer c__2 = 2;
static integer c__9 = 9;
static integer c__0 = 0;
static doublereal c_b18 = 0.;
static doublereal c_b19 = 1.;
static integer c__1 = 1;
/* System generated locals */
integer z_dim1, z_offset, i__1, i__2;
doublereal d__1, d__2;
/* Builtin functions */
double log(doublereal);
integer pow_ii(integer *, integer *);
double sqrt(doublereal);
/* Local variables */
static doublereal tiny;
static integer i__, j, k, m;
static doublereal p;
extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
integer *, doublereal *, doublereal *, integer *, doublereal *,
integer *, doublereal *, doublereal *, integer *);
extern logical lsame_(char *, char *);
extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
doublereal *, integer *);
static integer lwmin;
extern /* Subroutine */ int dlaed0_(integer *, integer *, integer *,
doublereal *, doublereal *, doublereal *, integer *, doublereal *,
integer *, doublereal *, integer *, integer *);
static integer start, ii;
extern doublereal dlamch_(char *);
extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
doublereal *, doublereal *, integer *, integer *, doublereal *,
integer *, integer *), dlacpy_(char *, integer *, integer
*, doublereal *, integer *, doublereal *, integer *),
dlaset_(char *, integer *, integer *, doublereal *, doublereal *,
doublereal *, integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
extern /* Subroutine */ int xerbla_(char *, integer *);
extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);
extern /* Subroutine */ int dsterf_(integer *, doublereal *, doublereal *,
integer *), dlasrt_(char *, integer *, doublereal *, integer *);
static integer liwmin, icompz;
extern /* Subroutine */ int dsteqr_(char *, integer *, doublereal *,
doublereal *, doublereal *, integer *, doublereal *, integer *);
static doublereal orgnrm;
static logical lquery;
static integer smlsiz, dtrtrw, storez, end, lgn;
static doublereal eps;
#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
示例13: dlaqps_
int dlaqps_(int *m, int *n, int *offset, int
*nb, int *kb, double *a, int *lda, int *jpvt,
double *tau, double *vn1, double *vn2, double *auxv,
double *f, int *ldf)
{
/* System generated locals */
int a_dim1, a_offset, f_dim1, f_offset, i__1, i__2;
double d__1, d__2;
/* Builtin functions */
double sqrt(double);
int i_dnnt(double *);
/* Local variables */
int j, k, rk;
double akk;
int pvt;
double temp;
extern double dnrm2_(int *, double *, int *);
double temp2, tol3z;
extern int dgemm_(char *, char *, int *, int *,
int *, double *, double *, int *, double *,
int *, double *, double *, int *),
dgemv_(char *, int *, int *, double *, double *,
int *, double *, int *, double *, double *,
int *);
int itemp;
extern int dswap_(int *, double *, int *,
double *, int *);
extern double dlamch_(char *);
extern int idamax_(int *, double *, int *);
extern int dlarfp_(int *, double *, double *,
int *, double *);
int lsticc, lastrk;
/* -- LAPACK auxiliary routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLAQPS computes a step of QR factorization with column pivoting */
/* of a float 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (KB) */
/* 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) */
//.........这里部分代码省略.........
示例14: main
//.........这里部分代码省略.........
bli_randm( &a );
bli_randm( &b );
bli_randm( &c );
bli_setsc( (2.0/1.0), 0.0, &alpha );
bli_setsc( -(1.0/1.0), 0.0, &beta );
bli_copym( &c, &c_save );
dtime_save = 1.0e9;
for ( r = 0; r < n_repeats; ++r )
{
bli_copym( &c_save, &c );
dtime = bli_clock();
#ifdef PRINT
bli_printm( "a", &a, "%4.1f", "" );
bli_printm( "b", &b, "%4.1f", "" );
bli_printm( "c", &c, "%4.1f", "" );
#endif
#ifdef BLIS
//bli_error_checking_level_set( BLIS_NO_ERROR_CHECKING );
bli_gemm( &alpha,
&a,
&b,
&beta,
&c );
#else
f77_char transa = 'N';
f77_char transb = 'N';
f77_int mm = bli_obj_length( c );
f77_int kk = bli_obj_width_after_trans( a );
f77_int nn = bli_obj_width( c );
f77_int lda = bli_obj_col_stride( a );
f77_int ldb = bli_obj_col_stride( b );
f77_int ldc = bli_obj_col_stride( c );
double* alphap = bli_obj_buffer( alpha );
double* ap = bli_obj_buffer( a );
double* bp = bli_obj_buffer( b );
double* betap = bli_obj_buffer( beta );
double* cp = bli_obj_buffer( c );
dgemm_( &transa,
&transb,
&mm,
&nn,
&kk,
alphap,
ap, &lda,
bp, &ldb,
betap,
cp, &ldc );
#endif
#ifdef PRINT
bli_printm( "c after", &c, "%4.1f", "" );
exit(1);
#endif
dtime_save = bli_clock_min_diff( dtime_save, dtime );
}
gflops = ( 2.0 * m * k * n ) / ( dtime_save * 1.0e9 );
#ifdef BLIS
printf( "data_gemm_blis" );
#else
printf( "data_gemm_%s", BLAS );
#endif
printf( "( %2lu, 1:5 ) = [ %4lu %4lu %4lu %10.3e %6.3f ];\n",
( unsigned long )(p - p_begin + 1)/p_inc + 1,
( unsigned long )m,
( unsigned long )k,
( unsigned long )n, dtime_save, gflops );
bli_obj_free( &alpha );
bli_obj_free( &beta );
bli_obj_free( &a );
bli_obj_free( &b );
bli_obj_free( &c );
bli_obj_free( &c_save );
}
bli_finalize();
return 0;
}示例15: dspr_
/* Subroutine */ int dsbt21_(char *uplo, integer *n, integer *ka, integer *ks,
doublereal *a, integer *lda, doublereal *d__, doublereal *e,
doublereal *u, integer *ldu, doublereal *work, doublereal *result)
{
/* System generated locals */
integer a_dim1, a_offset, u_dim1, u_offset, i__1, i__2, i__3, i__4;
doublereal d__1, d__2;
/* Local variables */
integer j, jc, jr, lw, ika;
doublereal ulp, unfl;
extern /* Subroutine */ int dspr_(char *, integer *, doublereal *,
doublereal *, integer *, doublereal *), dspr2_(char *,
integer *, doublereal *, doublereal *, integer *, doublereal *,
integer *, doublereal *), dgemm_(char *, char *, integer *
, integer *, integer *, doublereal *, doublereal *, integer *,
doublereal *, integer *, doublereal *, doublereal *, integer *);
extern logical lsame_(char *, char *);
doublereal anorm;
char cuplo[1];
logical lower;
doublereal wnorm;
extern doublereal dlamch_(char *), dlange_(char *, integer *,
integer *, doublereal *, integer *, doublereal *),
dlansb_(char *, char *, integer *, integer *, doublereal *,
integer *, doublereal *), dlansp_(char *, char *,
integer *, doublereal *, doublereal *);
/* -- LAPACK test routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DSBT21 generally checks a decomposition of the form */
/* A = U S U' */
/* where ' means transpose, A is symmetric banded, U is */
/* orthogonal, and S is diagonal (if KS=0) or symmetric */
/* tridiagonal (if KS=1). */
/* Specifically: */
/* RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and* */
/* RESULT(2) = | I - UU' | / ( n ulp ) */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER */
/* If UPLO='U', the upper triangle of A and V will be used and */
/* the (strictly) lower triangle will not be referenced. */
/* If UPLO='L', the lower triangle of A and V will be used and */
/* the (strictly) upper triangle will not be referenced. */
/* N (input) INTEGER */
/* The size of the matrix. If it is zero, DSBT21 does nothing. */
/* It must be at least zero. */
/* KA (input) INTEGER */
/* The bandwidth of the matrix A. It must be at least zero. If */
/* it is larger than N-1, then max( 0, N-1 ) will be used. */
/* KS (input) INTEGER */
/* The bandwidth of the matrix S. It may only be zero or one. */
/* If zero, then S is diagonal, and E is not referenced. If */
/* one, then S is symmetric tri-diagonal. */
/* A (input) DOUBLE PRECISION array, dimension (LDA, N) */
/* The original (unfactored) matrix. It is assumed to be */
/* symmetric, and only the upper (UPLO='U') or only the lower */
/* (UPLO='L') will be referenced. */
/* LDA (input) INTEGER */
/* The leading dimension of A. It must be at least 1 */
/* and at least min( KA, N-1 ). */
/* D (input) DOUBLE PRECISION array, dimension (N) */
/* The diagonal of the (symmetric tri-) diagonal matrix S. */
/* E (input) DOUBLE PRECISION array, dimension (N-1) */
/* The off-diagonal of the (symmetric tri-) diagonal matrix S. */
/* E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and */
/* (3,2) element, etc. */
/* Not referenced if KS=0. */
/* U (input) DOUBLE PRECISION array, dimension (LDU, N) */
/* The orthogonal matrix in the decomposition, expressed as a */
/* dense matrix (i.e., not as a product of Householder */
/* transformations, Givens transformations, etc.) */
/* LDU (input) INTEGER */
//.........这里部分代码省略.........