C
C Subroutine to compute eigenvalues and eigenvectors of a real C symmetric matrix, using subroutines from Numerical Recipes. Form C of the subroutine is the same as the form of subroutine 'EIGEN' C from the IBM Scientific Subroutine Package. C C Version: 6.70 C Date: 2012, November 7 C C Coded by Petr Bulant C http://sw3d.cz/staff/bulant.htm C C----------------------------------------------------------------------- C Usage: C CALL EIGEN(A,R,N,MV) C C Input: C A... Original matrix (symmetric), destroyed in computation. C N... Order of matrices A and R. C MV... Input code: C 0... Compute eigenvalues and eigenvectors. C 1... Compute eigenvalues only (currently not possible). C Output: C A... Eigenvalues of the input matrix are developed in diagonal C of matrix A in descending order. C R... Resultant matrix of eigenvectors (stored columnwise, C in same sequence as eigenvalues). C C Remarks: C Original matrix A must be real symmetric (storage mode=1). C Matrix A cannot be in the same location as matrix R. C Value of MAXRAM must be set by calling program, C MRAM-MAXRAM.GE.2*N. C C----------------------------------------------------------------------- SUBROUTINE EIGEN(A,R,N,MV) INTEGER N,MV REAL A(N*(N+1)/2+1),R(N,N) C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C REAL WORK(MRAM) INTEGER I1,I2,I3,INDX(MRAM) EQUIVALENCE (WORK,RAM) EQUIVALENCE (INDX,IRAM) C Subroutines and external functions required: EXTERNAL ERROR,TRED2,TQLI,INDEXX C ERROR... File error.for. C TRED2... File tred2.for of Numerical Recipes. C TQLI... File tqli.for of Numerical Recipes. C INDEXX... File indexx.for of Numerical Recipes. C----------------------------------------------------------------------- IF (N.GT.(MRAM-MAXRAM)/2) THEN C EIGENNR-03 CALL ERROR('EIGENNR-03: Small dimension of auxiliary arrays.') C Value of MAXRAM must be set by the calling program, C MRAM-MAXRAM must be greater or equal 2*N. ENDIF C IF (MV.EQ.0) THEN C Computation of eigenvalues and eigenvectors. C C Preparing matrix R: I3=0 DO 12, I1=1,N DO 10, I2=1,I1 I3=I3+1 R(I1,I2)=A(I3) R(I2,I1)=A(I3) 10 CONTINUE 12 CONTINUE C C Preparing a tridiagonal matrix from R: CALL TRED2(R,N,N,A,A(N+1)) C C Computing eigenvalues and eigenvectors of matrix R: CALL TQLI(A,A(N+1),N,N,R) C C Sorting the eigenvalues into descending order: CALL INDEXX(N,A,INDX(MAXRAM+N+1)) C Reorganizing eigenvalues: DO 14, I2=1,N C Store the element to WORK: WORK(MAXRAM+I2)=A(I2) 14 CONTINUE DO 16, I2=1,N C Copy it back in the rearranged order: A(I2)=WORK(MAXRAM+INDX(MAXRAM+N+N+1-I2)) 16 CONTINUE C Reorganizing eigenvectors: C Loop over lines: DO 25, I1=1,N C For each element of the line: DO 18, I2=1,N C Store the element to WORK: WORK(MAXRAM+I2)=R(I1,I2) 18 CONTINUE DO 20, I2=1,N C Copy it back in the rearranged order: R(I1,I2)=WORK(MAXRAM+INDX(MAXRAM+N+N+1-I2)) 20 CONTINUE 25 CONTINUE C C Writing eigenvalues to the diagonal of A: DO 27, I1=N,2,-1 A((I1*(I1+1))/2)=A(I1) 27 CONTINUE ELSEIF (MV.EQ.1) THEN C Computation of eigenvalues only. C EIGENNR-01 CALL ERROR('EIGENNR-01: MV=1 currently not supported.') ELSE C EIGENNR-02 CALL ERROR('EIGENNR-02: Wrong value of MV.') ENDIF RETURN END