dnlaso.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,804 行 · 第 1/5 页
C
1,804 行
i13 = i12 + *maxj;
dnwla_(op, iovect, n, &nband, &nv, nfig, nperm, val, nmvec, vec, nblock,
maxop, maxj, &nop, work, &work[i1], &work[i2], &work[i3], &work[i4],
&work[i5], &work[i6], &work[i7], &work[i8], &work[i9], &work[i10],
&work[i11], &work[i12], &work[i13], ind, &small, &raritz, &delta, &eps, ierr);
/* ------------------------------------------------------------------ */
/* THIS SECTION CALLS DNPPLA (THE POST PROCESSOR). */
if (*nperm == 0) {
ind[0] = nop;
return;
}
i1 = *n * *nblock;
i2 = i1 + *nperm * *nperm;
i3 = i2 + *nperm * *nperm;
i4 = i3 + max(*n * *nblock, 2 * *nperm * *nperm);
i5 = i4 + *n * *nblock;
i6 = i5 + (*nperm << 1) + 4;
dnppla_(op, iovect, n, nperm, &nop, nmval, val, nmvec, vec, nblock,
&work[i1], &work[i2], &work[i3], &work[i4], &work[i5], &work[i6],
&delta, &small, &raritz, &eps);
ind[0] = nop;
return;
} /* dnlaso_ */
/* ------------------------------------------------------------------ */
/* Subroutine */ void dnwla_(op, iovect, n, nband, nval, nfig, nperm, val,
nmvec, vec, nblock, maxop, maxj, nop, p1, p0, res, tau, otau, t, alp,
bet, s, p2, bound, atemp, vtemp, d, ind, small, raritz, delta, eps, ierr)
/* Subroutine */ void (*op) (const integer*,const integer*,const doublereal*,doublereal*);
/* Subroutine */ void (*iovect) (const integer*,const integer*,doublereal*,const integer*,const integer*);
const integer *n, *nband, *nval, *nfig;
integer *nperm;
doublereal *val;
const integer *nmvec;
doublereal *vec;
const integer *nblock, *maxop, *maxj;
integer *nop;
doublereal *p1, *p0, *res, *tau, *otau, *t, *alp, *bet, *s, *p2, *bound, *atemp, *vtemp, *d;
integer *ind;
logical *small, *raritz;
doublereal *delta, *eps;
integer *ierr;
{
/* System generated locals */
integer i__1;
doublereal d__1;
/* Local variables */
static doublereal tola, temp, tolg, tmin, tmax, tarr;
static logical test;
static doublereal zero=0., utol;
static integer i, j, k, l, m;
static integer ngood, nleft;
static doublereal anorm;
static integer mtemp;
static integer i1;
static doublereal pnorm, epsrt, rnorm;
static integer ng;
static doublereal betmin, alpmin, betmax, alpmax;
static integer ntheta;
static logical enough;
static integer number, nstart;
/* DNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE */
/* ORTHOGONALIZATION. */
/* */
/* NBAND NBLOCK + 1 THE BAND WIDTH OF T. */
/* */
/* NVAL THE NUMBER OF DESIRED EIGENVALUES. */
/* */
/* NPERM THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS */
/* INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE */
/* ALGORITHM IS ITERATED). PERMANENT VECTORS ARE ORTHOGONAL */
/* TO THE CURRENT KRYLOV SUBSPACE. */
/* */
/* NOP THE NUMBER OF CALLS TO OP. */
/* */
/* P0, P1, AND P2 THE CURRENT BLOCKS OF LANCZOS VECTORS. */
/* */
/* RES THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS. */
/* */
/* TAU AND OTAU USED TO MONITOR THE NEED FOR ORTHOGONALIZATION. */
/* */
/* T THE BAND MATRIX. */
/* */
/* ALP THE CURRENT DIAGONAL BLOCK. */
/* */
/* BET THE CURRENT OFF DIAGONAL BLOCK. */
/* */
/* BOUND, ATEMP, VTEMP, D */
/* TEMPORARY STORAGE USED BY THE BAND EIGENVALUE SOLVER DLAEIG. */
/* */
/* S EIGENVECTORS OF T. */
/* */
/* SMALL .TRUE. IF THE SMALL EIGENVALUES ARE DESIRED. */
/* */
/* RARITZ RETURNED AS .TRUE. IF A FINAL RAYLEIGH-RITZ PROCEDURE */
/* IS TO BE DONE. */
/* */
/* DELTA RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE */
/* OF THE MATRIX. USED IN ESTIMATING THE ACCURACY OF THE */
/* COMPUTED EIGENVALUES. */
/* */
/* INTERNAL VARIABLES USED */
/* */
/* J THE CURRENT DIMENSION OF T. (THE DIMENSION OF THE CURRENT */
/* KRYLOV SUBSPACE. */
/* */
/* NGOOD THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS */
/* LIE IN THE CURRENT KRYLOV SUBSPACE). */
/* */
/* NLEFT THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED, */
/* I.E. NLEFT = NVAL - NPERM. */
/* */
/* NUMBER = NPERM + NGOOD. */
/* */
/* ANORM AN ESTIMATE OF THE NORM OF THE MATRIX. */
/* */
/* EPS THE RELATIVE MACHINE PRECISION. */
/* */
/* UTOL THE USER TOLERANCE. */
/* */
/* TARR AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN */
/* CALLS TO DLAEIG */
/* */
/* DZERO AN ARRAY OF LENGTH ONE CONTAINING DZERO, USED TO INSURE TYPE */
/* CONSISTENCY IN CALLS TO DCOPY */
rnorm = 0.;
if (*nperm != 0) {
rnorm = max(-val[0],val[*nperm-1]);
}
pnorm = rnorm;
*delta = 1e31;
epsrt = sqrt(*eps);
nleft = *nval - *nperm;
*nop = 0;
number = *nperm;
*raritz = FALSE_;
utol = max((*n) * *eps, pow_di(&c__10, nfig));
j = *maxj;
/* ------------------------------------------------------------------ */
/* ANY ITERATION OF THE ALGORITHM BEGINS HERE. */
L30:
for (i = 0; i < *nblock; ++i) {
temp = dnrm2_(n, &p1[i * *n], &c__1);
if (temp == 0.) {
dlaran_(n, &p1[i * *n]);
}
}
for (i = 0; i < *nperm; ++i) {
tau[i] = 1.;
otau[i] = 0.;
}
i__1 = *n * *nblock;
dcopy_(&i__1, &zero, &c__0, p0, &c__1);
i__1 = *nblock * *nblock;
dcopy_(&i__1, &zero, &c__0, bet, &c__1);
i__1 = j * *nband;
dcopy_(&i__1, &zero, &c__0, t, &c__1);
mtemp = *nval + 1;
for (i = 0; i < mtemp; ++i) {
dcopy_(&j, &zero, &c__0, &s[i * *maxj], &c__1);
}
ngood = 0;
tmin = 1e30;
tmax = -1e30;
test = TRUE_;
enough = FALSE_;
betmax = 0.;
j = 0;
/* ------------------------------------------------------------------ */
/* THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP. */
L80:
j += *nblock;
/* THIS IS THE SELECTIVE ORTHOGONALIZATION. */
for (i = 0; i < number; ++i) {
if (tau[i] < epsrt) {
continue;
}
test = TRUE_;
tau[i] = 0.;
if (otau[i] != 0.) {
otau[i] = 1.;
}
for (k = 0; k < *nblock; ++k) {
temp = -ddot_(n, &vec[i * *nmvec], &c__1, &p1[k * *n], &c__1);
daxpy_(n, &temp, &vec[i * *nmvec], &c__1, &p1[k * *n], &c__1);
/* THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A */
/* NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR. THE ALGORITHM IS */
/* TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST. */
if (abs(temp * bet[k + k * *nblock]) > (*n) * epsrt * anorm && i >= *nperm) {
goto L380;
}
}
}
/* IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET. */
if (test)
dortqr_(n, n, nblock, p1, alp);
if (test && j != *nblock)
for (i = 0; i < *nblock; ++i) {
if (alp[i + i * *nblock] > 0.) {
continue;
}
m = j - (*nblock << 1) + i;
l = *nblock;
for (k = i; k < *nblock; ++k, --l, ++m) {
bet[i + k * *nblock] = -bet[i + k * *nblock];
t[l + m * *nband] = -t[l + m * *nband];
}
}
test = FALSE_;
/* THIS IS THE LANCZOS STEP. */
(*op)(n, nblock, p1, p2);
++(*nop);
(*iovect)(n, nblock, p1, &j, &c__0);
/* THIS COMPUTES P2=P2-P0*BET(TRANSPOSE) */
for (i = 0; i < *nblock; ++i) {
for (k = i; k < *nblock; ++k) {
d__1 = -bet[i + k * *nblock];
daxpy_(n, &d__1, &p0[k * *n], &c__1, &p2[i * *n], &c__1);
}
}
/* THIS COMPUTES ALP AND P2=P2-P1*ALP. */
for (i = 0; i < *nblock; ++i) {
for (k = 0; k <= i; ++k) {
i1 = i - k;
alp[i1 + k * *nblock] = ddot_(n, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
d__1 = -alp[i1 + k * *nblock];
daxpy_(n, &d__1, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
if (k != i) {
d__1 = -alp[i1 + k * *nblock];
daxpy_(n, &d__1, &p1[k * *n], &c__1, &p2[i * *n], &c__1);
}
}
}
/* REORTHOGONALIZATION OF THE SECOND BLOCK */
if (j == *nblock)
for (i = 0; i < *nblock; ++i) {
for (k = 0; k <= i; ++k) {
temp = -ddot_(n, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
daxpy_(n, &temp, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
if (k != i) {
daxpy_(n, &temp, &p1[k * *n], &c__1, &p2[i * *n], &c__1);
}
i1 = i - k;
alp[i1 + k * *nblock] += temp;
}
}
/* THIS ORTHONORMALIZES THE NEXT BLOCK */
dortqr_(n, n, nblock, p2, bet);
/* THIS STORES ALP AND BET IN T. */
for (i = 0; i < *nblock; ++i) {
m = j - *nblock + i;
for (k = i; k < *nblock; ++k) {
l = k - i;
t[l + m * *nband] = alp[l + i * *nblock];
}
for (k = 0; k <= i; ++k) {
l = *nblock - i + k;
t[l + m * *nband] = bet[k + i * *nblock];
}
}
/* THIS NEGATES T IF SMALL IS FALSE. */
if (! *small)
for (i = j - *nblock; i < j; ++i) {
for (k = 0; k <= l; ++k) { /* FIXME *** This must be an error! (already in the fortran code) -- l is undefined *** */
t[k + i * *nband] = -t[k + i * *nband];
}
}
/* THIS SHIFTS THE LANCZOS VECTORS */
i__1 = *nblock * *n;
dcopy_(&i__1, p1, &c__1, p0, &c__1);
dcopy_(&i__1, p2, &c__1, p1, &c__1);
i__1 = j - *nblock + 1;
dlager_(&j, nband, &i__1, t, &tmin, &tmax);
anorm = max(max(rnorm,tmax),-tmin);
/* THIS COMPUTES THE EXTREME EIGENVALUES OF ALP. */
if (number != 0) {
dcopy_(nblock, &zero, &c__0, p2, &c__1);
dlaeig_(nblock, nblock, &c__1, &c__1, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &tmin, &tmax);
alpmin = tarr;
dcopy_(nblock, &zero, &c__0, p2, &c__1);
dlaeig_(nblock, nblock, nblock, nblock, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &tmin, &tmax);
alpmax = tarr;
}
/* THIS COMPUTES ALP = BET(TRANSPOSE)*BET. */
for (i = 0; i < *nblock; ++i) {
for (k = 0; k <= i; ++k) {
l = i - k;
i__1 = *nblock - i;
alp[l + k * *nblock] = ddot_(&i__1, &bet[i + i * *nblock], nblock, &bet[k + i * *nblock], nblock);
}
}
if (number == 0) {
goto L330;
}
/* THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET. */
dcopy_(nblock, &zero, &c__0, p2, &c__1);
d__1 = anorm * anorm;
dlaeig_(nblock, nblock, &c__1, &c__1, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &c__00, &d__1);
betmin = sqrt(tarr);
/* THIS UPDATES TAU AND OTAU. */
for (i = 0; i < number; ++i) {
temp = (tau[i] * max(alpmax-val[i],val[i]-alpmin) + otau[i] * betmax + *eps * anorm) / betmin;
if (i < *nperm) {
temp += res[i] / betmin;
}
otau[i] = tau[i];
tau[i] = temp;
}
/* THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET. */
L330:
dcopy_(nblock, &zero, &c__0, p2, &c__1);
d__1 = anorm * anorm;
dlaeig_(nblock, nblock, nblock, nblock, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &c__00, &d__1);
betmax = sqrt(tarr);
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?