本征矢必须跟着本征值一块走，用的哪个本征值，就要用对应的本征矢，不知楼主程序考虑了没有，似乎没有！

2.2.5 Complex Hermitian Matrices
All eigenvalues .....................................................................EVLHF 502
All eigenvalues and eigenvectors........................................ EVCHF

这两个或许对楼主有用！

我查了一下EVCHF的说明，不知是否有用。

Computes all of the eigenvalues and eigenvectors of a complex Hermitian matrix.

Required Arguments
A — Complex Hermitian matrix of order N.   (Input) 
Only the upper triangle is used.

EVAL — Real vector of length N containing the eigenvalues of A in decreasing order of magnitude.   (Output)

EVEC — Complex matrix of order N.   (Output) 
The J-th eigenvector, corresponding to EVAL(J), is stored in the J-th column. Each vector is normalized to have Euclidean length equal to the value one.

Optional Arguments 
N — Order of the matrix A.   (Input)
Default: N = size (A,2).

LDA — Leading dimension of A exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDA = size (A,1).

LDEVEC — Leading dimension of EVEC exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDEVEC = size (EVEC,1).

FORTRAN 90 Interface
Generic:          CALL EVCHF (A, EVAL, EVEC [,…])

Specific:         The specific interface names are S_EVCHF and D_EVCHF.

FORTRAN 77 Interface
Single:            CALL EVCHF (N, A, LDA, EVAL, EVEC, LDEVEC)

Double:          The double precision name is DEVCHF.

恩。我试下

Description
Routine EVCHF computes the eigenvalues and eigenvectors of a complex Hermitian matrix. Unitary similarity transformations are used to reduce the matrix to an equivalent real symmetric tridiagonal matrix. The implicit QL algorithm is used to compute the eigenvalues and eigenvectors of this tridiagonal matrix. These eigenvectors and the transformations used to reduce the matrix to tridiagonal form are combined to obtain the eigenvectors for the user’s problem. The underlying code is based on either EISPACK or LAPACK code depending upon which supporting libraries are used during linking. For a detailed explanation, see “Using ScaLAPACK, LAPACK, LINPACK, and EISPACK” in the Introduction section of this manual.

Comments
1.         Workspace may be explicitly provided, if desired, by use of E5CHF/DE5CHF. The reference is:

CALL E5CHF (N, A, LDA, EVAL, EVEC, LDEVEC, ACOPY, RWK, CWK, IWK)

The additional arguments are as follows:

ACOPY — Complex work array of length N2. A and ACOPY may be the same, in which case A will be destroyed.

RWK — Work array of length N2 + N.

CWK — Complex work array of length 2N.

IWK — Integer work array of length N.

2.         Informational error

Type   Code

3           1                  The matrix is not Hermitian. It has a diagonal entry with a small imaginary part.

4           1                  The iteration for an eigenvalue failed to converge.

4           2                  The matrix is not Hermitian. It has a diagonal entry with an imaginary part.

3.         The success of this routine can be checked using EPIHF.

4.         Integer Options with Chapter 11 Options Manager

1          This option uses eight values to solve memory bank conflict (access inefficiency) problems. In routine E5CHF, the internal or working leading dimensions of ACOPY and ECOPY are both increased by IVAL(3) when N is a multiple of IVAL(4). The values IVAL(3) and IVAL(4) are temporarily replaced by IVAL(1) and IVAL(2), respectively, in routine EVCHF. Additional memory allocation and option value restoration are automatically done in EVCHF. There is no requirement that users change existing applications that use EVCHF or E5CHF. Default values for the option are IVAL(*) = 1, 16, 0, 1, 1, 16, 0, 1. Items 5-8 in IVAL(*) are for the generalized eigenvalue problem and are not used in EVCHF.

Example
In this example, a DATA statement is used to set A to a complex Hermitian matrix. The eigenvalues and eigenvectors of this matrix are computed and printed. The performance index is also computed and printed. This serves as a check on the computations, for more details, see routine EPIHF.

 

      USE IMSL_libraries

 

      IMPLICIT   NONE

 

!                                 Declare variables

      INTEGER    LDA, LDEVEC, N

      PARAMETER  (N=3, LDA=N, LDEVEC=N)

!

      INTEGER    NOUT

      REAL       EVAL(N), PI

      COMPLEX    A(LDA,N), EVEC(LDEVEC,N)

!                                 Set values of A

!

!                                 A = ((1, 0)  (  1,-7i)  ( 0,- i))

!                                     ((1,7i)  (  5,  0)  (10,-3i))

!                                     ((0, i)  ( 10, 3i)  (-2,  0))

!

      DATA A/(1.0,0.0), (1.0,7.0), (0.0,1.0), (1.0,-7.0), (5.0,0.0), &

          (10.0, 3.0), (0.0,-1.0), (10.0,-3.0), (-2.0,0.0)/

!

!                                 Find eigenvalues and vectors of A

      CALL EVCHF (A, EVAL, EVEC)

!                                 Compute performance index

      PI = EPIHF(N,A,EVAL,EVEC)

!                                 Print results

      CALL UMACH (2, NOUT)

      CALL WRRRN ('EVAL', EVAL, 1, N, 1)

      CALL WRCRN ('EVEC', EVEC)

      WRITE (NOUT,'(/,A,F6.3)') ' Performance index = ', PI

      END

Output

          EVAL

      1       2       3

  15.38  -10.63   -0.75


                            EVEC

                    1                  2                  3

 1  ( 0.0631,-0.4075)  (-0.0598,-0.3117)  ( 0.8539, 0.0000)

 2  ( 0.7703, 0.0000)  (-0.5939, 0.1841)  (-0.0313,-0.1380)

 3  ( 0.4668, 0.1366)  ( 0.7160, 0.0000)  ( 0.0808,-0.4942)

 

 Performance index =  0.093

有头绪了，先吃饱再说。早饭还没吃那

[em2][em2][em2][em2][em11][em11][em11][em11][em20][em20][em20][em20][em68][em68][em68][em68][em68][em68][em72][em72][em72][em72][em72]

A — Complex Hermitian matrix of order N.   (Input) 
Only the upper triangle is used.
A是N阶矩阵，A(N,N)，下面那句说只用上三角矩阵，我看你给我的例子好像任何复数矩阵都可以用啊。

EVAL — Real vector of length N containing the eigenvalues of A in decreasing order of magnitude.   (Output)
长度为N，包含以降序大小排列的A的特征值的实矢量

EVEC — Complex matrix of order N.   (Output) 
The J-th eigenvector, corresponding to EVAL(J), is stored in the J-th column. Each vector is normalized to have Euclidean length equal to the value one.

这句应该描述特征向量吧。N阶复数矩阵。
第J个特征向量，对应于EVAL(J)，保存在第J列。

Optional Arguments 
N — Order of the matrix A.   (Input)
Default: N = size (A,2).

N是矩阵A的阶数，2维数组的长度。
下面这两句话说的是什么呢？是fortran77的语句吧？上面的是90的
LDA — Leading dimension of A exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDA = size (A,1).

Leading dimension ，是主维的意思吧？就是说一个2维数组（5,8），8就是主维啦，长嘛。LDA取矩阵A最长的维数，是用来给出矩阵特征值大小的吗？

LDEVEC — Leading dimension of EVEC exactly as specified in the dimension statement in the calling program.   (Input)
Default: LDEVEC = size (EVEC,1).

LDEVEC是用来规定矩阵特征向量维数大小的吧？

我求的矩阵是厄米矩阵。大姐给我那程序在我机器上运行出现点小问题，我还是把例子拿出来说下吧：
USE IMSL_libraries ！Error: Error in opening the Library module file.   
      IMPLICIT   NONE

!     Declare variables变量说明
      INTEGER    LDA, LDEVEC, N
      PARAMETER  (N=3, LDA=N, LDEVEC=N)！设置矩阵及其特征值、特征向量的维数
      INTEGER    NOUT ！？干嘛用呢
      REAL       EVAL(N), PI ！矩阵特征值，不知道PI怎么用
      COMPLEX    A(LDA,N), EVEC(LDEVEC,N) ！矩阵A、特征向量矩阵

!           Set values of A 设置A的值

!                                 A = ((1, 0)  (  1,-7i)  ( 0,- i))

!                                     ((1,7i)  (  5,  0)  (10,-3i))

!                                     ((0, i)  ( 10, 3i)  (-2,  0))

!

      DATA A/(1.0,0.0), (1.0,7.0), (0.0,1.0), (1.0,-7.0), (5.0,0.0), &

          (10.0, 3.0), (0.0,-1.0), (10.0,-3.0), (-2.0,0.0)/
!      Find eigenvalues and vectors of A找A的本征值和本征矢量
      CALL EVCHF (A, EVAL, EVEC)

!    Compute performance index？这下面干嘛的？

      PI = EPIHF(N,A,EVAL,EVEC)

!     Print results输出结果，下面的语句都是干嘛用啊？

      CALL UMACH (2, NOUT)

      CALL WRRRN ('EVAL', EVAL, 1, N, 1)

      CALL WRCRN ('EVEC', EVEC)

      WRITE (NOUT,'(/,A,F6.3)') ' Performance index = ', PI

      END

两个错误：
Error: Error in opening the Library module file.   [IMSL_LIBRARIES]
 USE IMSL_libraries 

Error: This name does not have a type, and must have an explicit type.   [EPIHF]
          PI = EPIHF(N,A,EVAL,EVEC)
下面那错误还好所，上面那个错误跟我的fortran版本有关吗？我的是带imsl的，6.6啊

修改代码如下，加入eig和EVCHF对比：
!求矩阵特征值和特征向量
 eigenvalue=eig(c,w=eigenvector) 
 write(*,*) "this is reasult form eig function!"
do i = 1,4
 write(*,'(i4,9e20.10)')i,eigenvalue(i),eigenvector(:,i)
enddo 
 
 call hmatr(c)
 call EVCHF(4, c, 4, eigenvalue, eigenvector, 4)
 write(*,*) "this is reasult form EVCHF subroutine!"
 do i = 1,4
 write(*,'(i4,9e20.10)')i, eigenvalue(i),eigenvector(:,i)
enddo
结果如下：
the reasult format for :  i,eigenvalue(i),eigenvector(:,i)
 this is reasult form eig function!
   1    0.1618034005E+01    0.3717480302E+00    0.0000000000E+00    0.4772028029
E+00   -0.3661706150E+00    0.4772028029E+00   -0.3661706150E+00    0.2628655434
E+00   -0.2628655732E+00
   2   -0.1618033886E+01   -0.3717480004E+00    0.0000000000E+00    0.4772028029
E+00   -0.3661706150E+00   -0.4772028029E+00    0.3661706150E+00    0.2628655434
E+00   -0.2628655732E+00
   3   -0.6180340052E+00    0.6015009880E+00    0.0000000000E+00   -0.2949275374
E+00    0.2263058573E+00   -0.2949275076E+00    0.2263058424E+00    0.4253253639
E+00   -0.4253254235E+00
   4    0.6180339456E+00    0.6015009880E+00    0.0000000000E+00    0.2949275374
E+00   -0.2263058573E+00   -0.2949275076E+00    0.2263058424E+00   -0.4253254235
E+00    0.4253254533E+00
 this is reasult form EVCHF subroutine!
   1   -0.1618034005E+01   -0.2949274480E+00   -0.2263058126E+00    0.6015009880
E+00    0.0000000000E+00   -0.6015008092E+00    0.1562612262E-07    0.3685676754
E+00   -0.4852283746E-01
   2    0.1618033886E+01    0.2949274778E+00    0.2263058424E+00    0.6015009880
E+00    0.0000000000E+00    0.6015009284E+00   -0.2662865306E-07    0.3685675859
E+00   -0.4852282628E-01
   3    0.6180338860E+00    0.6015009880E+00    0.0000000000E+00    0.2949275374
E+00   -0.2263058722E+00   -0.2949275076E+00    0.2263058573E+00   -0.4253254235
E+00    0.4253254235E+00
   4   -0.6180338264E+00    0.6015010476E+00    0.0000000000E+00   -0.2949275076
E+00    0.2263058424E+00   -0.2949275076E+00    0.2263058573E+00    0.4253253639
E+00   -0.4253253639E+00
Press any key to continue

可以看出，其结果是：特征值基本一致，但对应的特征向量在0.1618033886E+01这个特征值上结果完全不同。

于是继续测试H*EV=EL*EV
!记录矩阵与特征向量的乘积
call hmatr(c)
open(12,file='ceshi1.txt')
   write(12,*) matmul(c,eigenvector(:,4))
 close(12)
!记录特征值与特征向量的乘积
 open(13,file='ceshi2.txt')
   write(13,*) eigenvalue(4)*eigenvector(:,4)
 close(13)
!此时的EV为-0.618...，对应特征向量为(0.6015009880E+00 ,0.0000000000E+00),( 0.2949275374
E+00, -0.2263058722E+00),(-0.2949275076E+00,0.2263058573E+00),(-0.4253254235
E+00,0.4253254235E+00)
结果完全一致，于是测试特征值为0.1618..e+1:
call hmatr(c)
open(12,file='ceshi1.txt')
   write(12,*) matmul(c,eigenvector(:,2))
 close(12)
!记录特征值与特征向量的乘积
 open(13,file='ceshi2.txt')
   write(13,*) eigenvalue(2)*eigenvector(:,2)
 close(13)

结果也完全一致，于是继续对比eig函数产生的特征值和特征向量，结果完全一致。
归根结底：排序了特征值，而对应特征向量没有排序出了问题。