C C $Id: diag.F,v 1.5 1998/07/16 16:39:34 jjv5 Exp arjan $ C C------------------------------------------------------------------------ #ifdef CRAY subroutine dummy end #else SUBROUTINE DIAG(NDIM,A,NEVEC1,TOLERA,V,EVAL1,IDEGEN1,EVEC1, & IERROR) C C DRIVER ROUTINE FOR DIAGONALIZATION OF THE REAL, SYMMETRIC, C MATRIX A. C C VARIABLES REQUIRED: C ------------------ C C NDIM = ORDER OF THE MATRIX A (I.E., A IS NDIM BY NDIM); C C A = REAL SYMMETRIC MATRIX TO BE DIAGONALIZED. ONLY THE LOWER C HALF OF A NEED BE FILLED WHEN CALLING. A IS DESTROYED BY C THIS ROUTINE. C C NEVEC1 = THE NUMBER OF EIGENVECTORS REQUIRED; C C TOLERA = TOLERANCE FACTOR USED IN THE QR ITERATION TO DETERMINE C WHEN OFF-DIAGONAL ENTRIES ARE ESSENTIALLY ZERO. (DEFAULT C IS 1.0D-8). C C V = 3 BY NDIM WORKSPACE. C C C VARIABLES RETURNED: C ------------------ C C EVAL1 = EIGENVALUES OF A (SORTED IN INCREASING ALGEBRAIC VALUE); C C IDEGEN1 = DEGENERACIES (NUMBER OF TIMES REPEATED) FOR EIGENVALUES; C C EVEC1 = EIGENVECTORS OF A (IN COLUMNS OF EVEC1); C C ERROR CODES: IERROR=0 - SUCCESSFUL CALCULATION. C IERROR=1 - NO CONVERGENCE IN QR ITERATION. C C C PROGRAMMED BY S. L. DIXON, OCT., 1991. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" DIMENSION A(NDIM,*),V(3,*),EVAL1(*),IDEGEN1(*),EVEC1(NDIM,*) C C FLAG FOR WHETHER OR NOT TO COMPUTE EIGENVALUES: C C IF(NDIM.EQ.1)THEN EVAL1(1) = A(1,1) IDEGEN1(1) = 1 EVEC1(1,1) = 1.0D0 RETURN ENDIF C C-RDC C TRIDIAGONALIZE THE MATRIX A. THIS WILL OVERWRITE THE DIAGONAL C AND SUBDIAGONAL OF A WITH THE TRIDIAGONALIZED VERSION. THE C HOUSEHOLDER VECTORS ARE RETURNED IN THE ROWS ABOVE THE DIAGONAL, C AND THE BETAHS ARE RETURNED BELOW THE SUBDIAGONAL. C CALL TRIDI(NDIM,V,A) C C COMPUTE NORM OF TRIDIAGONAL MATRIX FROM THE "LARGEST" COLUMN. C ANORM = ABS(A(1,1)) + ABS(A(2,1)) IF(NDIM.GT.2)THEN DO 20 I=2,NDIM-1 AICOL = ABS(A(I-1,I)) + ABS(A(I,I)) + ABS(A(I+1,I)) ANORM = MAX(ANORM,AICOL) 20 CONTINUE ENDIF ANCOL = ABS(A(NDIM-1,NDIM)) + ABS(A(NDIM,NDIM)) ANORM = MAX(ANORM,ANCOL) C C GET EIGENVALUES AND DEGENERACIES OF THE TRIDIAGONAL MATRIX A. C IF THE CALLING ROUTINE HAS NOT SUPPLIED A TOLERANCE FACTOR FOR C OFF-DIAGONAL ENTRIES IN THE QR ITERATION, A DEFAULT OF 1.0D-8 C WILL BE USED. C TOLTMP = TOLERA IF(TOLTMP.LE.0.0D0) TOLTMP = 1.0D-8 IF(DOEVAL) CALL EIGVAL(NDIM,A,V,TOLTMP,ANORM,EVAL1,IERROR) IF(IERROR.NE.0) RETURN C C DETERMINE DEGENERACIES OF EIGENVALUES. C CALL DEGEN(NDIM,EVAL1,TOLTMP,ANORM,IDEGEN1) C C GET EIGENVECTORS OF TRIDIAGONALIZED VERSION OF A. C IF(NEVEC1.LE.0) RETURN CALL EIGVEC(NDIM,NEVEC1,A,V,TOLTMP,ANORM,EVAL1,IDEGEN1,EVEC1) C C PREMULTIPLY EVEC1 BY THE HOUSEHOLDER MATRIX USED TO TRIDIAGONALIZE C A. THIS TRANSFORMS EIGENVECTORS OF THE TRIDIAGONAL MATRIX TO C THOSE OF THE ORIGINAL MATRIX A. SEE SUBROUTINE TRIDI FOR C STORAGE OF HOUSEHOLDER TRANSFORMATION. C IF(NDIM.GT.2)THEN DO 30 K=1,NDIM-2 V(1,K) = A(K+2,K) 30 CONTINUE C C SWAP STORAGE SO THAT THE EXPENSIVE TRIPLE LOOP BELOW DOESN'T C HAVE TO JUMP ACROSS COLUMNS TO GET ENTRIES OF A. C DO 50 I=2,NDIM DO 40 J=1,I-1 A(I,J) = A(J,I) 40 CONTINUE 50 CONTINUE DO 160 J=1,NEVEC1 DO 140 M=2,NDIM-1 K = NDIM - M BETAH = V(1,K) IF(ABS(BETAH).LT.1.0D-50) GO TO 140 SUM = 0.0D0 DO 100 I=K+1,NDIM SUM = SUM + A(I,K)*EVEC1(I,J) 100 CONTINUE BSUM = BETAH*SUM DO 120 I=K+1,NDIM EVEC1(I,J) = EVEC1(I,J) - A(I,K)*BSUM 120 CONTINUE 140 CONTINUE 160 CONTINUE ENDIF RETURN END C C C SUBROUTINE TRIDI(NDIM,V,A) C C TRIDIAGONALIZES A REAL, SYMMETRIC MATRIX A BY THE METHOD OF C HOUSEHOLDER (J. H. WILKINSON, THE COMPUTER JOURNAL, VOL. 3, C P. 23 (1960)). NDIM IS THE ORDER OF A. THE DIAGONAL AND C SUBDIAGONAL OF A ARE OVERWRITTEN WITH THE TRIDIAGONALIZED C VERSION OF A. THE VECTORS USED IN EACH HOUSEHOLDER C TRANSFORMATION ARE STORED ABOVE THE DIAGONAL IN THE FIRST C NDIM-2 ROWS OF A. THE BETAHS ARE RETURNED BELOW THE SUBDIAGONAL C OF A. V IS A WORKSPACE ARRAY. C C PROGRAMMED BY S. L. DIXON, OCT., 1991. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NDIM,*),V(3,*) C C THRESH WILL BE USED AS A THRESHOLD TO DETERMINE IF A VALUE SHOULD C BE CONSIDERED TO BE ZERO. THIS CAN BE CHANGED BY THE USER. C THRESH = 1.0D-50 C C IF A IS 2 BY 2 OR SMALLER, THEN IT IS ALREADY TRIDIAGONAL -- NO C NEED TO CONTINUE. C IF(NDIM.LE.2) GO TO 1000 DO 500 K=1,NDIM-2 C C DETERMINE THE VECTOR V USED IN THE HOUSEHOLDER TRANSFORMATION P. C FOR EACH VALUE OF K THE HOUSEHOLDER MATRIX P IS DEFINED AS: C C P = I - BETAH*V*V' C C C CONSTRUCT A HOUSEHOLDER TRANSFORMATION ONLY IF THERE IS A NONZERO C OFF-DIAGONAL ELEMENT BELOW A(K,K). C ALPHA2 = 0.0D0 DO 60 I=K+1,NDIM V(1,I) = A(I,K) ALPHA2 = ALPHA2 + V(1,I)**2 60 CONTINUE APTEMP = ALPHA2 - V(1,K+1)**2 ALPHA = DSQRT(ALPHA2) IF(ALPHA.GE.THRESH)THEN BETAH = 1.0D0/(ALPHA*(ALPHA + ABS(V(1,K+1)))) SGN = SIGN(1.0D0,V(1,K+1)) V(1,K+1) = V(1,K+1) + SGN*ALPHA C C-RDC C NOW OVERWRITE A WITH P'*A*P. THE ENTRIES BELOW THE SUBDIAGONAL C IN THE KTH COLUMN ARE ZEROED BY THE PREMULTIPLICATION BY P'. C THESE ENTRIES WILL BE LEFT ALONE TO SAVE TIME. C AKV = APTEMP + A(K+1,K)*V(1,K+1) S = BETAH*AKV A(K+1,K) = A(K+1,K) - S*V(1,K+1) C C NOW THE SUBMATRIX CONSISTING OF ROWS K+1,NDIM AND COLUMNS K+1,NDIM C MUST BE OVERWRITTEN WITH THE TRANSFORMATION. C DOT12 = 0.0D0 BHALF = BETAH*0.5D0 DO 220 I=K+1,NDIM SUM = 0.0D0 DO 100 J=K+1,I SUM = SUM + A(I,J)*V(1,J) 100 CONTINUE IF(I.LT.NDIM)THEN DO 180 J=I+1,NDIM C C AN UPPER TRIANGULAR ENTRY OF A WILL BE REQUIRED. MUST USE C THE SYMMETRIC ENTRY IN THE LOWER TRIANGULAR PART OF A. C SUM = SUM + A(J,I)*V(1,J) 180 CONTINUE ENDIF V(2,I) = BETAH*SUM DOT12 = DOT12 + V(1,I)*V(2,I) 220 CONTINUE BH12 = BHALF*DOT12 DO 300 I=K+1,NDIM V(2,I) = V(2,I) - BH12*V(1,I) 300 CONTINUE DO 350 J=K+1,NDIM DO 310 I=J,NDIM A(I,J) = A(I,J) - V(1,I)*V(2,J) - V(2,I)*V(1,J) 310 CONTINUE C C STORE V(1,J) ABOVE THE DIAGONAL IN ROW K OF A C A(K,J) = V(1,J) 350 CONTINUE C C STORE BETAH BELOW THE SUBDIAGONAL OF A. C A(K+2,K) = BETAH ELSE C C NO HOUSEHOLDER TRANSFORMATION IS NECESSARY BECAUSE THE OFF- C DIAGONALS ARE ALL ESSENTIALLY ZERO. C A(K+2,K) = 0.0D0 DO 460 J=K+1,NDIM A(K,J) = 0.0D0 460 CONTINUE ENDIF 500 CONTINUE 1000 RETURN END C C C SUBROUTINE EIGVAL(NDIM,A,BETAH,TOLERA,ANORM,EVAL1,IERROR) C C QR ROUTINE FOR THE DETERMINATION OF ALL THE EIGENVALUES C OF THE NDIM BY NDIM SYMMETRIC, TRIDIAGONAL MATRIX A. C C INPUT: C C NDIM = SIZE OF MATRIX A. C A = NDIM BY NDIM SYMMETRIC TRIDIAGONAL MATRIX. C BETAH = 3 BY NDIM WORKSPACE. C TOLERA = SMALL NUMBER USED TO DETERMINE WHEN OFF-DIAGONAL C ELEMENTS ARE ESSENTIALLY ZERO. C ANORM = ABSOLUTE COLUMN NORM OF TRIDIAGONAL MATRIX A. C C C RETURNED: C C EVAL1 = EIGENVALUES OF A IN ASCENDING ORDER. C IERROR = 1 IF QR ITERATION DID NOT CONVERGE; 0 OTHERWISE. C C PROGRAMMED BY S. L. DIXON. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NDIM,*),BETAH(3,*),EVAL1(*) IERROR = 0 ITMAX = 20 C C TOLERANCE FOR OFF-DIAGONAL ELEMENTS: C EPSLON = TOLERA*ANORM C C COPY DIAGONAL ELEMENTS OF A TO EVAL1, AND SUBDIAGONAL ELEMENTS C TO BETAH. C EVAL1(1) = A(1,1) BETAH(1,1) = A(2,1) IF(NDIM.GT.2)THEN DO 50 I=2,NDIM-1 EVAL1(I) = A(I,I) BETAH(1,I) = A(I+1,I) 50 CONTINUE ENDIF EVAL1(NDIM) = A(NDIM,NDIM) C C EACH QR ITERATION WILL OPERATE ON THE UNREDUCED TRIDIAGONAL C SUBMATRIX WITH UPPER LEFT ELEMENT (L,L) AND LOWER RIGHT ELEMENT C (N,N). C L = 1 N = NDIM ITER = 0 C C FIND THE SMALLEST UNREDUCED SUBMATRIX WITH LOWER RIGHT CORNER AT C (N,N). I.E., SEARCH UPWARD FOR A BETAH THAT IS ZERO. C 80 KUPPER = N-L DO 100 K=1,KUPPER I = N-K IF(ABS(BETAH(1,I)).LE.EPSLON)THEN L = I+1 GO TO 150 ENDIF 100 CONTINUE C C IF WE GET TO THE NEXT STATEMENT, THEN THERE ARE NO ZERO OFF-DIAGONALS C FOR THE SUBMATRIX WITH UPPER LEFT A(L,L) AND LOWER RIGHT A(N,N). C WE CAN STILL GET EIGENVALUES IF THE MATRIX IS 2 BY 2 OR 1 BY 1. C OTHERWISE, DO ANOTHER QR ITERATION PROVIDED ITMAX CYCLES HAVE C NOT OCCURRED. C IF(L.EQ.N.OR.L.EQ.N-1)THEN GO TO 150 ELSE IF(ITER.EQ.ITMAX)THEN IERROR = 1 GO TO 1000 ELSE GO TO 200 ENDIF ENDIF C C IF WE GET TO 150 THEN A(L,L-1) IS ZERO AND THE UNREDUCED SUBMATRIX C HAS UPPER LEFT AT A(L,L) AND LOWER RIGHT AT A(N,N). WE CAN C EXTRACT ONE EIGENVALUE IF THIS MATRIX IS 1 BY 1 AND 2 EIGENVALUES C IF IT IS 2 BY 2. C 150 IF(L.EQ.N)THEN C C IT'S A 1 BY 1 AND EVAL1(N) IS AN EIGENVALUE. IF L=2 OR 1 WE ARE C DONE. OTHERWISE, UPDATE N, RESET L AND ITER, AND REPEAT THE C SEARCH. C IF(L.LE.2)THEN GO TO 500 ELSE N = L-1 L = 1 ITER = 0 GO TO 80 ENDIF ELSEIF(L.EQ.N-1)THEN C C-RDC C THE UNREDUCED SUBMATRIX IS A 2 BY 2. OVERWRITE EVAL1(N-1) C AND EVAL1(N) WITH THE EIGENVALUES OF THE LOWER RIGHT 2 BY 2. C BTERM = EVAL1(N-1) + EVAL1(N) ROOT1 = BTERM*0.5D0 ROOT2 = ROOT1 DISCR = BTERM**2 - 4.0D0*(EVAL1(N-1)*EVAL1(N)-BETAH(1,N-1)**2) IF(DISCR.GT.0.0D0)THEN D = DSQRT(DISCR)*0.5D0 ROOT1 = ROOT1 - D ROOT2 = ROOT2 + D ENDIF EVAL1(N-1) = ROOT1 EVAL1(N) = ROOT2 C C SEE IF WE ARE DONE. IF NOT, RESET N, L, AND ITER AND LOOK C FOR NEXT UNREDUCED SUBMATRIX. C IF(L.LE.2)THEN GO TO 500 ELSE N = L-1 L = 1 ITER = 0 GO TO 80 ENDIF ELSE C C AN EIGENVALUE WAS FOUND AND THE NEW UNREDUCED MATRIX LIMITS C N AND L ARE SET. DO A QR ITERATION ON NEW MATRIX. C ITER = 0 GO TO 200 ENDIF C C QR ITERATION BEGINS HERE. C 200 ITER = ITER + 1 C C USE EIGENVALUES OF THE LOWER RIGHT 2 BY 2 TO COMPUTE SHIFT. SHIFT C BY THE EIGENVALUE CLOSEST TO EVAL1(N). C D = (EVAL1(N-1) - EVAL1(N))*0.5D0 SIGND = 1.0D0 IF(D.LT.0.0D0) SIGND = -1.0D0 SHIFT = EVAL1(N) + D - SIGND*DSQRT(D*D + BETAH(1,N-1)**2) P = EVAL1(L) - SHIFT R = BETAH(1,L) T = EVAL1(L) W = BETAH(1,L) C C-RDC C OVERWRITE A WITH Q'*A*Q. C DO 250 K=L,N-1 D = DSQRT(P*P + R*R) C = P/D S = R/D IF(K.NE.L) BETAH(1,K-1) = D CC = C*C SS = 1.0D0 - CC CS = C*S CSW = 2.0D0*CS*W AK1 = EVAL1(K+1) EVAL1(K) = CC*T + CSW + SS*AK1 P = (CC - SS)*W + CS*(AK1 - T) T = SS*T - CSW + CC*AK1 R = S*BETAH(1,K+1) W = C*BETAH(1,K+1) 250 CONTINUE BETAH(1,N-1) = P EVAL1(N) = T C C GO BACK AND SEE IF L AND N NEED TO BE UPDATED. C GO TO 80 C C SORT EIGENVALUES IN ASCENDING ALGEBRAIC ORDER. C 500 DO 600 I=2,NDIM JMAX = NDIM-I+1 ISORT = 0 DO 550 J=1,JMAX IF(EVAL1(J).GT.EVAL1(J+1))THEN ETEMP = EVAL1(J) EVAL1(J) = EVAL1(J+1) EVAL1(J+1) = ETEMP ISORT = 1 ENDIF 550 CONTINUE IF(ISORT.EQ.0) GO TO 1000 600 CONTINUE 1000 RETURN END C C C SUBROUTINE DEGEN(NDIM,EVAL1,TOLERA,ANORM,IDEGEN1) C C DETERMINES DEGENERACIES OF THE EIGENVALUES. C C INPUT: C C NDIM = SIZE OF MATRIX BEING DIAGONALIZED. C EVAL1 = SORTED EIGENVALUES (INCREASING VALUE). C TOLERA = SAME TOLERANCE USED TO DETERMINE EIGENVALUES. C ANORM = ABSOLUTE COLUMN NORM OF TRIDIAGONAL MATRIX. C C C RETURNED: C C IDEGEN1 = DEGENERACIES OF EIGENVALUES. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION EVAL1(*),IDEGEN1(*) C C DETERMINE DEGENERACIES OF EIGENVALUES. ADJACENT EIGENVALUES C WILL BE CONSIDERED TO BE DEGENERATE WHEN THEY DIFFER BY LESS C THAN DTOLER. C DTOLER = MAX(ANORM*DSQRT(TOLERA),1.0D-8) NSAME = 1 DO 200 I=2,NDIM DIFF = ABS(EVAL1(I-1) - EVAL1(I)) IF(DIFF.LE.DTOLER)THEN C C EIGENVALUES I-1 AND I ARE DEGENERATE. C NSAME = NSAME + 1 IF(I.EQ.NDIM)THEN C C WE'VE COME TO THE LAST REQUESTED EIGENVALUE, AND IT'S TIME C TO ASSIGN DEGENERACIES FOR THE BLOCK ENDING WITH THE ITH C EIGENVALUE. C DO 100 J=I-NSAME+1,I IDEGEN1(J) = NSAME 100 CONTINUE ENDIF C C GO TO THE NEXT EIGENVALUE (IF THERE ARE ANY LEFT) AND SEE IF C IT'S DEGENERATE WITH THE NSAME EIGENVALUES WE'VE ALREADY C FOUND. C GO TO 200 ELSE C C EITHER EIGENVALUE I-1 IS NONDEGENERATE OR IT'S THE LAST C EIGENVALUE IN A DEGENERATE BLOCK. CORRESPONDINGLY, ASSIGN THE C PROPER DEGENERACY TO I-1 OR TO EACH EIGENVALUE IN THE BLOCK. C DO 150 J=I-NSAME,I-1 IDEGEN1(J) = NSAME 150 CONTINUE NSAME = 1 C C IF I=NDIM THEN IT MUST BE THE CASE THAT THIS LAST EIGENVALUE C IS NONDEGENERATE. C IF(I.EQ.NDIM) IDEGEN1(I) = 1 ENDIF 200 CONTINUE RETURN END C C C SUBROUTINE EIGVEC(NDIM,NEVEC1,A,AWORK,TOLERA,ANORM,EVAL1,IDEGEN1, . EVEC1) C C INVERSE ITERATION ROUTINE FOR EIGENVECTOR DETERMINATION. C CALCULATES THE EIGENVECTORS OF AN NDIM BY NDIM SYMMETRIC, C TRIDIAGONAL MATRIX A. C C INPUT: C C NDIM = SIZE OF MATRIX A. C NEVEC1 = NUMBER OF EIGENVECTORS REQUIRED. C A = NDIM BY NDIM TRIDIAGONAL MATRIX. C TOLERA = SAME TOLERANCE USED TO DETERMINE EIGENVALUES. C ANORM = ABSOLUTE COLUMN NORM OF TRIDIAGONAL MATRIX A. C EVAL1 = SORTED EIGENVALUES OF A. C IDEGEN1 = DEGENERACIES OF EIGENVALUES. C C C RETURNED: C C EVEC1 = EIGENVECTORS OF TRIDIAGONAL MATRIX (IN COLUMNS). C C PROGRAMMED BY S. L. DIXON, OCT., 1991. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NDIM,*),AWORK(3,*),EVAL1(*),IDEGEN1(*),EVEC1(NDIM,*) LOGICAL ORTH IRAND = 13876532 C C COMPUTE THRESHOLD EPSLON WHICH WILL BE USED IF THE INVERSE ITERATION C MATRIX IS SINGULAR. C EPSLON = ANORM*TOLERA C C WHEN DEGENERACIES OCCUR, THERE ARE RARE INSTANCES WHEN THE C DEGENERATE BLOCK OF EIGENVECTORS ARE NOT LINEARLY INDEPENDENT. C IN THESE CASES, AN ADDITIONAL PASS THROUGH THE INVERSE ITERATION C (WITH A NEW SET OF RANDOM NUMBERS) IS CARRIED OUT, I.E., CONTROL C PASSES TO STATEMENT 40. IVECT WILL KEEP TRACK OF THE CURRENT C STARTING EIGENVECTOR WHEN ADDITIONAL PASSES ARE NECESSARY. C IVECT = 1 NFAIL = 0 NPRTRB = 0 C C DO ONE ITERATION FOR EACH EIGENVECTOR. C 40 ORTH = .TRUE. NDEGEN = 0 MSTART = IVECT DO 380 M=MSTART,NEVEC1 IF(IDEGEN1(M).GT.1) NDEGEN = NDEGEN + 1 Z = EVAL1(M) C C IF THE INVERSE ITERATION HAS FAILED TWICE DUE TO NON-ORTHOGONALITY C OF DEGENERATE EIGENVECTORS, PERTURB THE EIGENVALUE BY A SMALL C AMOUNT. C IF(NFAIL.GE.2)THEN NPRTRB = NPRTRB + 1 Z = Z + 0.001D0*DFLOAT(NPRTRB)*TOLERA ENDIF C C STORE THE TRIDIAGONAL ENTRIES OF THE INVERSE ITERATION MATRIX IN C THE 3 BY NDIM WORKSPACE AWORK. C AWORK(1,1) = 0.0D0 AWORK(2,1) = A(1,1) - Z AWORK(3,1) = A(2,1) IF(NDIM.GT.2)THEN DO 80 I=2,NDIM-1 AWORK(1,I) = A(I,I-1) AWORK(2,I) = A(I,I) - Z AWORK(3,I) = A(I+1,I) 80 CONTINUE ENDIF AWORK(1,NDIM) = A(NDIM,NDIM-1) AWORK(2,NDIM) = A(NDIM,NDIM) - Z AWORK(3,NDIM) = 0.0D0 C C ASSIGN INVERSE ITERATION VECTOR FROM RANDOM NUMBERS. C DO 120 I=1,NDIM CALL RANDOM1(IRAND,RNDOM) RNDOM = 2.0D0*(RNDOM - 0.5D0) EVEC1(I,M) = RNDOM*TOLERA 120 CONTINUE C C CARRY OUT FORWARD GAUSSIAN ELIMINATION WITH ROW PIVOTING C ON THE INVERSE ITERATION MATRIX. C DO 160 K=1,NDIM-1 ADIAG = ABS(AWORK(2,K)) ASUB = ABS(AWORK(1,K+1)) IF(ADIAG.GE.ASUB)THEN C C USE PIVOTAL ELEMENT FROM ROW K. C IF(AWORK(2,K).EQ.0.0D0) AWORK(2,K) = EPSLON T = AWORK(1,K+1)/AWORK(2,K) AWORK(1,K+1) = 0.0D0 AWORK(2,K+1) = AWORK(2,K+1) - T*AWORK(3,K) C C LEFT-JUSTIFY EQUATION K SO THAT DIAGONAL ENTRY IS STORED C IN AWORK(1,K). C AWORK(1,K) = AWORK(2,K) AWORK(2,K) = AWORK(3,K) AWORK(3,K) = 0.0D0 C C OPERATE ON VECTOR AS WELL. C EVEC1(K+1,M) = EVEC1(K+1,M) - T*EVEC1(K,M) ELSE C C USE PIVOTAL ELEMENT FROM ROW K+1 AND SWAP ROWS K AND K+1. C IF(AWORK(1,K+1).EQ.0.0D0) AWORK(1,K+1) = EPSLON T = AWORK(2,K)/AWORK(1,K+1) ATEMP = AWORK(3,K) - T*AWORK(2,K+1) AWORK(1,K) = AWORK(1,K+1) AWORK(2,K) = AWORK(2,K+1) AWORK(3,K) = AWORK(3,K+1) AWORK(1,K+1) = 0.0D0 AWORK(2,K+1) = ATEMP AWORK(3,K+1) = -T*AWORK(3,K+1) C C OPERATE ON VECTOR AND SWAP ENTRIES. C ETEMP = EVEC1(K+1,M) EVEC1(K+1,M) = EVEC1(K,M) - ETEMP*T EVEC1(K,M) = ETEMP ENDIF 160 CONTINUE C C FORWARD ELIMINATION COMPLETE. BACK SUBSTITUTE TO GET SOLUTION. C-RDC C OVERWRITE COLUMN M OF EVEC1 WITH SOLUTION. C IF(AWORK(2,NDIM).EQ.0.0D0) AWORK(2,NDIM) = EPSLON EVEC1(NDIM,M) = EVEC1(NDIM,M)/AWORK(2,NDIM) ETEMP = EVEC1(NDIM-1,M) - AWORK(2,NDIM-1)*EVEC1(NDIM,M) EVEC1(NDIM-1,M) = ETEMP/AWORK(1,NDIM-1) ENORM = EVEC1(NDIM,M)**2 + EVEC1(NDIM-1,M)**2 IF(NDIM.GT.2)THEN C C CAUTION: PROBLEM LOOP FOR SOME IBM RS/6000 COMPILERS. VALUE C OF K CAN GET LOST WHEN OPTIMIZE FLAG IS USED. C DO 200 L=1,NDIM-2 K = NDIM-L-1 ETEMP = EVEC1(K,M) - AWORK(2,K)*EVEC1(K+1,M) . - AWORK(3,K)*EVEC1(K+2,M) EVEC1(K,M) = ETEMP/AWORK(1,K) ENORM = ENORM + EVEC1(K,M)**2 200 CONTINUE ENDIF EINV = 1.0D0/DSQRT(ENORM) C C NORMALIZE EIGENVECTOR. C DO 240 I=1,NDIM EVEC1(I,M) = EVEC1(I,M)*EINV 240 CONTINUE C C IF WE HAVE COME TO THE END OF A DEGENERATE BLOCK OF EIGENVECTORS, C ORTHOGONALIZE THE BLOCK. C IF(NDEGEN.GT.1)THEN IF(NDEGEN.EQ.IDEGEN1(M).OR.M.EQ.NEVEC1)THEN JSTART = M-NDEGEN+1 CALL ORTHOG1(NDIM,NDEGEN,JSTART,EVEC1,ORTH) IF(ORTH)THEN NFAIL = 0 NPRTRB = 0 C C THE DEGENERATE VECTORS WERE LINEARLY INDEPENDENT AND WERE C SUCCESSFULLY ORTHOGONALIZED. C IVECT = IVECT + NDEGEN NDEGEN = 0 ELSE C C THE BLOCK IS APPARENTLY NOT LINEARLY INDEPENDENT. GO BACK C AND REPEAT THE INVERSE ITERATION FOR THESE VECTORS. AFTER C AN INDEPENDENT SET HAS BEEN FOUND, ANY ADDITIONAL EIGENVECTORS C WILL BE DETERMINED. C NFAIL = NFAIL + 1 GO TO 40 ENDIF ENDIF ENDIF C C THE CURRENT EIGENVECTOR SHOULD BE OKAY IF IT IS NONDEGENERATE. C IF(IDEGEN1(M).EQ.1) IVECT = IVECT + 1 380 CONTINUE RETURN END C C C SUBROUTINE ORTHOG1(NDIM,NVECT,JSTART,VECT,ORTH) C C CONSTRUCTS A SET OF ORTHONORMAL VECTORS FROM THE NVECT LINEARLY C INDEPENDENT, NORMALIZED VECTORS IN THE ARRAY VECT. THE VECTORS C SHOULD BE STORED COLUMNWISE, STARTING IN COLUMN JSTART. VECT IS C OVERWRITTEN WITH THE ORTHONORMAL SET. ALL VECTORS ARE NDIM BY 1. C ORTH IS RETURNED WITH A VALUE OF .TRUE. IF THE SET WAS LINEARLY C INDEPENDENT AND .FALSE. OTHERWISE. C C PROGRAMMED BY S. L. DIXON. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION VECT(NDIM,*) LOGICAL ORTH C ORTH = .TRUE. ORTEST = 1.0D-8 C C BEGIN ORTHOGONALIZATION. C JSTOP = JSTART + NVECT - 1 DO 120 J=JSTART,JSTOP IF(J.GT.JSTART)THEN C C SUBTRACT OFF COMPONENTS OF PREVIOUSLY DETERMINED ORTHOGONAL C VECTORS FROM THE VECTOR IN COLUMN J. C DO 60 JPREV=JSTART,J-1 DOT = 0.0D0 DO 20 I=1,NDIM DOT = DOT + VECT(I,JPREV)*VECT(I,J) 20 CONTINUE DO 40 I=1,NDIM VECT(I,J) = VECT(I,J) - DOT*VECT(I,JPREV) 40 CONTINUE 60 CONTINUE ENDIF C C NORMALIZE COLUMN J. C VJNORM = 0.0D0 DO 80 I=1,NDIM VJNORM = VJNORM + VECT(I,J)**2 80 CONTINUE VJNORM = DSQRT(VJNORM) C C IF THE NORM OF THIS VECTOR IS TOO SMALL THEN THE VECTORS ARE C NOT LINEARLY INDEPENDENT. C IF(VJNORM.LT.ORTEST)THEN ORTH = .FALSE. GO TO 1000 ENDIF DO 100 I=1,NDIM VECT(I,J) = VECT(I,J)/VJNORM 100 CONTINUE 120 CONTINUE 1000 RETURN END C C C SUBROUTINE RANDOM1(IX,X) C C GENERATES INTEGER*4 AND DOUBLE PRECISION RANDOM NUMBERS ON C THE INTERVALS: C C 0 < IX < (2**31)-1 C AND 0.0 < X < 1.0 C C ON THE FIRST CALL IX SHOULD SATISFY THE TOP INEQUALITY. C C NUMBERS ARE GENERATED USING THE RELATION, C C IX = IX*IC (MODULO (2**31)-1), WHERE IC = 7**5 C C DOUBLE PRECISION X C C C INITIALIZE: I15 = 2**15 C I16 = 2**16 C I31_1 = 2**31 - 1 C IC = 7**5 C DATA I15 /32768/ DATA I16 /65536/ DATA I31_1 /2147483647/ DATA IC /16807/ C SAVE I15,I16,I31_1,IC C IX16 = IX/I16 I16RMD = IX - IX16*I16 C C NOTE THAT IX = IX16*I16 + I16RMD (I16RMD = 16-BIT REMAINDER) C IX16IC = IX16*IC IXIC31 = IX16IC/I15 C C NOTE: IX*IC = (IX16*I16 + I16RMD)*IC C = (IX16*IC*I16) + (I16RMD*IC) C = (IX16IC*I16 ) + (I16RMD*IC) C = ( TERM 1 ) + ( TERM 2 ) C AND, C C IX16IC = ((IX16IC/I15)*I15) + (IX16IC - (IX16IC/I15)*I15)) C = ( IXIC31 )*I15) + (IX16IC - ( IXIC31 )*I15 ) C ( 15-BIT REMAINDER ) C C THEREFORE, TERM 1 = ((IXIC31*I15) + (IX16IC - IXIC31*I15))*I16 C C THEN, C C ( TERM A ) ( TERM B ) ( TERM C ) C IX*IC = ((IXIC31*I16*I15) + (IX16IC - IXIC31*I15)*I16) + (I16RMD*IC) C = ( TERM 1 ) + ( TERM 2 ) C C C NOTE THAT TERM B AND TERM C ARE BOTH LESS THAN 2**31 - 1. ONLY C TERM A HAS THE POSSIBILITY OF EXCEEDING 2**31 - 1. BUT SINCE C I16*I15 = 2**31, THE FACTOR IXIC31 INDICATES EXACTLY HOW MANY TIMES C TERM A "WRAPS" AROUND THE INTERVAL (0,2**31 - 1). THUS, FOR THE C MODULO OPERATION, TERM A MAY BE REPLACED BY IXIC31. THE SUM OF C TERM A AND TERM B MIGHT EXCEED 2**31 - 1, BUT WE CAN SUBSTRACT C 2**31 - 1 FROM ONE OF THEM TO PREVENT THIS FROM HAPPENING. C IX = IXIC31 + ((IX16IC-IXIC31*I15)*I16 - I31_1) + I16RMD*IC C C ADD I31_1 BACK IN IF THE SUBTRACTION MADE IX NEGATIVE. C IF(IX.LT.0) IX = IX + I31_1 C C MAKE X RANDOM ON (0.0,1.0) BY MULTIPLYING IX BY 1.0/I31_1 C X = DFLOAT(IX)*4.6566128752458D-10 RETURN END #endif