1C VERSION 2 DOES NOT USE EISPACK 2C 3C ------------------------------------------------------------------ 4C 5 SUBROUTINE DNLASO(OP, IOVECT, N, NVAL, NFIG, NPERM, 6 * NMVAL, VAL, NMVEC, VEC, NBLOCK, MAXOP, MAXJ, WORK, 7 * IND, IERR) 8C 9 INTEGER N, NVAL, NFIG, NPERM, NMVAL, NMVEC, NBLOCK, 10 * MAXOP, MAXJ, IND(1), IERR 11 DOUBLE PRECISION VEC(NMVEC,1), VAL(NMVAL,4), WORK(1) 12 EXTERNAL OP, IOVECT 13C 14C AUTHOR/IMPLEMENTER D.S.SCOTT-B.N.PARLETT/D.S.SCOTT 15C 16C COMPUTER SCIENCES DEPARTMENT 17C UNIVERSITY OF TEXAS AT AUSTIN 18C AUSTIN, TX 78712 19C 20C VERSION 2 ORIGINATED APRIL 1982 21C 22C CURRENT VERSION JUNE 1983 23C 24C DNLASO FINDS A FEW EIGENVALUES AND EIGENVECTORS AT EITHER END OF 25C THE SPECTRUM OF A LARGE SPARSE SYMMETRIC MATRIX. THE SUBROUTINE 26C DNLASO IS PRIMARILY A DRIVER FOR SUBROUTINE DNWLA WHICH IMPLEMENTS 27C THE LANCZOS ALGORITHM WITH SELECTIVE ORTHOGONALIZATION AND 28C SUBROUTINE DNPPLA WHICH POST PROCESSES THE OUTPUT OF DNWLA. 29C HOWEVER DNLASO DOES CHECK FOR INCONSISTENCIES IN THE CALLING 30C PARAMETERS AND DOES PREPROCESS ANY USER SUPPLIED EIGENPAIRS. 31C DNLASO ALWAYS LOOKS FOR THE SMALLEST (LEFTMOST) EIGENVALUES. IF 32C THE LARGEST EIGENVALUES ARE DESIRED DNLASO IMPLICITLY USES THE 33C NEGATIVE OF THE MATRIX. 34C 35C 36C ON INPUT 37C 38C 39C OP A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE 40C OP(N,M,P,Q). P AND Q ARE N X M MATRICES AND Q IS 41C RETURNED AS THE MATRIX TIMES P. 42C 43C IOVECT A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE 44C IOVECT(N,M,Q,J,K). Q IS AN N X M MATRIX. IF K = 0 45C THE COLUMNS OF Q ARE STORED AS THE (J-M+1)TH THROUGH 46C THE JTH LANCZOS VECTORS. IF K = 1 THEN Q IS RETURNED 47C AS THE (J-M+1)TH THROUGH THE JTH LANCZOS VECTORS. SEE 48C DOCUMENTATION FOR FURTHER DETAILS AND EXAMPLES. 49C 50C N THE ORDER OF THE MATRIX. 51C 52C NVAL NVAL SPECIFIES THE EIGENVALUES TO BE FOUND. 53C DABS(NVAL) IS THE NUMBER OF EIGENVALUES DESIRED. 54C IF NVAL < 0 THE ALGEBRAICALLY SMALLEST (LEFTMOST) 55C EIGENVALUES ARE FOUND. IF NVAL > 0 THE ALGEBRAICALLY 56C LARGEST (RIGHTMOST) EIGENVALUES ARE FOUND. NVAL MUST NOT 57C BE ZERO. DABS(NVAL) MUST BE LESS THAN MAXJ/2. 58C 59C NFIG THE NUMBER OF DECIMAL DIGITS OF ACCURACY DESIRED IN THE 60C EIGENVALUES. NFIG MUST BE GREATER THAN OR EQUAL TO 1. 61C 62C NPERM AN INTEGER VARIABLE WHICH SPECIFIES THE NUMBER OF USER 63C SUPPLIED EIGENPAIRS. IN MOST CASES NPERM WILL BE ZERO. SEE 64C DOCUMENTAION FOR FURTHER DETAILS OF USING NPERM GREATER 65C THAN ZERO. NPERM MUST NOT BE LESS THAN ZERO. 66C 67C NMVAL THE ROW DIMENSION OF THE ARRAY VAL. NMVAL MUST BE GREATER 68C THAN OR EQUAL TO DABS(NVAL). 69C 70C VAL A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW 71C DIMENSION NMVAL AND COLUMN DIMENSION AT LEAST 4. IF NPERM 72C IS GREATER THAN ZERO THEN CERTAIN INFORMATION MUST BE STORED 73C IN VAL. SEE DOCUMENTATION FOR DETAILS. 74C 75C NMVEC THE ROW DIMENSION OF THE ARRAY VEC. NMVEC MUST BE GREATER 76C THAN OR EQUAL TO N. 77C 78C VEC A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW 79C DIMENSION NMVEC AND COLUMN DIMENSION AT LEAST DABS(NVAL). IF 80C NPERM > 0 THEN THE FIRST NPERM COLUMNS OF VEC MUST 81C CONTAIN THE USER SUPPLIED EIGENVECTORS. 82C 83C NBLOCK THE BLOCK SIZE. SEE DOCUMENTATION FOR CHOOSING 84C AN APPROPRIATE VALUE FOR NBLOCK. NBLOCK MUST BE GREATER 85C THAN ZERO AND LESS THAN MAXJ/6. 86C 87C MAXOP AN UPPER BOUND ON THE NUMBER OF CALLS TO THE SUBROUTINE 88C OP. DNLASO TERMINATES WHEN MAXOP IS EXCEEDED. SEE 89C DOCUMENTATION FOR GUIDELINES IN CHOOSING A VALUE FOR MAXOP. 90C 91C MAXJ AN INDICATION OF THE AVAILABLE STORAGE (SEE WORK AND 92C DOCUMENTATION ON IOVECT). FOR THE FASTEST CONVERGENCE MAXJ 93C SHOULD BE AS LARGE AS POSSIBLE, ALTHOUGH IT IS USELESS TO HAVE 94C MAXJ LARGER THAN MAXOP*NBLOCK. 95C 96C WORK A DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST AS 97C LARGE AS 98C 99C 2*N*NBLOCK + MAXJ*(NBLOCK+NV+2) + 2*NBLOCK*NBLOCK + 3*NV 100C 101C + THE MAXIMUM OF 102C N*NBLOCK 103C AND 104C MAXJ*(2*NBLOCK+3) + 2*NV + 6 + (2*NBLOCK+2)*(NBLOCK+1) 105C 106C WHERE NV = DABS(NVAL) 107C 108C THE FIRST N*NBLOCK ELEMENTS OF WORK MUST CONTAIN THE DESIRED 109C STARTING VECTORS. SEE DOCUMENTATION FOR GUIDELINES IN 110C CHOOSING STARTING VECTORS. 111C 112C IND AN INTEGER ARRAY OF DIMENSION AT LEAST DABS(NVAL). 113C 114C IERR AN INTEGER VARIABLE. 115C 116C 117C ON OUTPUT 118C 119C 120C NPERM THE NUMBER OF EIGENPAIRS NOW KNOWN. 121C 122C VEC THE FIRST NPERM COLUMNS OF VEC CONTAIN THE EIGENVECTORS. 123C 124C VAL THE FIRST COLUMN OF VAL CONTAINS THE CORRESPONDING 125C EIGENVALUES. THE SECOND COLUMN CONTAINS THE RESIDUAL NORMS OF 126C THE EIGENPAIRS WHICH ARE BOUNDS ON THE ACCURACY OF THE EIGEN- 127C VALUES. THE THIRD COLUMN CONTAINS MORE DOUBLE PRECISIONISTIC ESTIMATES 128C OF THE ACCURACY OF THE EIGENVALUES. THE FOURTH COLUMN CONTAINS 129C ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS. SEE 130C DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES. 131C 132C WORK IF WORK IS TERMINATED BEFORE COMPLETION (IERR = -2) 133C THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS 134C FOR RESTARTING THE ALGORITHM AND DNLASO CAN BE IMMEDIATELY 135C RECALLED TO CONTINUE WORKING ON THE PROBLEM. 136C 137C IND IND(1) CONTAINS THE ACTUAL NUMBER OF CALLS TO OP. ON SOME 138C OCCASIONS THE NUMBER OF CALLS TO OP MAY BE SLIGHTLY LARGER 139C THAN MAXOP. 140C 141C IERR AN ERROR COMPLETION CODE. THE NORMAL COMPLETION CODE IS 142C ZERO. SEE THE DOCUMENTATION FOR INTERPRETATIONS OF NON-ZERO 143C COMPLETION CODES. 144C 145C 146C INTERNAL VARIABLES. 147C 148C 149 INTEGER I, I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11, 150 * I12, I13, M, NBAND, NOP, NV, IABS, MAX0 151 LOGICAL RARITZ, SMALL 152 DOUBLE PRECISION DELTA, EPS, TEMP, DNRM2, DABS, TARR(1) 153 EXTERNAL DNPPLA, DNWLA, DORTQR, DCOPY, DNRM2, DVSORT 154C 155C NOP RETURNED FROM DNWLA AS THE NUMBER OF CALLS TO THE 156C SUBROUTINE OP. 157C 158C NV SET EQUAL TO DABS(NVAL), THE NUMBER OF EIGENVALUES DESIRED, 159C AND PASSED TO DNWLA. 160C 161C SMALL SET TO .TRUE. IF THE SMALLEST EIGENVALUES ARE DESIRED. 162C 163C RARITZ RETURNED FROM DNWLA AND PASSED TO DNPPLA. RARITZ IS .TRUE. 164C IF A FINAL RAYLEIGH-RITZ PROCEDURE IS NEEDED. 165C 166C DELTA RETURNED FROM DNWLA AS THE EIGENVALUE OF THE MATRIX 167C WHICH IS CLOSEST TO THE DESIRED EIGENVALUES. 168C 169C DNPPLA A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED 170C BY DNWLA. 171C 172C DNWLA A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM 173C WITH SELECTIVE ORTHOGONALIZATION. 174C 175C DMVPC A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND 176C ORTHOGONALITY COEFFICIENT OF GIVEN RITZ VECTORS. 177C 178C DORTQR A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS 179C USING HOUSEHOLDER REFLECTIONS. 180C 181C DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP A SUBSET OF THE BASIC LINEAR 182C ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION. 183C 184C DLARAN A SUBROUTINE TO GENERATE RANDOM VECTORS 185C 186C DLAEIG, DLAGER, DLABCM, DLABFC SUBROUTINES FOR BAND EIGENVALUE 187C CALCULATIONS. 188C 189C ------------------------------------------------------------------ 190C 191C THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS. 192C 193 NV = IABS(NVAL) 194 IND(1) = 0 195 IERR = 0 196 IF (N.LT.6*NBLOCK) IERR = 1 197 IF (NFIG.LE.0) IERR = IERR + 2 198 IF (NMVEC.LT.N) IERR = IERR + 4 199 IF (NPERM.LT.0) IERR = IERR + 8 200 IF (MAXJ.LT.6*NBLOCK) IERR = IERR + 16 201 IF (NV.LT.MAX0(1,NPERM)) IERR = IERR + 32 202 IF (NV.GT.NMVAL) IERR = IERR + 64 203 IF (NV.GT.MAXOP) IERR = IERR + 128 204 IF (NV.GE.MAXJ/2) IERR = IERR + 256 205 IF (NBLOCK.LT.1) IERR = IERR + 512 206 IF (IERR.NE.0) RETURN 207C 208 SMALL = NVAL.LT.0 209C 210C ------------------------------------------------------------------ 211C 212C THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS. 213C IF A USER SUPPLIED VECTOR IS ZERO OR IF DSIGNIFICANT CANCELLATION 214C OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO -1 215C AND DNLASO TERMINATES. 216C 217 IF (NPERM.EQ.0) GO TO 110 218C 219C THIS NEGATES THE USER SUPPLIED EIGENVALUES WHEN THE LARGEST 220C EIGENVALUES ARE DESIRED, SINCE DNWLA WILL IMPLICITLY USE THE 221C NEGATIVE OF THE MATRIX. 222C 223 IF (SMALL) GO TO 20 224 DO 10 I=1,NPERM 225 VAL(I,1) = -VAL(I,1) 226 10 CONTINUE 227C 228C THIS SORTS THE USER SUPPLIED VALUES AND VECTORS. 229C 230 20 CALL DVSORT(NPERM, VAL, VAL(1,2), 0, TARR, NMVEC, N, VEC) 231C 232C THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON. 233C IT ALSO INSURES THAT THE RESIDUAL NORMS ARE POSITIVE. 234C 235 DO 60 I=1,NPERM 236 VAL(I,2) = DABS(VAL(I,2)) 237 VAL(I,3) = DNRM2(N,VEC(1,I),1) 238 60 CONTINUE 239C 240C THIS PERFORMS THE ORTHONORMALIZATION. 241C 242 M = N*NBLOCK + 1 243 CALL DORTQR(NMVEC, N, NPERM, VEC, WORK(M)) 244 M = N*NBLOCK - NPERM 245 DO 70 I = 1, NPERM 246 M = M + NPERM + 1 247 IF(DABS(WORK(M)) .GT. 0.9*VAL(I,3)) GO TO 70 248 IERR = -1 249 RETURN 250C 251 70 CONTINUE 252C 253C THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN 254C THE ARRAY WORK FOR LATER REFERENCE IN DNWLA. 255C 256 M = 2*N*NBLOCK + 1 257 CALL DCOPY(NPERM, VAL(1,2), 1, WORK(M), 1) 258C 259C THIS SETS EPS TO AN APPROXIMATION OF THE RELATIVE MACHINE 260C PRECISION 261C 262C ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT 263C ***IN A PRODUCTION CODE 264C 265 110 EPS = 1.0D0 266 DO 120 I = 1,1000 267 EPS = 0.5D0*EPS 268 TEMP = 1.0D0 + EPS 269 IF(TEMP.EQ.1.0D0) GO TO 130 270 120 CONTINUE 271C 272C ------------------------------------------------------------------ 273C 274C THIS SECTION CALLS DNWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM 275C WITH SELECTIVE ORTHOGONALIZATION. 276C 277 130 NBAND = NBLOCK + 1 278 I1 = 1 + N*NBLOCK 279 I2 = I1 + N*NBLOCK 280 I3 = I2 + NV 281 I4 = I3 + NV 282 I5 = I4 + NV 283 I6 = I5 + MAXJ*NBAND 284 I7 = I6 + NBLOCK*NBLOCK 285 I8 = I7 + NBLOCK*NBLOCK 286 I9 = I8 + MAXJ*(NV+1) 287 I10 = I9 + NBLOCK 288 I11 = I10 + 2*NV + 6 289 I12 = I11 + MAXJ*(2*NBLOCK+1) 290 I13 = I12 + MAXJ 291 CALL DNWLA(OP, IOVECT, N, NBAND, NV, NFIG, NPERM, VAL, NMVEC, 292 * VEC, NBLOCK, MAXOP, MAXJ, NOP, WORK(1), WORK(I1), 293 * WORK(I2), WORK(I3), WORK(I4), WORK(I5), WORK(I6), 294 * WORK(I7), WORK(I8), WORK(I9), WORK(I10), WORK(I11), 295 * WORK(I12), WORK(I13), IND, SMALL, RARITZ, DELTA, EPS, IERR) 296C 297C ------------------------------------------------------------------ 298C 299C THIS SECTION CALLS DNPPLA (THE POST PROCESSOR). 300C 301 IF (NPERM.EQ.0) GO TO 140 302 I1 = N*NBLOCK + 1 303 I2 = I1 + NPERM*NPERM 304 I3 = I2 + NPERM*NPERM 305 I4 = I3 + MAX0(N*NBLOCK,2*NPERM*NPERM) 306 I5 = I4 + N*NBLOCK 307 I6 = I5 + 2*NPERM + 4 308 CALL DNPPLA(OP, IOVECT, N, NPERM, NOP, NMVAL, VAL, NMVEC, 309 * VEC, NBLOCK, WORK(I1), WORK(I2), WORK(I3), WORK(I4), 310 * WORK(I5), WORK(I6), DELTA, SMALL, RARITZ, EPS) 311C 312 140 IND(1) = NOP 313 RETURN 314 END 315C 316C ------------------------------------------------------------------ 317C 318 SUBROUTINE DNWLA(OP, IOVECT, N, NBAND, NVAL, NFIG, NPERM, VAL, 319 * NMVEC, VEC, NBLOCK, MAXOP, MAXJ, NOP, P1, P0, 320 * RES, TAU, OTAU, T, ALP, BET, S, P2, BOUND, ATEMP, VTEMP, D, 321 * IND, SMALL, RARITZ, DELTA, EPS, IERR) 322C 323 INTEGER N, NBAND, NVAL, NFIG, NPERM, NMVEC, NBLOCK, MAXOP, MAXJ, 324 * NOP, IND(1), IERR 325 LOGICAL RARITZ, SMALL 326 DOUBLE PRECISION VAL(1), VEC(NMVEC,1), P0(N,1), P1(N,1), 327 * P2(N,1), RES(1), TAU(1), OTAU(1), T(NBAND,1), 328 * ALP(NBLOCK,1), BET(NBLOCK,1), BOUND(1), ATEMP(1), 329 * VTEMP(1), D(1), S(MAXJ,1), DELTA, EPS 330 EXTERNAL OP, IOVECT 331C 332C DNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE 333C ORTHOGONALIZATION. 334C 335C NBAND NBLOCK + 1 THE BAND WIDTH OF T. 336C 337C NVAL THE NUMBER OF DESIRED EIGENVALUES. 338C 339C NPERM THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS 340C INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE 341C ALGORITHM IS ITERATED). PERMANENT VECTORS ARE ORTHOGONAL 342C TO THE CURRENT KRYLOV SUBSPACE. 343C 344C NOP THE NUMBER OF CALLS TO OP. 345C 346C P0, P1, AND P2 THE CURRENT BLOCKS OF LANCZOS VECTORS. 347C 348C RES THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS. 349C 350C TAU AND OTAU USED TO MONITOR THE NEED FOR ORTHOGONALIZATION. 351C 352C T THE BAND MATRIX. 353C 354C ALP THE CURRENT DIAGONAL BLOCK. 355C 356C BET THE CURRENT OFF DIAGONAL BLOCK. 357C 358C BOUND, ATEMP, VTEMP, D TEMPORARY STORAGE USED BY THE BAND 359C EIGENVALUE SOLVER DLAEIG. 360C 361C S EIGENVECTORS OF T. 362C 363C SMALL .TRUE. IF THE SMALL EIGENVALUES ARE DESIRED. 364C 365C RARITZ RETURNED AS .TRUE. IF A FINAL RAYLEIGH-RITZ PROCEDURE 366C IS TO BE DONE. 367C 368C DELTA RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE 369C OF THE MATRIX. USED IN ESTIMATING THE ACCURACY OF THE 370C COMPUTED EIGENVALUES. 371C 372C 373C INTERNAL VARIABLES USED 374C 375 INTEGER I, I1, II, J, K, L, M, NG, NGOOD, 376 * NLEFT, NSTART, NTHETA, NUMBER, MIN0, MTEMP 377 LOGICAL ENOUGH, TEST 378 DOUBLE PRECISION ALPMAX, ALPMIN, ANORM, BETMAX, BETMIN, 379 * EPSRT, PNORM, RNORM, TEMP, 380 * TMAX, TMIN, TOLA, TOLG, UTOL, DABS, 381 * DMAX1, DMIN1, DSQRT, DDOT, DNRM2, TARR(1), DZERO(1) 382 EXTERNAL DMVPC, DORTQR, DAXPY, DCOPY, DDOT, 383 * DNRM2, DSCAL, DLAEIG, DLAGER, DLARAN, DVSORT 384C 385C J THE CURRENT DIMENSION OF T. (THE DIMENSION OF THE CURRENT 386C KRYLOV SUBSPACE. 387C 388C NGOOD THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS 389C LIE IN THE CURRENT KRYLOV SUBSPACE). 390C 391C NLEFT THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED, 392C I.E. NLEFT = NVAL - NPERM. 393C 394C NUMBER = NPERM + NGOOD. 395C 396C ANORM AN ESTIMATE OF THE NORM OF THE MATRIX. 397C 398C EPS THE RELATIVE MACHINE PRECISION. 399C 400C UTOL THE USER TOLERANCE. 401C 402C TARR AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN CALLS TO 403C DLAEIG 404C 405C DZERO AN ARRAY OF LENGTH ONE CONTAINING DZERO, USED TO INSURE TYPE 406C CONSISTENCY IN CALLS TO DCOPY 407C 408 DZERO(1) = 0.0D0 409 RNORM = 0.0D0 410 IF (NPERM.NE.0) RNORM = DMAX1(-VAL(1),VAL(NPERM)) 411 PNORM = RNORM 412 DELTA = 10.D30 413 EPSRT = DSQRT(EPS) 414 NLEFT = NVAL - NPERM 415 NOP = 0 416 NUMBER = NPERM 417 RARITZ = .FALSE. 418 UTOL = DMAX1(DBLE(FLOAT(N))*EPS,10.0D0**DBLE((-FLOAT(NFIG)))) 419 J = MAXJ 420C 421C ------------------------------------------------------------------ 422C 423C ANY ITERATION OF THE ALGORITHM BEGINS HERE. 424C 425 30 DO 50 I=1,NBLOCK 426 TEMP = DNRM2(N,P1(1,I),1) 427 IF (TEMP.EQ.0D0) CALL DLARAN(N, P1(1,I)) 428 50 CONTINUE 429 IF (NPERM.EQ.0) GO TO 70 430 DO 60 I=1,NPERM 431 TAU(I) = 1.0D0 432 OTAU(I) = 0.0D0 433 60 CONTINUE 434 70 CALL DCOPY(N*NBLOCK, DZERO, 0, P0, 1) 435 CALL DCOPY(NBLOCK*NBLOCK, DZERO, 0, BET, 1) 436 CALL DCOPY(J*NBAND, DZERO, 0, T, 1) 437 MTEMP = NVAL + 1 438 DO 75 I = 1, MTEMP 439 CALL DCOPY(J, DZERO, 0, S(1,I), 1) 440 75 CONTINUE 441 NGOOD = 0 442 TMIN = 1.0D30 443 TMAX = -1.0D30 444 TEST = .TRUE. 445 ENOUGH = .FALSE. 446 BETMAX = 0.0D0 447 J = 0 448C 449C ------------------------------------------------------------------ 450C 451C THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP. 452C 453 80 J = J + NBLOCK 454C 455C THIS IS THE SELECTIVE ORTHOGONALIZATION. 456C 457 IF (NUMBER.EQ.0) GO TO 110 458 DO 100 I=1,NUMBER 459 IF (TAU(I).LT.EPSRT) GO TO 100 460 TEST = .TRUE. 461 TAU(I) = 0.0D0 462 IF (OTAU(I).NE.0.0D0) OTAU(I) = 1.0D0 463 DO 90 K=1,NBLOCK 464 TEMP = -DDOT(N,VEC(1,I),1,P1(1,K),1) 465 CALL DAXPY(N, TEMP, VEC(1,I), 1, P1(1,K), 1) 466C 467C THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A 468C NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR. THE ALGORITHM IS 469C TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST. 470C 471 IF (DABS(TEMP*BET(K,K)).GT.DBLE(FLOAT(N))*EPSRT* 472 * ANORM .AND. I.GT.NPERM) GO TO 380 473 90 CONTINUE 474 100 CONTINUE 475C 476C IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET. 477C 478 110 IF(.NOT. TEST) GO TO 160 479 CALL DORTQR(N, N, NBLOCK, P1, ALP) 480 TEST = .FALSE. 481 IF(J .EQ. NBLOCK) GO TO 160 482 DO 130 I = 1,NBLOCK 483 IF(ALP(I,I) .GT. 0.0D0) GO TO 130 484 M = J - 2*NBLOCK + I 485 L = NBLOCK + 1 486 DO 120 K = I,NBLOCK 487 BET(I,K) = -BET(I,K) 488 T(L,M) = -T(L,M) 489 L = L - 1 490 M = M + 1 491 120 CONTINUE 492 130 CONTINUE 493C 494C THIS IS THE LANCZOS STEP. 495C 496 160 CALL OP(N, NBLOCK, P1, P2) 497 NOP = NOP + 1 498 CALL IOVECT(N, NBLOCK, P1, J, 0) 499C 500C THIS COMPUTES P2=P2-P0*BET(TRANSPOSE) 501C 502 DO 180 I=1,NBLOCK 503 DO 170 K=I,NBLOCK 504 CALL DAXPY(N, -BET(I,K), P0(1,K), 1, P2(1,I), 1) 505 170 CONTINUE 506 180 CONTINUE 507C 508C THIS COMPUTES ALP AND P2=P2-P1*ALP. 509C 510 DO 200 I=1,NBLOCK 511 DO 190 K=1,I 512 II = I - K + 1 513 ALP(II,K) = DDOT(N,P1(1,I),1,P2(1,K),1) 514 CALL DAXPY(N, -ALP(II,K), P1(1,I), 1, P2(1,K), 1) 515 IF (K.NE.I) CALL DAXPY(N, -ALP(II,K), P1(1,K), 516 * 1, P2(1,I), 1) 517 190 CONTINUE 518 200 CONTINUE 519C 520C REORTHOGONALIZATION OF THE SECOND BLOCK 521C 522 IF(J .NE. NBLOCK) GO TO 220 523 DO 215 I=1,NBLOCK 524 DO 210 K=1,I 525 TEMP = DDOT(N,P1(1,I),1,P2(1,K),1) 526 CALL DAXPY(N, -TEMP, P1(1,I), 1, P2(1,K), 1) 527 IF (K.NE.I) CALL DAXPY(N, -TEMP, P1(1,K), 528 * 1, P2(1,I), 1) 529 II = I - K + 1 530 ALP(II,K) = ALP(II,K) + TEMP 531 210 CONTINUE 532 215 CONTINUE 533C 534C THIS ORTHONORMALIZES THE NEXT BLOCK 535C 536 220 CALL DORTQR(N, N, NBLOCK, P2, BET) 537C 538C THIS STORES ALP AND BET IN T. 539C 540 DO 250 I=1,NBLOCK 541 M = J - NBLOCK + I 542 DO 230 K=I,NBLOCK 543 L = K - I + 1 544 T(L,M) = ALP(L,I) 545 230 CONTINUE 546 DO 240 K=1,I 547 L = NBLOCK - I + K + 1 548 T(L,M) = BET(K,I) 549 240 CONTINUE 550 250 CONTINUE 551C 552C THIS NEGATES T IF SMALL IS FALSE. 553C 554 IF (SMALL) GO TO 280 555 M = J - NBLOCK + 1 556 DO 270 I=M,J 557 DO 260 K=1,L 558 T(K,I) = -T(K,I) 559 260 CONTINUE 560 270 CONTINUE 561C 562C THIS SHIFTS THE LANCZOS VECTORS 563C 564 280 CALL DCOPY(NBLOCK*N, P1, 1, P0, 1) 565 CALL DCOPY(NBLOCK*N, P2, 1, P1, 1) 566 CALL DLAGER(J, NBAND, J-NBLOCK+1, T, TMIN, TMAX) 567 ANORM = DMAX1(RNORM, TMAX, -TMIN) 568 IF (NUMBER.EQ.0) GO TO 305 569C 570C THIS COMPUTES THE EXTREME EIGENVALUES OF ALP. 571C 572 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) 573 CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK, 574 1 P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX) 575 ALPMIN = TARR(1) 576 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) 577 CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR, 578 1 NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX) 579 ALPMAX = TARR(1) 580C 581C THIS COMPUTES ALP = BET(TRANSPOSE)*BET. 582C 583 305 DO 310 I = 1, NBLOCK 584 DO 300 K = 1, I 585 L = I - K + 1 586 ALP(L,K) = DDOT(NBLOCK-I+1, BET(I,I), NBLOCK, BET(K,I), 587 1 NBLOCK) 588 300 CONTINUE 589 310 CONTINUE 590 IF(NUMBER .EQ. 0) GO TO 330 591C 592C THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET. 593C 594 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) 595 CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK, 596 1 P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0, ANORM*ANORM) 597 BETMIN = DSQRT(TARR(1)) 598C 599C THIS UPDATES TAU AND OTAU. 600C 601 DO 320 I=1,NUMBER 602 TEMP = (TAU(I)*DMAX1(ALPMAX-VAL(I),VAL(I)-ALPMIN) 603 * +OTAU(I)*BETMAX+EPS*ANORM)/BETMIN 604 IF (I.LE.NPERM) TEMP = TEMP + RES(I)/BETMIN 605 OTAU(I) = TAU(I) 606 TAU(I) = TEMP 607 320 CONTINUE 608C 609C THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET. 610C 611 330 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) 612 CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP, TARR, 613 1 NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0, 614 2 ANORM*ANORM) 615 BETMAX = DSQRT(TARR(1)) 616 IF (J.LE.2*NBLOCK) GO TO 80 617C 618C ------------------------------------------------------------------ 619C 620C THIS SECTION COMPUTES AND EXAMINES THE SMALLEST NONGOOD AND 621C LARGEST DESIRED EIGENVALUES OF T TO SEE IF A CLOSER LOOK 622C IS JUSTIFIED. 623C 624 TOLG = EPSRT*ANORM 625 TOLA = UTOL*RNORM 626 IF(MAXJ-J .LT. NBLOCK .OR. (NOP .GE. MAXOP .AND. 627 1 NLEFT .NE. 0)) GO TO 390 628 GO TO 400 629 630C 631C ------------------------------------------------------------------ 632C 633C THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO 634C SEE IF FURTHER ACTION IS INDICATED, ENTRY IS AT 380 OR 390 IF AN 635C ITERATION (OR TERMINATION) IS KNOWN TO BE NEEDED, OTHERWISE ENTRY 636C IS AT 400. 637C 638 380 J = J - NBLOCK 639 IERR = -8 640 390 IF (NLEFT.EQ.0) RETURN 641 TEST = .TRUE. 642 400 NTHETA = MIN0(J/2,NLEFT+1) 643 CALL DLAEIG(J, NBAND, 1, NTHETA, T, VAL(NUMBER+1), 644 1 MAXJ, S, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX) 645 CALL DMVPC(NBLOCK, BET, MAXJ, J, S, NTHETA, ATEMP, VTEMP, D) 646C 647C THIS CHECKS FOR TERMINATION OF A CHECK RUN 648C 649 IF(NLEFT .NE. 0 .OR. J .LT. 6*NBLOCK) GO TO 410 650 IF(VAL(NUMBER+1)-ATEMP(1) .GT. VAL(NPERM) - TOLA) GO TO 790 651C 652C THIS UPDATES NLEFT BY EXAMINING THE COMPUTED EIGENVALUES OF T 653C TO DETERMINE IF SOME PERMANENT VALUES ARE NO LONGER DESIRED. 654C 655 410 IF (NTHETA.LE.NLEFT) GO TO 470 656 IF (NPERM.EQ.0) GO TO 430 657 M = NUMBER + NLEFT + 1 658 IF (VAL(M).GE.VAL(NPERM)) GO TO 430 659 NPERM = NPERM - 1 660 NGOOD = 0 661 NUMBER = NPERM 662 NLEFT = NLEFT + 1 663 GO TO 400 664C 665C THIS UPDATES DELTA. 666C 667 430 M = NUMBER + NLEFT + 1 668 DELTA = DMIN1(DELTA,VAL(M)) 669 ENOUGH = .TRUE. 670 IF(NLEFT .EQ. 0) GO TO 80 671 NTHETA = NLEFT 672 VTEMP(NTHETA+1) = 1 673C 674C ------------------------------------------------------------------ 675C 676C THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL. 677C 678C THIS CHECKS FOR ENOUGH ACCEPTABLE VALUES. 679C 680 IF (.NOT.(TEST .OR. ENOUGH)) GO TO 470 681 DELTA = DMIN1(DELTA,ANORM) 682 PNORM = DMAX1(RNORM,DMAX1(-VAL(NUMBER+1),DELTA)) 683 TOLA = UTOL*PNORM 684 NSTART = 0 685 DO 460 I=1,NTHETA 686 M = NUMBER + I 687 IF (DMIN1(ATEMP(I)*ATEMP(I)/(DELTA-VAL(M)),ATEMP(I)) 688 * .GT.TOLA) GO TO 450 689 IND(I) = -1 690 GO TO 460 691C 692 450 ENOUGH = .FALSE. 693 IF (.NOT.TEST) GO TO 470 694 IND(I) = 1 695 NSTART = NSTART + 1 696 460 CONTINUE 697C 698C COPY VALUES OF IND INTO VTEMP 699C 700 DO 465 I = 1,NTHETA 701 VTEMP(I) = DBLE(FLOAT(IND(I))) 702 465 CONTINUE 703 GO TO 500 704C 705C THIS CHECKS FOR NEW GOOD VECTORS. 706C 707 470 NG = 0 708 DO 490 I=1,NTHETA 709 IF (VTEMP(I).GT.TOLG) GO TO 480 710 NG = NG + 1 711 VTEMP(I) = -1 712 GO TO 490 713C 714 480 VTEMP(I) = 1 715 490 CONTINUE 716C 717 IF (NG.LE.NGOOD) GO TO 80 718 NSTART = NTHETA - NG 719C 720C ------------------------------------------------------------------ 721C 722C THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS. 723C IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED. 724C 725 500 TEST = TEST .AND. .NOT.ENOUGH 726 NGOOD = NTHETA - NSTART 727 NSTART = NSTART + 1 728 NTHETA = NTHETA + 1 729C 730C THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND 731C EIGENVECTORS OF T. THE OTHER EIGENVECTORS ARE SAVED FOR 732C FORMING STARTING VECTORS, IF NECESSARY. IT ALSO SHIFTS THE 733C EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS 734C PAUSE. 735C 736 CALL DCOPY(NTHETA, VAL(NUMBER+1), 1, VAL(NPERM+1), 1) 737 IF (NSTART.EQ.0) GO TO 580 738 IF (NSTART.EQ.NTHETA) GO TO 530 739 CALL DVSORT(NTHETA, VTEMP, ATEMP, 1, VAL(NPERM+1), MAXJ, 740 * J, S) 741C 742C THES ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING 743C VECTORS. 744C 745 530 IF (.NOT.TEST) NSTART = 0 746 IF (.NOT.TEST) GO TO 580 747C 748C FIND MINIMUM ATEMP VALUE TO AVOID POSSIBLE OVERFLOW 749C 750 TEMP = ATEMP(1) 751 DO 535 I = 1, NSTART 752 TEMP = DMIN1(TEMP, ATEMP(I)) 753 535 CONTINUE 754 M = NGOOD + 1 755 L = NGOOD + MIN0(NSTART,NBLOCK) 756 DO 540 I=M,L 757 CALL DSCAL(J, TEMP/ATEMP(I), S(1,I), 1) 758 540 CONTINUE 759 M = (NSTART-1)/NBLOCK 760 IF (M.EQ.0) GO TO 570 761 L = NGOOD + NBLOCK 762 DO 560 I=1,M 763 DO 550 K=1,NBLOCK 764 L = L + 1 765 IF (L.GT.NTHETA) GO TO 570 766 I1 = NGOOD + K 767 CALL DAXPY(J, TEMP/ATEMP(L), S(1,L), 1, S(1,I1), 1) 768 550 CONTINUE 769 560 CONTINUE 770 570 NSTART = MIN0(NSTART,NBLOCK) 771C 772C THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS. 773C 774 580 IF (NGOOD.EQ.0 .OR. .NOT.(TEST .OR. ENOUGH)) GO TO 600 775 DO 590 I=1,NGOOD 776 M = NPERM + I 777 RES(M) = ATEMP(I) 778 590 CONTINUE 779C 780C THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE 781C LANCZOS VECTORS. 782C 783 600 NUMBER = NPERM + NGOOD 784 IF (TEST .OR. ENOUGH) CALL DCOPY(N*NBLOCK, DZERO, 0, P1, 1) 785 IF (NGOOD.EQ.0) GO TO 620 786 M = NPERM + 1 787 DO 610 I=M,NUMBER 788 CALL DCOPY(N, DZERO, 0, VEC(1,I), 1) 789 610 CONTINUE 790 620 DO 670 I=NBLOCK,J,NBLOCK 791 CALL IOVECT(N, NBLOCK, P2, I, 1) 792 DO 660 K=1,NBLOCK 793 M = I - NBLOCK + K 794 IF (NSTART.EQ.0) GO TO 640 795 DO 630 L=1,NSTART 796 I1 = NGOOD + L 797 CALL DAXPY(N, S(M,I1), P2(1,K), 1, P1(1,L), 1) 798 630 CONTINUE 799 640 IF (NGOOD.EQ.0) GO TO 660 800 DO 650 L=1,NGOOD 801 I1 = L + NPERM 802 CALL DAXPY(N, S(M,L), P2(1,K), 1, VEC(1,I1), 1) 803 650 CONTINUE 804 660 CONTINUE 805 670 CONTINUE 806 IF (TEST .OR. ENOUGH) GO TO 690 807C 808C THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE 809C TAU RECURRENCE. 810C 811 M = NPERM + 1 812 DO 680 I=M,NUMBER 813 TEMP = 1.0D0/DNRM2(N,VEC(1,I),1) 814 CALL DSCAL(N, TEMP, VEC(1,I), 1) 815 TAU(I) = 1.0D0 816 OTAU(I) = 1.0D0 817 680 CONTINUE 818C 819C SHIFT S VECTORS TO ALIGN FOR LATER CALL TO DLAEIG 820C 821 CALL DCOPY(NTHETA, VAL(NPERM+1), 1, VTEMP, 1) 822 CALL DVSORT(NTHETA, VTEMP, ATEMP, 0, TARR, MAXJ, J, S) 823 GO TO 80 824C 825C ------------------------------------------------------------------ 826C 827C THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE 828C PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING 829C THE PERMANENT VECTORS. 830C 831 690 IF (NGOOD.EQ.0 .AND. NOP.GE.MAXOP) GO TO 810 832 IF (NGOOD.EQ.0) GO TO 30 833C 834C THIS ORTHONORMALIZES THE VECTORS 835C 836 CALL DORTQR(NMVEC, N, NPERM+NGOOD, VEC, S) 837C 838C THIS SORTS THE VALUES AND VECTORS. 839C 840 IF(NPERM .NE. 0) CALL DVSORT(NPERM+NGOOD, VAL, RES, 0, TEMP, 841 * NMVEC, N, VEC) 842 NPERM = NPERM + NGOOD 843 NLEFT = NLEFT - NGOOD 844 RNORM = DMAX1(-VAL(1),VAL(NPERM)) 845C 846C THIS DECIDES WHERE TO GO NEXT. 847C 848 IF (NOP.GE.MAXOP .AND. NLEFT.NE.0) GO TO 810 849 IF (NLEFT.NE.0) GO TO 30 850 IF (VAL(NVAL)-VAL(1).LT.TOLA) GO TO 790 851C 852C THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED 853C TO LOOK FOR UNDISCLOSED MULTIPLICITIES. 854C 855 M = NPERM - NBLOCK + 1 856 IF (M.LE.0) RETURN 857 DO 780 I=1,M 858 L = I + NBLOCK - 1 859 IF (VAL(L)-VAL(I).LT.TOLA) GO TO 30 860 780 CONTINUE 861C 862C THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ 863C PROCEDURE IS NEEDED. 864C 865 790 M = NPERM - NBLOCK 866 IF (M.LE.0) RETURN 867 DO 800 I=1,M 868 L = I + NBLOCK 869 IF (VAL(L)-VAL(I).GE.TOLA) GO TO 800 870 RARITZ = .TRUE. 871 RETURN 872 800 CONTINUE 873C 874 RETURN 875C 876C ------------------------------------------------------------------ 877C 878C THIS REPORTS THAT MAXOP WAS EXCEEDED. 879C 880 810 IERR = -2 881 GO TO 790 882C 883 END 884