Direct Eigenvalue solver, which retains the the order of eigenvectors

[alyflex]

Hi guys
I have a problem that I need some help with, I am looking for a direct eigenvalue solver algorithm.

The problem is that all the eigenvalue solvers I can find seems to reorder the final matrix after the size of the eigenvalues.

The matrix it shall calculate is very small (5-10), so speed should not be an issue.

I am writing this in Fortran, and have the lapack library available, but thus far I have found no routine in it which can do what I need.

 

[Simon_Tyler]

I don't think that there is a natural order to eigenvalues. As long as the eigenvalues and eigenvectors are provided in the same order, you can reconstruct the original matrix...

Can you give a simple example of what you want/expect to happen?

 

[alyflex]

What I'm actually doing:
I'm doing imaginary time propagation in density functional theory.
So after each iteration I have to othogonalize my orbitals, this is done by subspace othogonalization, in which I create the overlap matrix M_ij, to get the new orbitals I then solve the coresponding eigenvalue problem to this matrix.

A simple example of what I'm doing mathematically:
I get a matrix
A=[1 0 0;0 0.7 0.5; 0 0.5 0.8]
Now this tells me that the first orbital does not overlap with any of the two other orbitals.
However if I solve the eigenvalue problem, I get the following:
[V,D]=eig[A]
D=[0.065 0 0; 0 0.8; 0 0 1.545]
V=[0 1 0;0.7047 0 -0.7095;-0.7095 0 -0.7047]
Now when I look at it, it is clear that the one that does not mix is now the second orbital. But how do I in general make this relation? such that I can permute the columns back to their right place?

I hope that made it more clear

 

[Simon_Tyler]

Are you sure that example is correct?
The diagonal matrix of eigenvalues don't match what I calculate - see the Mathematica code below.
Maybe you just made a mistake copying the example?

In some matrix algorithms, you might also be returned information about pivoting. It could be that maybe the you are returned this data and didn't realize it. This could explain how rows (or columns) are changing places.

Note that pivoting won't affect the values of the eigenvalues - so does not account for for the difference in eigenvalues below...

[PLAIN]http://www.calerga.com/doc/LME_matr.htm#eig[/PLAIN] 
[PLAIN]http://www.calerga.com/doc/LME_lapk.htm#eig[/PLAIN] 

In[1]:= importMatrix[m_String] := ImportString[StringTrim[m, "["|"]"], 
          "Table", "FieldSeparators" -> {" "}, "LineSeparators" -> {";"}]

In[2]:= a = importMatrix["[1 0 0;0 0.7 0.5; 0 0.5 0.8]"]

Out[2]= {{1, 0, 0}, {0, 0.7, 0.5}, {0, 0.5, 0.8}}

In[3]:= d = importMatrix["[0.065 0 0; 0 0.8 0; 0 0 1.545]"]
        v = importMatrix["[0 1 0;0.7047 0 -0.7095;-0.7095 0 -0.7047]"]

Out[3]= {{0.065, 0, 0}, {0, 0.8, 0}, {0, 0, 1.545}}

Out[4]= {{0, 1, 0}, {0.7047, 0, -0.7095}, {-0.7095, 0, -0.7047}}

In[5]:= a.Transpose[v]
        Transpose[v].d

Out[5]= {{0., 0.7047, -0.7095}, {0.7, -0.35475, -0.35235}, {0.5, -0.5676, -0.56376}}

Out[6]= {{0., 0.56376, -1.09618}, {0.065, 0., 0.}, {0., -0.5676, -1.08876}}

In[7]:= {val, vec} = Eigensystem[a]

Out[7]= {{1.25249, 1., 0.247506}, 
 {{0., 0.671005, 0.741453}, {1., 0., 0.}, {0., 0.741453, -0.671005}}}

In[8]:= a.Transpose[vec] == Transpose[vec].DiagonalMatrix[val]

Out[8]= True

In[9]:= Table[a.vec[[i]] == val[[i]] vec[[i]], {i, 3}]

Out[9]= {True, True, True}

 

[AlephZero]

I don't quite understand what you are asking here, but any algorithm that finds eigenvalues and vectors will return the eigenvectors (i.e. the columns of the matrix) in the same order as the eigenvalues. Otherwise, you wouldn't know which vector corresponds to which value!

If you want to permute them into a different order, when you interchange eigenvalues e_i and e_j, just interchange columns i and j of the eigenvector matrix.

 

[DuncanM]

. . .

However if I solve the eigenvalue problem, I get the following:
[V,D]=eig[A]
D=[0.065 0 0; 0 0.8; 0 0 1.545]
V=[0 1 0;0.7047 0 -0.7095;-0.7095 0 -0.7047]

. . .

You seem to have access to MATLAB; why do you have to write a program in FORTRAN?

Are you sure you cannot find anything in the LAPACK libraries? They have just about everything there. And the eigensolvers usually print out eigenvectors in the same order as the eigenvalues.

Alternately, there is the http://www.netlib.org/eispack/" library, which is older. But it still works. And I know for sure there is a routine there that prints out eigenvector1 for eigenvalue1, eigenvector2 for eigenvalue2, etc.

In fact, there is also an online solver that is a Javascript translation of an EISPACK routine that you may find useful (perhaps to compare results):

http://www.akiti.ca/EigR12Solver.html"

(Keep in mind that some routines might normalize the eigenvectors; others may not.)

Hope this helps.