snlaso.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,841 行 · 第 1/5 页
C
1,841 行
#include "f2c.h"
#include "netlib.h"
extern double sqrt(double); /* #include <math.h> */
static void slabax_(const integer *n, const integer *nband, real *a, real *x, real *y);
static void slabcm_(const integer *n, const integer *nband, const integer *nl, const integer *nr,
real *a, real *eigval, const integer *lde, real *eigvec,
real *atol, real *artol, real *bound, real *atemp, real *d, real *vtemp);
static void slabfc_(const integer *n, const integer *nband, real *a, real *sigma, const integer *number,
const integer *lde, real *eigvec, integer *numl, integer *ldad,
real *atemp, real *d, real *atol);
static void slaeig_(const integer *n, const integer *nband, const integer *nl, const integer *nr,
real *a, real *eigval, const integer *lde,
real *eigvec, real *bound, real *atemp, real *d,
real *vtemp, real *eps, real *tmin, real *tmax);
static void slager_(const integer *n, const integer *nband, const integer *nstart,
real *a, real *tmin, real *tmax);
static void slaran_(const integer *n, real *x);
static void smvpc_(const integer *nblock, const real *bet, const integer *maxj, const integer *j,
const real *s, const integer *number, real *resnrm, real *orthcf, real *rv);
static void snppla_(void (*op)(const integer*,const integer*,const real*,real*),
void (*iovect)(const integer*,const integer*,real*,const integer*,const integer*),
const integer *n, const integer *nperm, integer *nop, const integer *nmval,
real *val, const integer *nmvec, real *vec, const integer *nblock,
real *h, real *hv, real *p, real *q, real *bound,
real *d, real *delta, logical *small, logical *raritz, real *eps);
static void snwla_(void (*op)(const integer*,const integer*,const real*,real*),
void (*iovect)(const integer*,const integer*,real*,const integer*,const integer*),
const integer *n, const integer *nband, const integer *nval,
const integer *nfig, integer *nperm, real *val, const integer *nmvec, real *vec,
const integer *nblock, const integer *maxop, const integer *maxj, integer *nop,
real *p1, real *p0, real *res, real *tau, real *otau,
real *t, real *alp, real *bet, real *s, real *p2,
real *bound, real *atemp, real *vtemp,
real *d, integer *ind, logical *small, logical *raritz,
real *delta, real *eps, integer *ierr);
static void sortqr_(const integer *nz, const integer *n, const integer *nblock, real *z, real *b);
static void svsort_(const integer *num, real *val, real *res, const integer *iflag,
real *v, const integer *nmvec, const integer *n, real *vec);
/* Table of constant values */
static integer c__0 = 0;
static integer c__1 = 1;
static real c__10 = 0.1f;
static real c__00 = 0.0f;
/* VERSION 2 DOES NOT USE EISPACK */
/* ------------------------------------------------------------------ */
/* Subroutine */ void snlaso_(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 real* p,real* q);
void (*iovect) (const integer* n,const integer* m,real* q,const integer* j,const integer* k);
const integer *n, *nval, *nfig, *nmval;
integer *nperm;
real *val;
const integer *nmvec;
real *vec;
const integer *nblock, *maxop, *maxj;
real *work;
integer *ind, *ierr;
{
/* Local variables */
static real temp, tarr;
static integer i, m, nband;
static real delta;
static logical small;
static integer i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, nv;
static logical raritz;
static real 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 */
/* SNLASO FINDS A FEW EIGENVALUES AND EIGENVECTORS AT EITHER END OF */
/* THE SPECTRUM OF A LARGE SPARSE SYMMETRIC MATRIX. THE SUBROUTINE */
/* SNLASO IS PRIMARILY A DRIVER FOR SUBROUTINE SNWLA WHICH IMPLEMENTS */
/* THE LANCZOS ALGORITHM WITH SELECTIVE ORTHOGONALIZATION AND */
/* SUBROUTINE SNPPLA WHICH POST PROCESSES THE OUTPUT OF SNWLA. */
/* HOWEVER SNLASO DOES CHECK FOR INCONSISTENCIES IN THE CALLING */
/* PARAMETERS AND DOES PREPROCESS ANY USER SUPPLIED EIGENPAIRS. */
/* SNLASO ALWAYS LOOKS FOR THE SMALLEST (LEFTMOST) EIGENVALUES. IF */
/* THE LARGEST EIGENVALUES ARE DESIRED SNLASO 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. */
/* ABS(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. ABS(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 ABS(NVAL). */
/* */
/* VAL A TWO DIMENSIONAL REAL 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 REAL ARRAY OF ROW */
/* DIMENSION NMVEC AND COLUMN DIMENSION AT LEAST ABS(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. SNLASO 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 REAL 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 = ABS(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 ABS(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 SNLASO 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 SNWLA AS THE NUMBER OF CALLS TO THE */
/* SUBROUTINE OP. */
/* */
/* NV SET EQUAL TO ABS(NVAL), THE NUMBER OF EIGENVALUES DESIRED, */
/* AND PASSED TO SNWLA. */
/* */
/* SMALL SET TO .TRUE. IF THE SMALLEST EIGENVALUES ARE DESIRED. */
/* */
/* RARITZ RETURNED FROM SNWLA AND PASSED TO SNPPLA. RARITZ IS .TRUE. */
/* IF A FINAL RAYLEIGH-RITZ PROCEDURE IS NEEDED. */
/* */
/* DELTA RETURNED FROM SNWLA AS THE EIGENVALUE OF THE MATRIX */
/* WHICH IS CLOSEST TO THE DESIRED EIGENVALUES. */
/* */
/* SNPPLA A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED */
/* BY SNWLA. */
/* */
/* SNWLA A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM */
/* WITH SELECTIVE ORTHOGONALIZATION. */
/* */
/* SMVPC A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND */
/* ORTHOGONALITY COEFFICIENT OF GIVEN RITZ VECTORS. */
/* */
/* SORTQR A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS */
/* USING HOUSEHOLDER REFLECTIONS. */
/* */
/* SAXPY,SCOPY,SDOT,SNRM2,SSCAL,SSWAP A SUBSET OF THE BASIC LINEAR */
/* ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION. */
/* */
/* SLARAN A SUBROUTINE TO GENERATE RANDOM VECTORS */
/* */
/* SLAEIG, SLAGER, SLABCM, SLABFC 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 SNLASO TERMINATES. */
if (*nperm == 0) {
goto L110;
}
/* THIS NEGATES THE USER SUPPLIED EIGENVALUES WHEN THE LARGEST */
/* EIGENVALUES ARE DESIRED, SINCE SNWLA 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. */
svsort_(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] = snrm2_(n, &vec[i * *nmvec], &c__1);
}
/* THIS PERFORMS THE ORTHONORMALIZATION. */
m = *n * *nblock;
sortqr_(nmvec, n, nperm, vec, &work[m]);
for (i = 0; i < *nperm; ++i, m += *nperm + 1) {
if (abs(work[m]) <= val[i + *nmval * 2] * .9f) {
*ierr = -1;
return;
}
}
/* THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN */
/* THE ARRAY WORK FOR LATER REFERENCE IN SNWLA. */
m = (*n << 1) * *nblock;
scopy_(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.f;
for (i = 0; i < 1000; ++i) {
eps *= .5f;
temp = eps + 1.f;
if (temp == 1.f) {
break;
}
}
/* ------------------------------------------------------------------ */
/* THIS SECTION CALLS SNWLA 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);
i13 = i12 + *maxj;
snwla_(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);
/* ------------------------------------------------------------------ */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?