SUBROUTINE MtxCnv( A, B, N, N1, W, IER ) EAB00100 C***********************************************************************EAB00200 C* *EAB00300 C* PURPOSE *EAB00400 C* SOLVES GENERALIZED EIGENVALUE PROBLEM *EAB00500 C* (A) X = LAMBDA (B) X *EAB00600 C* FOR REAL SYMMETRIC MATRICES (A) AND (B), *EAB00700 C* THE LATTER BEING POSITIVE DEFINITE. *EAB00800 C* DOUBLE PRECISION. *EAB00900 C* *EAB01000 C* USAGE *EAB01100 C* CALL MtxCnv( A, B, N, N1, W, IER ) *EAB01200 C* *EAB01300 C* DESCRIPTION OF PARAMETERS *EAB01400 C* (INPUT) *EAB01500 C* A(N1,N) - 2-DIM ARRAY CONTAINING THE REAL SYMMETRIC MATRIX. *EAB01600 C* B(N1,N) - 2-DIM ARRAY CONTAINING THE REAL SYMMETRIC *EAB01700 C* POSITIVE-DEFINITE MATRIX. *EAB01800 C* *EAB01900 C* N - ORDER OF THE MATRIX. *EAB02000 C* N1 - INTEGER WHICH SPECIFIES THE SIZE OF THE ARRAYS *EAB02100 C* A, B, AND V. *EAB02200 C* *EAB02600 C* (WORK AREA) *EAB03200 C* W(N,7) - WORK AREA USED IN COMPUTATION OF EIGENVALUES AND *EAB03300 C* EIGENVECTORS *EAB03400 C* IER - RETURN CODE *EAB03410 C* 0: NORMAL RETURN *EAB03420 C* NONZERO: ERROR RETURN *EAB03430 C* *EAB03500 C* REMARK *EAB03600 C* 1. THE ARRAYS A AND B WILL NOT BE CONSERVED. *EAB03700 C* 2. IT SUFFICES TO FILL A(I,J) AND B(I,J) FOR I .GE. J. *EAB03800 C* *EAB03900 C* SUBROUTINE REQUIRED *EAB04000 C* EIGRS *EAB04100 C* *EAB04200 C* COPYRIGHT *EAB04300 C* Y. ONODERA, AUGUST, 1989. *EAB04400 C* *EAB04500 C***********************************************************************EAB04600 IMPLICIT REAL*8 (A-H,O-Z) EAB04700 IMPLICIT INTEGER (I-N) EAB04800 real*8 A(N1,N), B(N1,N), W(N,7) EAB04900 IF ( N1 .LE. 0 ) GOTO 8 EAB05100 IF ( N .LE. 0 .OR. N .GT. N1 ) GOTO 8 EAB05200 GOTO 10 EAB05500 8 WRITE(6,9) N, N1 EAB05600 9 FORMAT(' (SUBR EIGAB) INVALID ARGUMENT. N, N1 =', EAB05700 & 4I5) EAB05800 IER = 2 EAB05850 RETURN EAB05900 10 CONTINUE EAB06000 IF ( N .NE. 1 ) GOTO 20 EAB06100 T = B(1,1) EAB06200 IF ( T .LE. 0 ) GOTO 15 EAB06300 c E(1) = A(1,1) / T EAB06400 c IF ( NV .EQ. 1 ) V(1,1) = 1 EAB06500 RETURN EAB06600 15 WRITE(6,16) EAB06700 16 FORMAT( ' (SUBR EIGAB) MATRIX B IS NOT POSITIVE DEFINITE.') EAB06800 IER = 1 EAB06850 RETURN EAB06900 C *************************************************************** EAB07000 C MODIFIED CHOLESKY FACTORIZATION STARTS: EAB07100 C (B) = (L TRANSPOSED) (D) (L), EAB07200 C WHERE (D) IS DIAGONAL MATRIX, EAB07300 C (L) IS LOWER TRIANGULAR MATRIX. EAB07400 C *************************************************************** EAB07500 20 DO 24 K=N,1,-1 EAB07600 SUM = 0 EAB07700 DO 21 I=K+1,N EAB07800 SUM = SUM + W(I,1) * B(I,K)**2 EAB07900 21 CONTINUE EAB08000 D = B(K,K) - SUM EAB08100 IF ( D .LE. 0 ) GOTO 15 EAB08200 W(K,1) = D EAB08300 D = 1 / D EAB08400 B(K,K) = SQRT(D) EAB08500 DO 23 J=1,K-1 EAB08600 SUM = 0 EAB08700 DO 22 I=K+1,N EAB08800 SUM = SUM + W(I,1) * B(I,K) * B(I,J) EAB08900 22 CONTINUE EAB09000 B(K,J) = ( B(K,J) - SUM ) * D EAB09100 23 CONTINUE EAB09200 24 CONTINUE EAB09300 C *************************************************************** EAB09400 C TRANSFORM THE MATRIX EAB09500 C (A NEW) = (L TRANSPOSED)**(-1) (A) (L)**(-1) EAB09600 C *************************************************************** EAB09700 DO 36 K=N-1,1,-1 EAB09800 DO 32 J=1,K-1 EAB09900 SUM = 0 EAB10000 DO 31 I=K+1,N EAB10100 SUM = SUM + B(I,K) * A(I,J) EAB10200 31 CONTINUE EAB10300 A(K,J) = A(K,J) - SUM EAB10400 32 CONTINUE EAB10500 SUMK = 0 EAB10600 DO 35 J=K+1,N EAB10700 SUM = 0 EAB10800 DO 33 I=K+1,J-1 EAB10900 SUM = SUM + B(I,K) * A(J,I) EAB11000 33 CONTINUE EAB11100 DO 34 I=J,N EAB11200 SUM = SUM + B(I,K) * A(I,J) EAB11300 34 CONTINUE EAB11400 T = A(J,K) - SUM EAB11500 SUMK = SUMK + B(J,K) * ( A(J,K) + T ) EAB11600 A(J,K) = T EAB11700 35 CONTINUE EAB11800 A(K,K) = A(K,K) - SUMK EAB11900 36 CONTINUE EAB12000 C *************************************************************** EAB12100 C PRE-MULTIPLY AND POST-MULTIPLY (A) WITH (D)**(-1/2) EAB12200 C *************************************************************** EAB12300 DO 42 I=1,N EAB12400 A(I,I) = A(I,I) / W(I,1) EAB12500 T=B(I,I) EAB12600 DO 41 J=1,I-1 EAB12700 A(I,J) = A(I,J) * B(J,J) * T EAB12800 a(j,i) = a(i,j) 41 CONTINUE EAB12900 42 CONTINUE EAB13000 RETURN EAB15300 END EAB15400