C C geigen.f ver 1.1 November 18, 2011 C C By Tom Irvine C Email: tomirvine@aol.com C C Generalized Eigenvalue Problem C C For a structural dynamics system, C C Matrix A = stiffness C Matrix B = mass C C C23456789 integer :: I,J,N,NL integer :: IDF,IW integer :: ichoice integer :: DeAllocateStatus double precision, allocatable :: A(:,:),B(:,:),W(:),fn(:) character :: A_filename*20,B_filename*20,File1*20,File2*20 double precision PI parameter (PI = 3.1415926535897932D0) intrinsic INT, MIN C * .. External Subroutines .. external DSYEV C print '(" Select analysis type")' print '(" 1=standard eigen problem")' print '(" 2=generalized eigen problem")' read(*,*) ichoice C if(ichoice.eq.1)then print '(" The input matrix must be square, real, & symmetric")' else print '(" The input matrices must be square, real, & symmetric")' endif C print '("")' print '( "Enter the number of rows for matrix A")' read(*,*)N C print '("")' print '( "Enter the filename name of matrix A (stiffness)")' read(*,*)A_filename C allocate(A(N,N)) allocate(W(N)) allocate(fn(N)) C open(unit=10,FILE=A_filename,STATUS='OLD',ACTION='read') read(10,*) ((A(I,J),I=1,N),J=1,N) C if(ichoice.eq.2)then print '( "Enter the filename of matrix B (mass)")' read(*,*)B_filename allocate(B(N,N)) open(unit=11,FILE=B_filename,STATUS='OLD',ACTION='read') read(11,*) ((B(I,J),I=1,N),J=1,N) endif C if(ichoice.eq.1)then C call STANDARD_EIGEN(A,N,W) C else C call GENERALIZED_EIGEN(A,B,N,W) C endif C C Print eigenvalues. C print '("")' print '(" i eigen fn(Hz) ")' NL=MIN(N,20) do i=1,NL if(W(i).GT.0.)then fn(i)=sqrt(W(i))/(2*PI) else fn(i)=0.; endif write(*,101)i,W(i),fn(i) 101 format(1X,I4,2X,E10.4,2X,E10.4) end do C C Print eigenvectors. C print '("")' print '( "Display Eigenvectors? 1=yes 2=no ")' read(*,*)IDR C if(IDR.eq.1)then call PRINT_MATRIX( 'Eigenvectors (columnwise)', N, N, A, N ) endif C print '("")' print '( "Write output data to files? 1=yes 2=no ")' read(*,*)IW C if(IW.eq.1)then C print '("")' print '( "Enter the eigenvalue filename")' read(*,*)File1 C print '("")' print '( "Enter the eigenvector")' read(*,*)File2 C open(unit=20,FILE=File1) open(unit=21,FILE=File2) C C write(20,333) 333 format(" i eigen fn(Hz) ") do i=1,N if(W(i).GT.0.)then fn(i)=sqrt(W(i))/(2*PI) else fn(i)=0.; endif write(20,101)i,W(i),fn(i) end do C do i = 1, N write(21,9998) ( A( i, j ), j = 1, N ) end do C 9998 format( 11(:,1X,E10.4) ) C endif C deallocate (A, STAT = DeAllocateStatus) deallocate (B, STAT = DeAllocateStatus) deallocate (W, STAT = DeAllocateStatus) deallocate (fn, STAT = DeAllocateStatus) C stop end C C***************************************************************************************************** C subroutine STANDARD_EIGEN(A,N,W) integer N integer IDR integer LDA integer LWMAX parameter ( LWMAX = 100000 ) double precision A(N,N),W(N),WORK( LWMAX ) integer INFO, LWORK intrinsic INT, MIN real :: START, FINISH C LDA=N C do i=1,N do j=i+1,N A(j,i)=0.; end do end do C C Query the optimal workspace. C call CPU_TIME(START) LWORK = -1 call DSYEV( 'Vectors', 'Upper', N, A, LDA, W, WORK, LWORK, INFO) LWORK = MIN( LWMAX, INT( WORK( 1 ) ) ) C C Solve eigenproblem. C call DSYEV( 'Vectors', 'Upper', N, A, LDA, W, WORK, LWORK, INFO) call CPU_TIME(FINISH) C C Check for convergence. C if( INFO.GT.0 ) then write(*,*)'The algorithm failed to compute eigenvalues.' stop endif C print '("")' print '("Time = ",F8.3," seconds")',FINISH-START C return end C C***************************************************************************************************** C subroutine GENERALIZED_EIGEN(A,B,N,W) integer N integer IDR integer LDA,LDB integer LWMAX parameter ( LWMAX = 100000 ) double precision A(N,N),B(N,N),W(N),WORK( LWMAX ) integer ITYPE, INFO, LWORK intrinsic INT, MIN real :: START, FINISH C ITYPE=1 LDA=N LDB=N C do i=1,N do j=i+1,N A(j,i)=0. B(j,i)=0. end do end do C C Query the optimal workspace. C call CPU_TIME(START) LWORK = -1 call DSYGV( ITYPE, 'V', 'UP', N, A, LDA, B, LDB, W, WORK, $ LWORK, INFO ) LWORK = MIN( LWMAX, INT( WORK( 1 ) ) ) C C Solve eigenproblem. C call DSYGV( ITYPE, 'V', 'UP', N, A, LDA, B, LDB, W, WORK, $ LWORK, INFO ) call CPU_TIME(FINISH) C C Check for convergence. C if( INFO.GT.0 ) then write(*,*)'The algorithm failed to compute eigenvalues.' stop endif C print '("")' print '("Time = ",F8.3," seconds")',FINISH-START C return end C C***************************************************************************************************** C C Auxiliary routine: printing a matrix. C subroutine PRINT_MATRIX( DESC, M, N, A, LDA ) character*(*) DESC integer M, N, LDA double precision A( LDA, * ) C integer I, J C write(*,*) write(*,*) DESC do I = 1, M write(*,9998) ( A( I, J ), J = 1, N ) end do * 9998 format( 11(:,1X,F7.3) ) return end