Hello,
I’m wondering if someone can help me with the following problem. I
have a complex matrix (see the input data “input2.txt” below). In the
matrix data set,
the first element is the real part and the second element is the
imaginary part in each cell. I would like to read this matrix and
calculate the igen values using ZHPEV subroutine of LAPACK. The
problem I’m facing is that the results that I get by including the
data directly in the program (program 1) are different from the
results that I get by reading the input data in my program(program
2)from a separate data file (input2.txt), which is very puzzling to
me. Could someone please look into my data, programs, and output and
help me with this? Another puzzling thing is that I’m using the input
data from an example and the neither of my results match with the
actual results of the example. Following is the link for the ZHPEV
subroutine of LAPACK and I’m trying to replicate the example 6 of this
document.
http://publib.boulder.ibm.com/doc_link/en_US/a_doc_lib/sp34/essl/essl192.html
I’m using Fortran 77 running on AIX. Following are my programs and
data file:
--------------------begin program1----------
-----Here Data are included in the program--
c finding the eigenvalues of a complex matrix using LAPACK
Implicit none
c declarations, notice double precision
complex*16 A(4,4), ap(10), aux(12,12)
real*8 w(4)
integer i, j, n, iopt, naux, k
open(2,file='test1004.out')
c define matrix A
A(1,1)=(3.0, 0.0)
A(1,2)=(1.0, 0.0)
A(1,3)=(0.0, 0.0)
A(1,4)=(0.0, 2.0)
A(2,1)=(1.0, 0.0)
A(2,2)=(3.0, 0.0)
A(2,3)=(0.0,-2.0)
A(2,4)=(0.0, 0.0)
A(3,1)=(0.0, 0.0)
A(3,2)=(0.0, 2.0)
A(3,3)=(1.0, 0.0)
A(3,4)=(1.0, 0.0)
A(4,1)=(0.0,-2.0)
A(4,2)=(0.0, 0.0)
A(4,3)=(1.0, 0.0)
A(4,4)=(1.0, 0.0)
do 10 i=1,4
write(2,*)(A(i,j),j=1,4)
10 continue
c find the solution using the LAPACK routine ZHPEV
iopt=20
n=4
naux=12
k = 0
do 1 j=1,n
do 2 i=1,j
k = k+1
ap(k)=A(i,j)
2 continue
1 continue
do 50 k=1,10
write(2,*)(ap(k))
50 continue
call ZHPEV(iopt,ap,w,n,aux,naux)
do 20 i=1,4
write(2,*)w(i)
20 continue
stop
end
----------------end program 1---------------
----------------begin program2----------
---Here program reading the data from input2.txt--
c finding the eigenvalues of a complex hermitian matrix using LAPACK
Implicit none
complex*16 A(8,8), ap(10), aux(12,12)
real*8 w(4)
integer i, j, n, iopt, ldz, naux, k ,nmax
open(1,file='input2.txt')
open(2,file='test1005.out')
n=4
nmax=8
if (n.le.nmax) then
read(1,*)((A(i,j),j=1,n),i=1,n)
do 10 i=1,n
write(2,*)(A(i,j),j=1,n)
10 continue
c find the solution using the LAPACK routine ZHPEV
iopt=20
ldz=1
naux=12
k = 0
do 1 j=1,n
do 2 i=1,j
k = k+1
ap(k)=A(i,j)
2 continue
1 continue
do 60 k=1,10
write(2,*)(ap(k))
60 continue
call ZHPEV(iopt,ap,w,n,aux,naux)
do 50 i=1,4
write(2,*)w(i)
50 continue
endif
stop
end
----------------end program 2---------------
------------begin input2.txt--------------
(3.0,0.0) (1.0,0.0) (0.0,0.0) (0.0,2.0)
(1.0,0.0) (3.0,0.0) (0.0,-2.0) (0.0,0.0)
(0.0,0.0) (0.0,2.0) (1.0,0.0) (1.0,0.0)
(0.0,-2.0) (0.0,0.0) (1.0,0.0) (1.0,0.0)
------------end input2.txt-----------------
--begin output from program1(test1004.out)------
(3.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,2.00000000000000000)
(1.00000000000000000,0.000000000000000000E+00)
(3.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,-2.00000000000000000)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,2.00000000000000000)
(1.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,-2.00000000000000000)
(0.000000000000000000E+00,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(3.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(3.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,-2.00000000000000000)
(1.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,2.00000000000000000)
(0.000000000000000000E+00,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
-4.89540187687792017
-1.42262165624283488
-0.842671673797785403
-0.979369184259774178E-01
--end output from program1--------------------
--begin output from program2(test1005.out)----
(3.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,2.00000000000000000)
(1.00000000000000000,0.000000000000000000E+00)
(3.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,-2.00000000000000000)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,2.00000000000000000)
(1.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,-2.00000000000000000)
(0.000000000000000000E+00,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(3.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(3.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,0.000000000000000000E+00)
(0.000000000000000000E+00,-2.00000000000000000)
(1.00000000000000000,0.000000000000000000E+00)
(0.000000000000000000E+00,2.00000000000000000)
(0.000000000000000000E+00,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
(1.00000000000000000,0.000000000000000000E+00)
-6.51063402820271531
-3.50424563545021117
-1.85595028161575293
-0.621359588192943679
--end output from program2--------------------
Note that the last four rows in the output are the eigen values which
should be exactly same resulting from both the programs and also, with
the example from the web. But they are not.
Please let me know if you have any clarifying question. |
