Hi,

Can any one recommend a reliable, fast, and easy to find fortran
routine to diagonalize a general real matrix?  I need routines for two
applications.  In the first I need only find one basis vector with an
eigenvalue close to one.  In the other I need only the lowest couple
(or many) eigenvectors. 

Thanks, Rudy

rjmagyar wrote:
> Hi,
>
> Can any one recommend a reliable, fast, and easy to find fortran
> routine to diagonalize a general real matrix?  I need routines for two
> applications.  In the first I need only find one basis vector with an
> eigenvalue close to one.  In the other I need only the lowest couple
> (or many) eigenvectors.

To get all eigenvectors, Lapack is reliable and fast, but not so easy
to install and use (but certainly doable). Eispack is not as fast as
Lapack but can be a little easier to work with. You can get a Fortran
90 version of Eispack from
http://www.csit.fsu.edu/~burkardt/f_src/eispack/eispack.html or get the
original F77 from Netlib, which also has Lapack.

To get the eigenvector with eigenvalue near one (or any specified
number) you can use the method of "inverse iteration". I am not sure if
there is a Fortran library that implements this method, but it is not
difficult to program.

LAPACK is the way to go

http://www.netlib.org/lapack/double/

Michael

rjmagyar wrote:
> Hi,
>
> Can any one recommend a reliable, fast, and easy to find fortran
> routine to diagonalize a general real matrix?  I need routines for two
> applications.  In the first I need only find one basis vector with an
> eigenvalue close to one.  In the other I need only the lowest couple
> (or many) eigenvectors. 
> 
> Thanks, Rudy

rjmagyar wrote:
> Hi,
> 
> Can any one recommend a reliable, fast, and easy to find fortran
> routine to diagonalize a general real matrix?  I need routines for two
> applications.  In the first I need only find one basis vector with an
> eigenvalue close to one.  In the other I need only the lowest couple
> (or many) eigenvectors. 
> 
> Thanks, Rudy
> 

What OS are you using?

raju

> To get the eigenvector with eigenvalue near one (or any specified
> number) you can use the method of "inverse iteration". I am not sure if
> there is a Fortran library that implements this method, but it is not
> difficult to program.

xHSEIN from LAPACK.

I would like something portable between many different OS's, but I am
mostly running on altix machines.

Kamaraju Kusumanchi wrote:
> rjmagyar wrote:
> > Hi,
> >
> > Can any one recommend a reliable, fast, and easy to find fortran
> > routine to diagonalize a general real matrix?  I need routines for two
> > applications.  In the first I need only find one basis vector with an
> > eigenvalue close to one.  In the other I need only the lowest couple
> > (or many) eigenvectors.
> > 
> > Thanks, Rudy
> > 
> 
> What OS are you using?
> 
> raju

Any suggestions as to the best way to take the general matrix to uppser
Hessenberg form?

Thanks, Rudy

highegg wrote:
> > To get the eigenvector with eigenvalue near one (or any specified
> > number) you can use the method of "inverse iteration". I am not sure if
> > there is a Fortran library that implements this method, but it is not
> > difficult to program.
> 
> xHSEIN from LAPACK.

rjmagyar wrote:
> 
> I would like something portable between many different OS's, but I am
> mostly running on altix machines.

LAPACK should be included, as part of ProPack, on your Altix machine.
Just link with -lscs (or -lscs_mp for the parallel version.)

An alternative version of LAPACK is contained in the Intel MKL library.
If your Altix machine has MKL, it will typically reside in /opt/intel/mkl/...

Walt

I was wondering if any of you have used to DGEEV routine from netlib.
I can't get the routine to properly diagnolize matrices of inputted
length.  The only way I can get the routine to work is when the array
is completely full.  For example,

real*8 M(max)

For I=1,N
For J=1,N
M(I,J)=i*j
enddo
enddo

Call DGEEV(...M,N...)

Only works when N=max
 
Any ideas how to work with N<MAX?

Thanks, Rudy

rjmagyar wrote:
> I was wondering if any of you have used to DGEEV routine from netlib.
> I can't get the routine to properly diagnolize matrices of inputted
> length.  The only way I can get the routine to work is when the array
> is completely full.  For example,

> real*8 M(max)

> For I=1,N
> For J=1,N
> M(I,J)=i*j
> enddo
> enddo

> Call DGEEV(...M,N...)

> Only works when N=max

I presume you mean N**2=max...

LDA tells GDEEV what the actual value is for the first dimension
of M if it isn't N.  This is needed so DGEEV can find the array
elements in memory.  If the array is larger than N*N set LDA
to the actual dimension if the array.

call DGEEV(...M,N,LDA,...)

*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*


http://www.netlib.org/lapack/double/dgeev.f

-- glen

I still get the same trouble.  Here's an example:

This works . . .
C     ****f* diag_test/main[0.1]
C
C     NAME
C     Diag_Test
C
***

C--------------------------------------
C     Declare and Define
C--------------------------------------
      implicit none

      integer maxbas
      parameter(maxbas=2)
      real*8 s(1:maxbas,1:maxbas)
      integer i,j,nbasis

C     just for the lapack routine
      real*8 WR(maxbas), WI(maxbas)
      real*8 VL(1,maxbas),vr(maxbas,maxbas)
      integer ldvl
C     parameter (ldvl=1)
      integer lwork
      parameter (lwork=8*maxbas)
      real*8 work(lwork)
      integer info

C--------------------------------------
C     Set up Trial Matrix
C--------------------------------------
      NBASIS=2
      S(1,1)=1.
      S(2,1)=2.
      S(1,2)=2.
      S(2,2)=4.
C     evalues should be 0,5
C     eigenvectors should be (-2,1)  and (1,2)

      do i=1,nbasis
         write(6,*) (s(i,j), j=1,nbasis)
      enddo

      write(6,*) s

C--------------------------------------
C     Diagnolatize Matrix
C--------------------------------------
C diagnonalize the matrix S
      call DGEEV( 'N','V', nbasis,s,nbasis,
     &     WR, WI, VL, nbasis, vr,
     &     nbasis, WORK, LWORK, INFO )

      write(6,*) 'done.'
      write(6,*) 'eigenvalues:'
      write(6,*) wr


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc

C     That's it.  We have reached the end.
      end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%%%%%%%%%%%END PROGRAM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

With a write statement added into DGEEV. . .

       write(6,*) 'in dgeev ',lda,n,ldvl,ldvr

       do i=1,lda
       write(6,*)  (A(i,j),j=1,n)
       enddo



The output is. . .

  1.  2.
  2.  4.
  1.  2.  2.  4.
 in dgeev  2 2 2 2
  1.  2.
  2.  4.
 done.
 eigenvalues:
  0.  5.



Which is correct.


But when I do . . .

C     ****f* diag_test/main[0.1]
C
C     NAME
C     Diag_Test
C
***

C--------------------------------------
C     Declare and Define
C--------------------------------------
      implicit none

      integer maxbas
      parameter(maxbas=3)
      real*8 s(1:maxbas,1:maxbas)
      integer i,j,nbasis

C     just for the lapack routine
      real*8 WR(maxbas), WI(maxbas)
      real*8 VL(1,maxbas),vr(maxbas,maxbas)
      integer ldvl
C     parameter (ldvl=1)
      integer lwork
      parameter (lwork=8*maxbas)
      real*8 work(lwork)
      integer info

C--------------------------------------
C     Set up Trial Matrix
C--------------------------------------
      NBASIS=2
      S(1,1)=1.
      S(2,1)=2.
      S(1,2)=2.
      S(2,2)=4.
C     evalues should be 0,5
C     eigenvectors should be (-2,1)  and (1,2)

      do i=1,nbasis
         write(6,*) (s(i,j), j=1,nbasis)
      enddo

      write(6,*) s

C--------------------------------------
C     Diagnolatize Matrix
C--------------------------------------
C diagnonalize the matrix S
      call DGEEV( 'N','V', nbasis,s,nbasis,
     &     WR, WI, VL, nbasis, vr,
     &     nbasis, WORK, LWORK, INFO )

      write(6,*) 'done.'
      write(6,*) 'eigenvalues:'
      write(6,*) wr


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc

C     That's it.  We have reached the end.
      end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%%%%%%%%%%%END PROGRAM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I get . . .

  1.  2.
  2.  4.
  1.  2.  5.69654448E-270  2.  4.  2.94053602  5.69654448E-270
5.6421269E-270
  6.6563259E-316
 in dgeev  2 2 2 2
  1.  5.69654448E-270
  2.  2.
 done.
 eigenvalues:
  1.  2.  2.10288763



Any ideas what's wrong?  

Thanks, Rudy

rjmagyar wrote:

(snip)

> But when I do . . .

(snip)

>       integer maxbas
>       parameter(maxbas=3)
>       real*8 s(1:maxbas,1:maxbas)
(snip)

-----
> C diagnonalize the matrix S
>       call DGEEV( 'N','V', nbasis,s,nbasis,
>      &     WR, WI, VL, nbasis, vr,
>      &     nbasis, WORK, LWORK, INFO )

The fifth argument, LDA, needs to be maxbas.

Your main program has maxbas=3, and so s is in memory as:

s(1,1) s(2,1) s(3,1) s(1,2) s(2,2) s(3,2) s(1,3) s(2,3) s(3,3)

LDA is used by DGEEV to find the appropriate elements.  If you
set LDA to 2, DGEEV will use the memory as

s(1,1) s(2,1) s(1,2) s(2,2)

such that s(1,2) and s(2,2) are stored in the wrong place.

In addition, the ninth and eleventh argument should match the
left dimension of the VL and VR arrays, respectively.
That seems to be 1 and maxbas in your example, though I
wonder if VL should also be (maxbas,maxbas).

-- glen