dnlaso.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,804 行 · 第 1/5 页
C
1,804 行
#include "f2c.h"
#include "netlib.h"
extern double sqrt(double); /* #include <math.h> */
static void dlabax_(const integer *n, const integer *nband, doublereal *a, doublereal *x, doublereal *y);
static void dlabcm_(const integer *n, const integer *nband, const integer *nl, const integer *nr,
doublereal *a, doublereal *eigval, const integer *lde, doublereal *eigvec,
doublereal *atol, doublereal *artol, doublereal *bound, doublereal *atemp, doublereal *d, doublereal *vtemp);
static void dlabfc_(const integer *n, const integer *nband, doublereal *a, doublereal *sigma, const integer *number,
const integer *lde, doublereal *eigvec, integer *numl, integer *ldad,
doublereal *atemp, doublereal *d, doublereal *atol);
static void dlaeig_(const integer *n, const integer *nband, const integer *nl, const integer *nr,
doublereal *a, doublereal *eigval, const integer *lde,
doublereal *eigvec, doublereal *bound, doublereal *atemp, doublereal *d,
doublereal *vtemp, doublereal *eps, doublereal *tmin, doublereal *tmax);
static void dlager_(const integer *n, const integer *nband, const integer *nstart,
doublereal *a, doublereal *tmin, doublereal *tmax);
static void dlaran_(const integer *n, doublereal *x);
static void dmvpc_(const integer *nblock, const doublereal *bet, const integer *maxj, const integer *j,
const doublereal *s, const integer *number, doublereal *resnrm, doublereal *orthcf, doublereal *rv);
static void dnppla_(void (*op)(const integer*,const integer*,const doublereal*,doublereal*),
void (*iovect)(const integer*,const integer*,doublereal*,const integer*,const integer*),
const integer *n, const integer *nperm, integer *nop, const integer *nmval,
doublereal *val, const integer *nmvec, doublereal *vec, const integer *nblock,
doublereal *h, doublereal *hv, doublereal *p, doublereal *q, doublereal *bound,
doublereal *d, doublereal *delta, logical *small, logical *raritz, doublereal *eps);
static void dnwla_(void (*op)(const integer*,const integer*,const doublereal*,doublereal*),
void (*iovect)(const integer*,const integer*,doublereal*,const integer*,const integer*),
const integer *n, const integer *nband, const integer *nval,
const integer *nfig, integer *nperm, doublereal *val, const integer *nmvec, doublereal *vec,
const integer *nblock, const integer *maxop, const integer *maxj, integer *nop,
doublereal *p1, doublereal *p0, doublereal *res, doublereal *tau, doublereal *otau,
doublereal *t, doublereal *alp, doublereal *bet, doublereal *s, doublereal *p2,
doublereal *bound, doublereal *atemp, doublereal *vtemp,
doublereal *d, integer *ind, logical *small, logical *raritz,
doublereal *delta, doublereal *eps, integer *ierr);
static void dortqr_(const integer *nz, const integer *n, const integer *nblock, doublereal *z, doublereal *b);
static void dvsort_(const integer *num, doublereal *val, doublereal *res, const integer *iflag,
doublereal *v, const integer *nmvec, const integer *n, doublereal *vec);
/* Table of constant values */
static integer c__0 = 0;
static integer c__1 = 1;
static doublereal c__10 = 0.1;
static doublereal c__00 = 0.0;
/* VERSION 2 DOES NOT USE EISPACK */
/* ------------------------------------------------------------------ */
/* Subroutine */ void dnlaso_(op, iovect, n, nval, nfig, nperm, nmval, val, nmvec, vec, nblock, maxop, maxj, work, ind, ierr)
void (*op) (const integer* n,const integer* m,const doublereal* p,doublereal* q);
void (*iovect) (const integer* n,const integer* m,doublereal* q,const integer* j,const integer* k);
const integer *n, *nval, *nfig, *nmval;
integer *nperm;
doublereal *val;
const integer *nmvec;
doublereal *vec;
const integer *nblock, *maxop, *maxj;
doublereal *work;
integer *ind, *ierr;
{
/* Local variables */
static doublereal temp, tarr;
static integer i, m, nband;
static doublereal delta;
static logical small;
static integer i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, nv;
static logical raritz;
static doublereal eps;
static integer nop;
/* AUTHOR/IMPLEMENTER D.S.SCOTT-B.N.PARLETT/D.S.SCOTT */
/* */
/* COMPUTER SCIENCES DEPARTMENT */
/* UNIVERSITY OF TEXAS AT AUSTIN */
/* AUSTIN, TX 78712 */
/* */
/* VERSION 2 ORIGINATED APRIL 1982 */
/* */
/* CURRENT VERSION JUNE 1983 */
/* DNLASO FINDS A FEW EIGENVALUES AND EIGENVECTORS AT EITHER END OF */
/* THE SPECTRUM OF A LARGE SPARSE SYMMETRIC MATRIX. THE SUBROUTINE */
/* DNLASO IS PRIMARILY A DRIVER FOR SUBROUTINE DNWLA WHICH IMPLEMENTS */
/* THE LANCZOS ALGORITHM WITH SELECTIVE ORTHOGONALIZATION AND */
/* SUBROUTINE DNPPLA WHICH POST PROCESSES THE OUTPUT OF DNWLA. */
/* HOWEVER DNLASO DOES CHECK FOR INCONSISTENCIES IN THE CALLING */
/* PARAMETERS AND DOES PREPROCESS ANY USER SUPPLIED EIGENPAIRS. */
/* DNLASO ALWAYS LOOKS FOR THE SMALLEST (LEFTMOST) EIGENVALUES. IF */
/* THE LARGEST EIGENVALUES ARE DESIRED DNLASO IMPLICITLY USES THE */
/* NEGATIVE OF THE MATRIX. */
/* */
/* ON INPUT */
/* */
/* OP A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE */
/* OP(N,M,P,Q). P AND Q ARE N X M MATRICES AND Q IS */
/* RETURNED AS THE MATRIX TIMES P. */
/* */
/* IOVECT A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE */
/* IOVECT(N,M,Q,J,K). Q IS AN N X M MATRIX. IF K = 0 */
/* THE COLUMNS OF Q ARE STORED AS THE (J-M+1)TH THROUGH */
/* THE JTH LANCZOS VECTORS. IF K = 1 THEN Q IS RETURNED */
/* AS THE (J-M+1)TH THROUGH THE JTH LANCZOS VECTORS. SEE */
/* DOCUMENTATION FOR FURTHER DETAILS AND EXAMPLES. */
/* */
/* N THE ORDER OF THE MATRIX. */
/* */
/* NVAL NVAL SPECIFIES THE EIGENVALUES TO BE FOUND. */
/* DABS(NVAL) IS THE NUMBER OF EIGENVALUES DESIRED. */
/* IF NVAL < 0 THE ALGEBRAICALLY SMALLEST (LEFTMOST) */
/* EIGENVALUES ARE FOUND. IF NVAL > 0 THE ALGEBRAICALLY */
/* LARGEST (RIGHTMOST) EIGENVALUES ARE FOUND. NVAL MUST NOT */
/* BE ZERO. DABS(NVAL) MUST BE LESS THAN MAXJ/2. */
/* */
/* NFIG THE NUMBER OF DECIMAL DIGITS OF ACCURACY DESIRED IN THE */
/* EIGENVALUES. NFIG MUST BE GREATER THAN OR EQUAL TO 1. */
/* */
/* NPERM AN INTEGER VARIABLE WHICH SPECIFIES THE NUMBER OF USER */
/* SUPPLIED EIGENPAIRS. IN MOST CASES NPERM WILL BE ZERO. SEE */
/* DOCUMENTAION FOR FURTHER DETAILS OF USING NPERM GREATER */
/* THAN ZERO. NPERM MUST NOT BE LESS THAN ZERO. */
/* */
/* NMVAL THE ROW DIMENSION OF THE ARRAY VAL. NMVAL MUST BE GREATER */
/* THAN OR EQUAL TO DABS(NVAL). */
/* */
/* VAL A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW */
/* DIMENSION NMVAL AND COLUMN DIMENSION AT LEAST 4. IF NPERM */
/* IS GREATER THAN ZERO THEN CERTAIN INFORMATION MUST BE STORED */
/* IN VAL. SEE DOCUMENTATION FOR DETAILS. */
/* */
/* NMVEC THE ROW DIMENSION OF THE ARRAY VEC. NMVEC MUST BE GREATER */
/* THAN OR EQUAL TO N. */
/* */
/* VEC A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW */
/* DIMENSION NMVEC AND COLUMN DIMENSION AT LEAST DABS(NVAL). IF */
/* NPERM > 0 THEN THE FIRST NPERM COLUMNS OF VEC MUST */
/* CONTAIN THE USER SUPPLIED EIGENVECTORS. */
/* */
/* NBLOCK THE BLOCK SIZE. SEE DOCUMENTATION FOR CHOOSING */
/* AN APPROPRIATE VALUE FOR NBLOCK. NBLOCK MUST BE GREATER */
/* THAN ZERO AND LESS THAN MAXJ/6. */
/* */
/* MAXOP AN UPPER BOUND ON THE NUMBER OF CALLS TO THE SUBROUTINE */
/* OP. DNLASO TERMINATES WHEN MAXOP IS EXCEEDED. SEE */
/* DOCUMENTATION FOR GUIDELINES IN CHOOSING A VALUE FOR MAXOP. */
/* */
/* MAXJ AN INDICATION OF THE AVAILABLE STORAGE (SEE WORK AND */
/* DOCUMENTATION ON IOVECT). FOR THE FASTEST CONVERGENCE MAXJ */
/* SHOULD BE AS LARGE AS POSSIBLE, ALTHOUGH IT IS USELESS TO HAVE */
/* MAXJ LARGER THAN MAXOP*NBLOCK. */
/* */
/* WORK A DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST AS */
/* LARGE AS */
/* */
/* 2*N*NBLOCK + MAXJ*(NBLOCK+NV+2) + 2*NBLOCK*NBLOCK + 3*NV */
/* */
/* + THE MAXIMUM OF */
/* N*NBLOCK */
/* AND */
/* MAXJ*(2*NBLOCK+3) + 2*NV + 6 + (2*NBLOCK+2)*(NBLOCK+1) */
/* */
/* WHERE NV = DABS(NVAL) */
/* */
/* THE FIRST N*NBLOCK ELEMENTS OF WORK MUST CONTAIN THE DESIRED */
/* STARTING VECTORS. SEE DOCUMENTATION FOR GUIDELINES IN */
/* CHOOSING STARTING VECTORS. */
/* */
/* IND AN INTEGER ARRAY OF DIMENSION AT LEAST DABS(NVAL). */
/* */
/* IERR AN INTEGER VARIABLE. */
/* */
/* ON OUTPUT */
/* */
/* NPERM THE NUMBER OF EIGENPAIRS NOW KNOWN. */
/* */
/* VEC THE FIRST NPERM COLUMNS OF VEC CONTAIN THE EIGENVECTORS. */
/* */
/* VAL THE FIRST COLUMN OF VAL CONTAINS THE CORRESPONDING */
/* EIGENVALUES. THE SECOND COLUMN CONTAINS THE RESIDUAL NORMS OF */
/* THE EIGENPAIRS WHICH ARE BOUNDS ON THE ACCURACY OF THE EIGEN- */
/* VALUES. THE THIRD COLUMN CONTAINS MORE REALISTIC ESTIMATES */
/* OF THE ACCURACY OF THE EIGENVALUES. THE FOURTH COLUMN CONTAINS */
/* ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS. SEE */
/* DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES. */
/* */
/* WORK IF WORK IS TERMINATED BEFORE COMPLETION (IERR = -2) */
/* THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS */
/* FOR RESTARTING THE ALGORITHM AND DNLASO CAN BE IMMEDIATELY */
/* RECALLED TO CONTINUE WORKING ON THE PROBLEM. */
/* */
/* IND IND(1) CONTAINS THE ACTUAL NUMBER OF CALLS TO OP. ON SOME */
/* OCCASIONS THE NUMBER OF CALLS TO OP MAY BE SLIGHTLY LARGER */
/* THAN MAXOP. */
/* */
/* IERR AN ERROR COMPLETION CODE. THE NORMAL COMPLETION CODE IS */
/* ZERO. SEE THE DOCUMENTATION FOR INTERPRETATIONS OF NON-ZERO */
/* COMPLETION CODES. */
/* */
/* INTERNAL VARIABLES. */
/* */
/* NOP RETURNED FROM DNWLA AS THE NUMBER OF CALLS TO THE */
/* SUBROUTINE OP. */
/* */
/* NV SET EQUAL TO DABS(NVAL), THE NUMBER OF EIGENVALUES DESIRED, */
/* AND PASSED TO DNWLA. */
/* */
/* SMALL SET TO .TRUE. IF THE SMALLEST EIGENVALUES ARE DESIRED. */
/* */
/* RARITZ RETURNED FROM DNWLA AND PASSED TO DNPPLA. RARITZ IS .TRUE. */
/* IF A FINAL RAYLEIGH-RITZ PROCEDURE IS NEEDED. */
/* */
/* DELTA RETURNED FROM DNWLA AS THE EIGENVALUE OF THE MATRIX */
/* WHICH IS CLOSEST TO THE DESIRED EIGENVALUES. */
/* */
/* DNPPLA A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED */
/* BY DNWLA. */
/* */
/* DNWLA A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM */
/* WITH SELECTIVE ORTHOGONALIZATION. */
/* */
/* DMVPC A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND */
/* ORTHOGONALITY COEFFICIENT OF GIVEN RITZ VECTORS. */
/* */
/* DORTQR A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS */
/* USING HOUSEHOLDER REFLECTIONS. */
/* */
/* DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP A SUBSET OF THE BASIC LINEAR */
/* ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION. */
/* */
/* DLARAN A SUBROUTINE TO GENERATE RANDOM VECTORS */
/* */
/* DLAEIG, DLAGER, DLABCM, DLABFC SUBROUTINES FOR BAND EIGENVALUE */
/* CALCULATIONS. */
/* */
/* ------------------------------------------------------------------ */
/* THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS. */
nv = abs(*nval);
ind[0] = 0;
*ierr = 0;
if (*n < *nblock * 6) {
*ierr = 1;
}
if (*nfig <= 0) {
*ierr += 2;
}
if (*nmvec < *n) {
*ierr += 4;
}
if (*nperm < 0) {
*ierr += 8;
}
if (*maxj < *nblock * 6) {
*ierr += 16;
}
if (nv < max(1,*nperm)) {
*ierr += 32;
}
if (nv > *nmval) {
*ierr += 64;
}
if (nv > *maxop) {
*ierr += 128;
}
if (nv >= *maxj / 2) {
*ierr += 256;
}
if (*nblock < 1) {
*ierr += 512;
}
if (*ierr != 0) {
return;
}
small = *nval < 0;
/* ------------------------------------------------------------------ */
/* THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS. */
/* IF A USER SUPPLIED VECTOR IS ZERO OR IF SIGNIFICANT CANCELLATION */
/* OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO -1 */
/* AND DNLASO TERMINATES. */
if (*nperm == 0) {
goto L110;
}
/* THIS NEGATES THE USER SUPPLIED EIGENVALUES WHEN THE LARGEST */
/* EIGENVALUES ARE DESIRED, SINCE DNWLA WILL IMPLICITLY USE THE */
/* NEGATIVE OF THE MATRIX. */
if (!small)
for (i = 0; i < *nperm; ++i) {
val[i] = -val[i];
}
/* THIS SORTS THE USER SUPPLIED VALUES AND VECTORS. */
dvsort_(nperm, val, &val[*nmval], &c__0, &tarr, nmvec, n, vec);
/* THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON. */
/* IT ALSO INSURES THAT THE RESIDUAL NORMS ARE POSITIVE. */
for (i = 0; i < *nperm; ++i) {
val[i + *nmval] = abs(val[i + *nmval]);
val[i + *nmval * 2] = dnrm2_(n, &vec[i * *nmvec], &c__1);
}
/* THIS PERFORMS THE ORTHONORMALIZATION. */
m = *n * *nblock;
dortqr_(nmvec, n, nperm, vec, &work[m]);
for (i = 0; i < *nperm; ++i, m += *nperm + 1) {
if (abs(work[m]) <= val[i + *nmval * 2] * .9) {
*ierr = -1;
return;
}
}
/* THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN */
/* THE ARRAY WORK FOR LATER REFERENCE IN DNWLA. */
m = (*n << 1) * *nblock;
dcopy_(nperm, &val[*nmval], &c__1, &work[m], &c__1);
/* THIS SETS EPS TO AN APPROXIMATION OF THE RELATIVE MACHINE */
/* PRECISION */
/* ***THIS SHOULD BE REPLACED BY AN ASSIGNMENT STATEMENT */
/* ***IN A PRODUCTION CODE */
L110:
eps = 1.;
for (i = 0; i < 1000; ++i) {
eps *= .5;
temp = eps + 1.;
if (temp == 1.) {
break;
}
}
/* ------------------------------------------------------------------ */
/* THIS SECTION CALLS DNWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM */
/* WITH SELECTIVE ORTHOGONALIZATION. */
nband = *nblock + 1;
i1 = *n * *nblock;
i2 = i1 + *n * *nblock;
i3 = i2 + nv;
i4 = i3 + nv;
i5 = i4 + nv;
i6 = i5 + *maxj * nband;
i7 = i6 + *nblock * *nblock;
i8 = i7 + *nblock * *nblock;
i9 = i8 + *maxj * (nv + 1);
i10 = i9 + *nblock;
i11 = i10 + (nv << 1) + 6;
i12 = i11 + *maxj * ((*nblock << 1) + 1);
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?