📄 dnlaso.c
字号:
/* laso/dnlaso.f -- translated by f2c (version 20050501).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"
/* Table of constant values */
static integer c__0 = 0;
static integer c__1 = 1;
static doublereal c_b15 = 10.;
static doublereal c_b88 = 0.;
/* VERSION 2 DOES NOT USE EISPACK */
/* ------------------------------------------------------------------ */
/*< >*/
/* Subroutine */ int dnlaso_(
void (*op)(integer*,integer*,doublereal*,doublereal*),
void (*iovect)(integer*,integer*,doublereal*,integer*,integer*),
integer *n, integer *nval,
integer *nfig, integer *nperm, integer *nmval, doublereal *val,
integer *nmvec, doublereal *vec, integer *nblock, integer *maxop,
integer *maxj, doublereal *work, integer *ind, integer *ierr)
{
/* System generated locals */
integer vec_dim1, vec_offset, val_dim1, val_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer i__, m, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13,
nv;
doublereal eps;
integer nop;
doublereal temp, tarr[1];
extern doublereal dnrm2_(integer *, doublereal *, integer *);
integer nband;
doublereal delta;
extern /* Subroutine */ int dnwla_(
void (*op)(integer*,integer*,doublereal*,doublereal*),
void (*iovect)(integer*,integer*,doublereal*,integer*,integer*),
integer *, integer *,
integer *, integer *, integer *, doublereal *, integer *,
doublereal *, integer *, integer *, integer *, integer *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, integer *, logical *, logical *,
doublereal *, doublereal *, integer *);
logical small;
extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
doublereal *, integer *),
dnppla_(void (*op)(integer*,integer*,doublereal*,doublereal*),
void (*iovect)(integer*,integer*,doublereal*,integer*,integer*),
integer *, integer *
, integer *, integer *, doublereal *, integer *, doublereal *,
integer *, doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, logical *, logical *,
doublereal *);
logical raritz;
extern /* Subroutine */ int dortqr_(integer *, integer *, integer *,
doublereal *, doublereal *), dvsort_(integer *, doublereal *,
doublereal *, integer *, doublereal *, integer *, integer *,
doublereal *);
/*< >*/
/*< DOUBLE PRECISION VEC(NMVEC,1), VAL(NMVAL,4), WORK(1) >*/
/*< EXTERNAL OP, IOVECT >*/
/* 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 DOUBLE PRECISIONISTIC 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. */
/*< >*/
/*< LOGICAL RARITZ, SMALL >*/
/*< DOUBLE PRECISION DELTA, EPS, TEMP, DNRM2, DABS, TARR(1) >*/
/*< EXTERNAL DNPPLA, DNWLA, DORTQR, DCOPY, DNRM2, DVSORT >*/
/* 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 = IABS(NVAL) >*/
/* Parameter adjustments */
val_dim1 = *nmval;
val_offset = 1 + val_dim1;
val -= val_offset;
vec_dim1 = *nmvec;
vec_offset = 1 + vec_dim1;
vec -= vec_offset;
--work;
--ind;
/* Function Body */
nv = abs(*nval);
/*< IND(1) = 0 >*/
ind[1] = 0;
/*< IERR = 0 >*/
*ierr = 0;
/*< IF (N.LT.6*NBLOCK) IERR = 1 >*/
if (*n < *nblock * 6) {
*ierr = 1;
}
/*< IF (NFIG.LE.0) IERR = IERR + 2 >*/
if (*nfig <= 0) {
*ierr += 2;
}
/*< IF (NMVEC.LT.N) IERR = IERR + 4 >*/
if (*nmvec < *n) {
*ierr += 4;
}
/*< IF (NPERM.LT.0) IERR = IERR + 8 >*/
if (*nperm < 0) {
*ierr += 8;
}
/*< IF (MAXJ.LT.6*NBLOCK) IERR = IERR + 16 >*/
if (*maxj < *nblock * 6) {
*ierr += 16;
}
/*< IF (NV.LT.MAX0(1,NPERM)) IERR = IERR + 32 >*/
if (nv < max(1,*nperm)) {
*ierr += 32;
}
/*< IF (NV.GT.NMVAL) IERR = IERR + 64 >*/
if (nv > *nmval) {
*ierr += 64;
}
/*< IF (NV.GT.MAXOP) IERR = IERR + 128 >*/
if (nv > *maxop) {
*ierr += 128;
}
/*< IF (NV.GE.MAXJ/2) IERR = IERR + 256 >*/
if (nv >= *maxj / 2) {
*ierr += 256;
}
/*< IF (NBLOCK.LT.1) IERR = IERR + 512 >*/
if (*nblock < 1) {
*ierr += 512;
}
/*< IF (IERR.NE.0) RETURN >*/
if (*ierr != 0) {
return 0;
}
/*< SMALL = NVAL.LT.0 >*/
small = *nval < 0;
/* ------------------------------------------------------------------ */
/* THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS. */
/* IF A USER SUPPLIED VECTOR IS ZERO OR IF DSIGNIFICANT CANCELLATION */
/* OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO -1 */
/* AND DNLASO TERMINATES. */
/*< IF (NPERM.EQ.0) GO TO 110 >*/
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) GO TO 20 >*/
if (small) {
goto L20;
}
/*< DO 10 I=1,NPERM >*/
i__1 = *nperm;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< VAL(I,1) = -VAL(I,1) >*/
val[i__ + val_dim1] = -val[i__ + val_dim1];
/*< 10 CONTINUE >*/
/* L10: */
}
/* THIS SORTS THE USER SUPPLIED VALUES AND VECTORS. */
/*< 20 CALL DVSORT(NPERM, VAL, VAL(1,2), 0, TARR, NMVEC, N, VEC) >*/
L20:
dvsort_(nperm, &val[val_offset], &val[(val_dim1 << 1) + 1], &c__0, tarr,
nmvec, n, &vec[vec_offset]);
/* THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON. */
/* IT ALSO INSURES THAT THE RESIDUAL NORMS ARE POSITIVE. */
/*< DO 60 I=1,NPERM >*/
i__1 = *nperm;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< VAL(I,2) = DABS(VAL(I,2)) >*/
val[i__ + (val_dim1 << 1)] = (d__1 = val[i__ + (val_dim1 << 1)], abs(
d__1));
/*< VAL(I,3) = DNRM2(N,VEC(1,I),1) >*/
val[i__ + val_dim1 * 3] = dnrm2_(n, &vec[i__ * vec_dim1 + 1], &c__1);
/*< 60 CONTINUE >*/
/* L60: */
}
/* THIS PERFORMS THE ORTHONORMALIZATION. */
/*< M = N*NBLOCK + 1 >*/
m = *n * *nblock + 1;
/*< CALL DORTQR(NMVEC, N, NPERM, VEC, WORK(M)) >*/
dortqr_(nmvec, n, nperm, &vec[vec_offset], &work[m]);
/*< M = N*NBLOCK - NPERM >*/
m = *n * *nblock - *nperm;
/*< DO 70 I = 1, NPERM >*/
i__1 = *nperm;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< M = M + NPERM + 1 >*/
m = m + *nperm + 1;
/*< IF(DABS(WORK(M)) .GT. 0.9*VAL(I,3)) GO TO 70 >*/
if ((d__1 = work[m], abs(d__1)) > val[i__ + val_dim1 * 3] * (float).9)
{
goto L70;
}
/*< IERR = -1 >*/
*ierr = -1;
/*< RETURN >*/
return 0;
/*< 70 CONTINUE >*/
L70:
;
}
/* THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN */
/* THE ARRAY WORK FOR LATER REFERENCE IN DNWLA. */
/*< M = 2*N*NBLOCK + 1 >*/
m = (*n << 1) * *nblock + 1;
/*< CALL DCOPY(NPERM, VAL(1,2), 1, WORK(M), 1) >*/
dcopy_(nperm, &val[(val_dim1 << 1) + 1], &c__1, &work[m], &c__1);
/* THIS SETS EPS TO AN APPROXIMATION OF THE RELATIVE MACHINE */
/* PRECISION */
/* ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT */
/* ***IN A PRODUCTION CODE */
/*< 110 EPS = 1.0D0 >*/
L110:
eps = 1.;
/*< DO 120 I = 1,1000 >*/
for (i__ = 1; i__ <= 1000; ++i__) {
/*< EPS = 0.5D0*EPS >*/
eps *= .5;
/*< TEMP = 1.0D0 + EPS >*/
temp = eps + 1.;
/*< IF(TEMP.EQ.1.0D0) GO TO 130 >*/
if (temp == 1.) {
goto L130;
}
/*< 120 CONTINUE >*/
/* L120: */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -