1CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2C this subroutine is moved to the top of the 3C source file in order to resolve the mis- 4C matching parameter types later in code. 5C 2005-10-11 Tommi Hassinen 6CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 7 SUBROUTINE EC08C(A,B,VALUE,VEC,N,IV,W) 8C 9C TO FIND THE EIGENVALUES AND VECTORS OF A TRI-DIAGONAL 10C HERMITIAN MATRIX. 11 REAL VALUE(*),W(*),PCK(2),ONE,ZERO,VEC(*) 12 COMPLEX A(*),B(*),DN,UPCK 13 EQUIVALENCE (PCK(1),UPCK) 14C WE TREAT VEC AS IF IT IS DEFINED AS COMPLEX VEC(IV,N) 15C IN THE CALLING PROGRAM. 16 DATA ONE, ZERO/1.0,0.0/ 17 IV2=IV+IV 18 N2=N+N 19 IL=IV2*(N-1)+1 20 W(1)=A(1) 21C 22C THE HERMITIAN FORM IS TRANSFORMED INTO A REAL FORM 23 IF(N-1)80,80,10 24 10 DO 20 I=2,N 25 W(I)=A(I) 26 20 W(I+N)=CABS(B(I)) 27C 28C FIND THE EIGENVALUES AND VECTORS OF THE REAL FORM 29 30 CALL EA08C(W,W(N+1),VALUE,VEC,N,IV2,W(N2+1)) 30C 31C THE VECTORS IN VEC AT THIS POINT ARE REAL,WE NOW EXPAND THEM 32C INTO VEC AS THOUGH THEY WERE COMPLEX. 33 DO 50 I=1,IL,IV2 34 DO 40 J=1,N 35 K=N-J 36 L=K+K 37 40 VEC(I+L)=VEC(I+K) 38 50 VEC(I+1)=ZERO 39 IF(N.LE.1)GO TO 80 40 DN=ONE 41 K=1 42C 43C TRANSFORM VECTORS OF REAL FORM TO THOSE OF COMPLEX FORM. 44 DO 70 I=4,N2,2 45 K=K+1 46 UPCK=ONE 47 IF(W(K+N).NE.ZERO)UPCK=DN*CONJG(B(K))/W(K+N) 48 I1=IL-1+I 49 DO 60 J=I,I1,IV2 50 VEC(J)=VEC(J-1)*PCK(2) 51 60 VEC(J-1)=VEC(J-1)*PCK(1) 52 70 DN=UPCK 53 80 RETURN 54 END 55CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 56C this subroutine is moved to the top of the 57C source file in order to resolve the mis- 58C matching parameter types later in code. 59C 2005-10-11 Tommi Hassinen 60CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 61 SUBROUTINE ME08B (A,Q,B,N,IA) 62 REAL A(IA,*),Q(2,*),B(IA,*) 63 DO 10 I=1,N 64 A(1,I)=A(1,I) -Q(1,1)*B(1,I)+Q(2,1)*B(2,I) 65 1 -Q(1,I)*B(1,1)+Q(2,I)*B(2,1) 66 10 A(2,I)=A(2,I)-Q(2,1)*B(1,I)-Q(1,1)*B(2,I) 67 1 +Q(2,I)*B(1,1)+Q(1,I)*B(2,1) 68 RETURN 69 END 70 SUBROUTINE CDIAG(A,VALUE,VEC,N, NEED) 71C 72C TO FIND THE EIGENVALUES AND EIGENVECTORS OF A HERMITIAN MATRIX. 73 REAL VALUE(*),H 74 COMPLEX W(6000) 75 COMPLEX A(N,*),VEC(N,*) 76 COMPLEX FM06AS 77 IA=N 78 IV=N 79C 80C REDUCE MATRIX TO A TRI-DIAGONAL HERMITIAN MATRIX. 81 CALL ME08A(A,W,W(N+1),N,IA,W(2*N+1)) 82C 83C FIND THE EIGENVALUES AND EIGENVECTORS OF THE TRI-DIAGONAL MATRIX 84 CALL EC08C(W,W(N+1),VALUE,VEC,N,IV,W(2*N+1)) 85 IF(NEED.EQ.0) GOTO 50 86 IF(N.LT.2)RETURN 87C 88C THE EIGENVECTORS OF THE ORIGINAL MATRIX ARE NOW FOUND BY 89C BACK TRANSFORMATION USING INFORMATION STORE IN THE UPPER 90C TRIANGLE OF MATRIX A (BY ME08) 91 DO 40 II=3,N 92 I=N-II+1 93 H=W(N+I+1)*CONJG(A(I,I+1)) 94 IF(H)10,40,10 95 10 DO30L=1,N 96 I1=I+1 97 S=FM06AS(N-I,A(I,I+1),IA,VEC(I+1,L),1) 98 S=S/H 99 DO20K=I1,N 100 20 VEC(K,L)=VEC(K,L)+CONJG(A(I,K))*S 101 30 CONTINUE 102 40 CONTINUE 103 50 CALL SORT(VALUE,VEC,N) 104 RETURN 105 END 106 SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W) 107C STANDARD FORTRAN 66 (A VERIFIED PFORT SUBROUTINE) 108 DIMENSION A(*),B(*),VALUE(*),VEC(*),W(*) 109 DATA EPS/1.E-6/,A34/0.0/ 110C THIS USES QR ITERATION TO FIND THE EIGENVALUES AND EIGENVECTORS 111C OF THE SYMMETRIC TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE 112C A(I),I=1,M AND OFF-DIAGONAL ELEMENTS ARE B(I),I=2,M. THE ARRAY 113C W IS USED FOR WORKSPACE AND MUST HAVE DIMENSION AT LEAST 2*M. 114C WE TREAT VEC AS IF IT HAD DIMENSIONS (IV,M). 115 SML=EPS*FLOAT(M) 116 CALL EA09C(A,B,W(M+1),M,W) 117C SET VEC TO THE IDENTITY MATRIX. 118 DO 20 I=1,M 119 VALUE(I)=A(I) 120 W(I)=B(I) 121 K=(I-1)*IV+1 122 L=K+M-1 123 DO 10 J=K,L 124 10 VEC(J)=0. 125 KI=K+I-1 126 20 VEC(KI)=1. 127 K=0 128 IF(M.EQ.1)RETURN 129 DO 100 N3=2,M 130 N2=M+2-N3 131C EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2 132 MN2=M+N2 133 ROOT=W(MN2) 134 DO 80 ITER=1,20 135 BB=(VALUE(N2)-VALUE(N2-1))*0.5 136 CC=W(N2)*W(N2) 137 A22=VALUE(N2) 138 IF(CC.NE.0.0)A22=A22+CC/(BB+SIGN(1.0,BB)*SQRT(BB*BB+CC)) 139 DO 30 I=1,N2 140 MI=M+I 141 IF(ABS(ROOT-A22).LE.ABS(W(MI)-A22))GO TO 30 142 ROOT=W(MI) 143 MN=M+N2 144 W(MI)=W(MN) 145 W(MN)=ROOT 146 30 CONTINUE 147 DO 40 II=2,N2 148 N1=2+N2-II 149 IF(ABS(W(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML)GO 150 1 TO 50 151 40 CONTINUE 152 N1=1 153 50 IF(N2.EQ.N1) GO TO 100 154 N2M1=N2-1 155 IF(ITER.GE.3)ROOT=A22 156 K=K+1 157 A22=VALUE(N1) 158 A12=A22-ROOT 159 A23=W(N1+1) 160 A13=A23 161 DO70 I=N1,N2M1 162 A33=VALUE(I+1) 163 IF(I.NE.N2M1)A34=W (I+2) 164 S=SIGN(SQRT(A12*A12+A13*A13),A12) 165 SI=A13/S 166 CO=A12/S 167 JK=I*IV+1 168 J1=JK-IV 169 J2=J1+MIN0(M,I+K)-1 170 DO 60 JI=J1,J2 171 V1=VEC(JI) 172 V2=VEC(JK) 173 VEC(JI)=V1*CO+V2*SI 174 VEC(JK)=V2*CO-V1*SI 175 60 JK=JK+1 176 IF(I.NE.N1) W(I)=S 177 A11=CO*A22+SI*A23 178 A12=CO*A23+SI*A33 179 A13=SI*A34 180 A21=CO*A23-SI*A22 181 A22=CO*A33-SI*A23 182 A23=CO*A34 183 VALUE(I)=A11*CO+A12*SI 184 A12=-A11*SI+A12*CO 185 W(I+1)=A12 186 70 A22=A22*CO-A21*SI 187 80 VALUE(N2)=A22 188 WRITE(6,90) 189 90 FORMAT(48H1CYCLE DETECTED IN SUBROUTINE EA08 -STOPPING NOW) 190 STOP 191 100 CONTINUE 192C RAYLEIGH QUOTIENT 193 DO 120 J=1,M 194 K=(J-1)*IV 195 XX=VEC(K+1)**2 196 XAX=XX*A(1) 197 DO 110 I=2,M 198 KI=K+I 199 XX=XX+VEC(KI)**2 200 110 XAX=XAX+VEC(KI)*(2.*B(I)*VEC(KI-1)+A(I)*VEC(KI)) 201 120 VALUE(J)=XAX/XX 202 RETURN 203 END 204 SUBROUTINE EA09C(A,B,VALUE,M,OFF) 205C STANDARD FORTRAN 66 (A VERFIED PFORT SUBROUTINE) 206 DIMENSION A(*),B(*),VALUE(*),OFF(*) 207 DATA A34/0.0/,EPS/1.0E-6/ 208 SML=EPS*FLOAT(M) 209 VALUE(1)=A(1) 210 IF(M.EQ.1)RETURN 211 DO 10 I=2,M 212 VALUE(I)=A(I) 213 10 OFF(I)=B(I) 214C EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2 215 MAXIT=10*M 216 DO 80 ITER=1,MAXIT 217 DO 40 N3=2,M 218 N2=M+2-N3 219 DO 20 II=2,N2 220 N1=2+N2-II 221 IF(ABS(OFF(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*SML) 222 1GO TO 30 223 20 CONTINUE 224 N1=1 225 30 IF(N2.NE.N1) GO TO 50 226 40 CONTINUE 227 RETURN 228C ROOT IS THE EIGENVALUE OF THE BOTTOM 2*2 MATRIX THAT IS NEAREST 229C TO THE LAST MATRIX ELEMENT AND IS USED TO ACCELERATE THE 230C CONVERGENCE 231 50 BB=(VALUE(N2)-VALUE(N2-1))*0.5 232 CC=OFF(N2)*OFF(N2) 233 SBB=1. 234 IF(BB.LT.0.)SBB=-1. 235 ROOT=VALUE(N2)+CC/(BB+SBB*SQRT(BB*BB+CC)) 236 N2M1=N2-1 237 60 A22=VALUE(N1) 238 A12=A22-ROOT 239 A23=OFF(N1+1) 240 A13=A23 241 DO 70 I=N1,N2M1 242 A33=VALUE(I+1) 243 IF(I.NE.N2M1)A34=OFF(I+2) 244 S=SQRT(A12*A12+A13*A13) 245 SI=A13/S 246 CO=A12/S 247 IF(I.NE.N1)OFF(I)=S 248 A11=CO*A22+SI*A23 249 A12=CO*A23+SI*A33 250 A13=SI*A34 251 A21=CO*A23-SI*A22 252 A22=CO*A33-SI*A23 253 A23=CO*A34 254 VALUE(I)=A11*CO+A12*SI 255 A12=-A11*SI+A12*CO 256 OFF(I+1)=A12 257 70 A22=A22*CO-A21*SI 258 80 VALUE(N2)=A22 259 WRITE(6,90) 260 90 FORMAT(39H1LOOPING DETECTED IN EA09-STOPPING NOW ) 261 STOP 262 END 263 COMPLEX FUNCTION FM06AS(N,A,IA,B,IB) 264 IMPLICIT COMPLEX (A-H,O-Z) 265 COMPLEX A(*), B(*) 266******************************************************* 267* 268* FM06AS - A FUNCTION ROUTINE TO COMPUTE THE VALUE OF THE 269* INNER PRODUCT, OR DOT PRODUCT, OF TWO SINGLE PRECISION 270* COMPLEX VECTORS, ACCUMULATING THE INTERMEDIATE PRODUCTS 271* DOUBLE PRECISION. THE ELEMENTS OF EACH VECTOR CAN BE 272* STORED IN ANY FIXED DISPLACEMENT FROM NEIGHBOURING 273* ELEMENTS. 274* 275* COMPUTES: SUM(J=1,N) A((J-1)*IA+1)*B((J-1)*IB+1) 276* 277* W = FM06AS(N,A,IA,B,IB) 278* 279* N INTEGER SCALAR; (USER:*); LENGTH OF THE VECTORS A AND B. 280* IF N <= 0 THE INNER PRODUCT VALUE IS DEFINED TO BE ZERO. 281* A COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY 282* CONTAINING THE 1ST VECTOR. THE FORTRAN CONVENTION OF STORING 283* REAL AND IMAGINARY PARTS IN ADJACENT WORDS IS ASSUMED. 284* IA INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF 285* AN ELEMENT IN THE ARRAY A TO ITS NEIGHBOUR, I.E. THE VECTOR 286* ELEMENTS ARE IN A(1), A(IA+1), A(2*IA+1),... 287* IF IA < 0 THE ELEMENTS ARE ASSUMED TO BE STORED IN 288* A(1-(N-1)*IA), A(1-(N-2)*IA),..., A(1-IA), A(1). 289* B COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY 290* CONTAINING THE SECOND VECTOR. TREAT LIKE A. 291* IB INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF 292* AN ELEMENT IN B TO ITS NEIGHBOUR. TREAT LIKE IA. 293* FM06AS COMPLEX FUNCTION; (*:FM06AS); THE INNER PRODUCT VALUE. 294* IT IS RETURNED DOUBLE PRECISION, THE REAL PART IN FLT PNT 295* REG 0 AND THE IMAGINARY PART IN FLT PNT REG 2. 296* 297* THIS ROUTINE IS WRITTEN IN FORTRAN. 298* 299*--------------------------------------------------------* 300 SUM=(0.0,0.0) 301 DO 10 I=1,N 302 10 SUM=SUM+A((I-1)*IA+1)*B((I-1)*IB+1) 303 FM06AS=SUM 304 RETURN 305 END 306 COMPLEX FUNCTION FM06BS(N,A,IA,B,IB) 307 IMPLICIT COMPLEX (A-H,O-Z) 308 COMPLEX A(*), B(*) 309******************************************************* 310* 311* FM06BS - A FUNCTION ROUTINE TO COMPUTE THE VALUE OF THE 312* INNER PRODUCT, OR DOT PRODUCT, OF A SIGLE PRECISION 313* COMPLEX VECTORS, ACCUMULATING THE INTERMEDIATE PRODUCTS 314* DOUBLE PRECISION. THE ELEMENTS OF EACH VECTOR CAN BE 315* STORED IN ANY FIXED DISPLACEMENT FROM NEIGHBOURING 316* ELEMENTS. 317* 318* COMPUTES: SUM(J=1,N) A((J-1)*IA+1)*B((J-1)*IB+1) 319* 320* W = FM06BS(N,A,IA,B,IB) 321* 322* N INTEGER SCALAR; (USER:*); LENGTH OF THE VECTORS A AND B. 323* IF N <= 0 THE INNER PRODUCT VALUE IS DEFINED TO BE ZERO. 324* A COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY 325* CONTAINING THE 1ST VECTOR. THE FORTRAN CONVENTION OF STORING 326* REAL AND IMAGINARY PARTS IN ADJACENT WORDS IS ASSUMED. 327* IA INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF 328* AN ELEMENT IN THE ARRAY A TO ITS NEIGHBOUR, I.E. THE VECTOR 329* ELEMENTS ARE IN A(1), A(IA+1), A(2*IA+1),... 330* IF IA < 0 THE ELEMENTS ARE ASSUMED TO BE STORED IN 331* A(1-(N-1)*IA), A(1-(N-2)*IA),..., A(1-IA), A(1). 332* B COMPLEX*8 ARRAY((N-1)*IABS(IA)+1); (USER:*); THE ARRAY 333* CONTAINING THE SECOND VECTOR. TREAT LIKE A. 334* IB INTEGER SCALAR; (USER:*); THE SUBSCRIPT DISPLACEMENT OF 335* AN ELEMENT IN B TO ITS NEIGHBOUR. TREAT LIKE IA. 336* FM06BS COMPLEX FUNCTION; (*:FM06BS); THE INNER PRODUCT VALUE. 337* IT IS RETURNED DOUBLE PRECISION, THE REAL PART IN FLT PNT 338* REG 0 AND THE IMAGINARY PART IN FLT PNT REG 2. 339* 340* THIS ROUTINE IS WRITTEN IN FORTRAN. 341* 342*--------------------------------------------------------* 343 SUM=(0.0,0.0) 344 DO 10 I=1,N 345 10 SUM=SUM+A((I-1)*IA+1)*CONJG(B((I-1)*IB+1)) 346 FM06BS=SUM 347 RETURN 348 END 349 SUBROUTINE ME08A(A,ALPHA,BETA,N,IA,Q) 350 COMPLEX A(IA,*),ALPHA(*),BETA(*),Q(*),CW,QJ 351 COMPLEX FM06AS,FM06BS 352 REAL PP,ZERO,P5,H 353 REAL S1,PP1 354 DATA ZERO/0.0/, P5 /0.50/ 355C************************************************************** 356C THE REDUCTION OF FULL HERMITIAN MATRIX INTO TRI-DIAGONAL HERMITIAN 357C FORM IS DONE IN N-2 STEPS.AT THE I TH STEP ZEROS ARE INTRODUCED IN 358C THE I TH ROW AND COLUMNS WITHOUT DESTROYING PREVIOUSLY PRODUCED ZEROS 359C 360C H 361C AT THE I TH STEP WE HAVE A =P A P WITH P =I-U U /K 362C I I I-1 I I I I I 363C 364C WHERE U =(0,0,...,A (1+S /T ),A ,...,A ) 365C I I,I+1 I I I,I+2 I,N 366C 2 N 2 2 2 367C S = SUM ] A ] K =S +T AND T = SQRT(]A ] S ) 368C I J=I+1 I,J I I I I I,I+1 I 369C 370C COMPUTATIONAL DETAILS AT THE I TH STAGE ARE (1) FORM S ,K THEN 371C I I 372C 373C H H H 374C (2) Q =A U /K (3) Q =Q -.5U (U Q /K ) (4) A =A -U Q -Q U 375C I I-1 I I I I I I I I I I-1 I I I I 376C 377C THE VECTORS U BEING APPROPIATELY IN A. 378 IF(N.LE.0)GO TO 90 379 DO 10 J=1,N 380 ALPHA(J)=A(J,J) 381 DO 10 I=1,J 382 10 A(I,J)=CONJG(A(J,I)) 383 IF(N-2)90,80,20 384 20 N2=N-2 385 DO 60 I=1,N2 386 I1=I+1 387C (1) 388 CW=FM06BS(N-I,A(I,I+1),IA,A(I,I+1),IA) 389 PP=CW 390 PP1=SQRT(PP) 391 BETA(I+1)=CMPLX(-PP1,ZERO) 392 S1=CABS(A(I,I+1)) 393 IF(S1.GT.ZERO)BETA(I+1)=BETA(I+1)*A(I,I+1)/S1 394 IF(PP.LE.1.D-15)GO TO 60 395 H=PP+PP1*S1 396 A(I,I+1)=A(I,I+1)-BETA(I+1) 397C (2) 398 DO 30 K=I1,N 399 QJ=FM06AS(-(I-K),A(I+1,K),1,A(I,I+1),IA) 400 QJ=FM06BS(N-K,A(K,K+1),IA,A(I,K+1),IA)+CONJG(QJ) 401 30 Q(K)=QJ/H 402C (3) 403 CW=FM06AS(N-I,A(I,I+1),IA,Q(I+1),1) 404 PP=CW*P5/H 405 DO 40 K=I1,N 406 40 Q(K)=Q(K)-PP*CONJG(A(I,K)) 407C (4) 408 DO 50 K=I1,N 409 50 CALL ME08B (A(K,K),Q(K),A(I,K),N-K+1,IA*2) 410 60 CONTINUE 411 DO 70 I=2,N 412 QJ=ALPHA(I) 413 ALPHA(I)=A(I,I) 414 70 A(I,I)=QJ 415 80 BETA(N)=A(N-1,N) 416 90 RETURN 417 END 418 SUBROUTINE SORT(VAL,VEC,N) 419 COMPLEX VEC(N,*), SUM 420 REAL VAL(*) 421 DO 30 I=1,N 422 X=1.E9 423 DO 10 J=I,N 424 IF(VAL(J).LT.X) THEN 425 K=J 426 X=VAL(J) 427 ENDIF 428 10 CONTINUE 429 DO 20 J=1,N 430 SUM=VEC(J,K) 431 VEC(J,K)=VEC(J,I) 432 20 VEC(J,I)=SUM 433 VAL(K)=VAL(I) 434 VAL(I)=X 435 30 CONTINUE 436 RETURN 437 END 438