UCL WIKI

UCL Logo
Child pages
  • diag_syev (LAPACK)
Skip to end of metadata
Go to start of metadata

This programs is for diagonalizing a real symmetric matrix using DSYEV from the LAPACK library (see http://www.netlib.org/lapack/explore-3.1.1-html/dsyev.f.html). The matrix is given in the file mat.dat. 

Create a Code::Blocks project with the following code and compile. You will need to activate the MKL option:

Project->Build options->Linker settings->Other linker options

Write -mkl in the right panel:

 

In order to compile from the command line the MKL libraries must be used as follows:

ifort -mkl diatest.90 -o diatest.x

or for the OpenMP parallel  version

ifort -mkl=parallel diatest.90 -o diatest.x

 

Exercise 1: Prepare mat.dat with a 2x2 symmetric matrix mat.dat in a suitable format (see the program) and run with diatest.  

Exercise 2: Use the 28x28 and 56x56 matrices from the data files mat28.dat  and mat56.dat.

Exercise 3: Use the LAPACK subroutine DPOTRI to invert the matrix: http://www.netlib.org/lapack/explore-3.1.1-html/dpotri.f.html

 

!------------------------------------------------------------------------!
! Example of matrix diagonalization using LAPACK routine dsyev.f and the ! 
! Fortran 90 interface diasym.f90.                                       !
!------------------------------------------------------------------------!
! Input from file 'mat.dat' (matrix to be diagonalized):                 !
! line 1      : order of the symmetric matrix [M]                        !
! lines 2-n+1 : rows of the matrix                                       !
!------------------------------------------------------------------------!
! Output in file 'dia.dat':                                              !
! - eigenvalues                                                          !
! - eigenvectors (diagonalizing matrix [D])                              !
! - the original matrix [M] transformed by [D]; [1/D][M][D]              !
!------------------------------------------------------------------------!
!---------------!
 program diatest
!---------------!
 implicit none
 integer :: i,j,n,info
 real(8), allocatable :: m0(:,:),m1(:,:),m2(:,:),eig(:)
 real(8) :: elem
 open(10,file='mat.dat',status='old')
 read(10,*)n
 allocate (m0(n,n))
 allocate (m1(n,n))
 allocate (m2(n,n))
 allocate (eig(n))
 !
 m0 = 0
 !
 do while(.true.)
   read( 10, *, iostat=info ) i, j, elem
   if ( info /= 0 ) exit   ! exit loop when reach EOF
   m0(i,j) = elem
 enddo
 close(10)
 m1(:,:)=m0(:,:)
 call diasym(m1,eig,n)
   
 open(10,file='dia.dat',status='replace')
 write(10,*)'Eigenvalues:'
 do i=1,n
    write(10,10)i,eig(i)
    10 format(I3,'   ',f14.8)
 enddo
 write(10,*)
 write(10,*)'Eigenvectors:'
 do i=1,n
    write(10,20)i,m1(:,i)
    20 format(i3,'   ',100f14.8)
 enddo
 write(10,*)
 
 m2=matmul(transpose(m1),m0)
 m0=matmul(m2,m1)
 write(10,*)'Transformed matrix (check):'
 do i=1,n
    write(10,30)m0(:,i)
    30 format(100f14.8)
 enddo
 write(10,*)
 close(10)
 deallocate(m0); deallocate(m1); deallocate(m2);  deallocate(eig)
 end program diatest
!-------------------!
!---------------------------------------------------------!
!Calls the LAPACK diagonalization subroutine DSYEV        !
!input:  a(n,n) = real symmetric matrix to be diagonalized!
!            n  = size of a                               !
!output: a(n,n) = orthonormal eigenvectors of a           !
!        eig(n) = eigenvalues of a in ascending order     !
!---------------------------------------------------------!
 subroutine diasym(a,eig,n)
 implicit none
 integer n,l,inf
 real*8  a(n,n),eig(n),work(n*(3+n/2))
 l=n*(3+n/2)
 call dsyev('V','U',n,a,n,eig,work,l,inf)
 end subroutine diasym
!---------------------!

 

 

 

 

 

  • No labels