Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

I have a sparse matrix coming from the discretization of a 3D elliptic PDE. The matrix is banded with seven non-zeros diagonals. The sparsity pattern of the matrix looks like this (the actual matrix is much larger):

Spy

It is worth noting that it is positive definite, but not symmetric.

I would like to compute the eigenvalues, right eigenvectors and left eigenvectors. I need at least, for a given integer k, the k first eigenvectors corresponding to the k eigenvalues closest to the origin.

What is the most efficient way to do this for such a matrix?

share|improve this question

3 Answers 3

up vote 8 down vote accepted

The simple answer is that you would use inverse iteration (subspace or with deflation). This is basically the power method (repeatedly multiplying the matrix by a vector and normalizing, singling out the eigenvector corresponding to the eigenvalue of largest magnitude) applied to $A^{-1}$. Since you desire the $k$ closest to the origin, you need to use some kind of block method or orthogonalization, otherwise the power iteration converges to the same eigenvector each time.

Practically, these simple algorithms are too primitive to be useful (they would take too long). More sophisticated variations of these include the (block) Lanczos or Arnoldi or Jacobi-Davidson methods. Be careful to note the restriction on the matrix properties (e.g. classical Lanczos only works on Hermitian matrices). Note that since you want the smallest eigenvalues, you need to be able to apply $A^{-1}x$ in these iterations. So each iteration requires a sparse matrix solve. Your 3D spatial grid structure might be amenable for multigrid techniques to speed up the solve, or if it's of moderate size you can do an (incomplete) LU or Cholesky decomposition. You can also use these as a preconditioner for the solver, which typically the eigenvalue methods take as a parameter.

Practically, Matlab has eigs, and this free book has a decent listing of existing code.

share|improve this answer

The easiest way is to use Matlab, Maple or Mathematica if the matrix is small (up to maybe a few 1,000 rows). If you have something bigger, use packages such as ARPACK.

I suggest this because writing eigenvalue solvers from scratch is a most definitely non-trivial endeavor -- and it is completely unnecessary because there are great packages for this task around already!

share|improve this answer
    
The MATLAB eigs function provides a convenient interface to ARPACK and can handle quite large matrices. –  Bill Greene May 31 at 12:32

There is EIGEN C++ library (http://eigen.tuxfamily.org/)

share|improve this answer
1  
Welcome to SciComp! Link only answers are discouraged. Please elaborate. –  Geoff Oxberry May 30 at 5:54

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.