9

I have a problem where I need to find all positive (as in the eigenvalue is positive) eigenpairs of a small (usually smaller than 60x60) nonsymmetric matrix. I can stop calculating when the eigenvalue is smaller than a certain threshold. I know that the eigenvalues are real. Any suggestions on algorithms I could use to try to squeeze out the best performance? I have to do several thousand of these decompositions, so speed is important.

Thank you in advance.

EDIT: I need to do this on the GPU in shared memory. The matrices are also not necessarily the same size. I'm not aware of any libraries that do this at the moment. Suggestions of algorithms that would be well suited to the problem would be appreciated.

Aron Ahmadia
  • 6,951
  • 4
  • 34
  • 54
Kantoku
  • 91
  • 3
  • 1
    If I got it right, you have a CUDA kernel that computes thousands of small matrices in shared memory, and you are not willing to copy them to global memory. Before trying to give an answer, there are some points to clarify. In CUDA shared memory lifetime is bound to block lifetime: how many threads you have for each matrix to decompose? Is extreme performance really important? (How expected eigenvalue extraction times compare to matrix generation times?) Based on what argument you know that the eigensystem is real? Can the eigensystem be defective? – Stefano M Sep 02 '12 at 10:15
  • Hello Stefano and thank you for your comment. For now, I will have the closest multiple of the warp size to the dimension of the matrix I'd like to decompose. Matrix generation times vary a lot, and there are cases where matrix generation time is more expensive, but there are many situations where the matrix generation time is less than the decomposition. I know the eigenvalues are real because of the way the matrix is generated. I'd rather not go into the details here, since it would detract from the original question. Finally, yes, the system can be defective. – Kantoku Sep 02 '12 at 16:32

2 Answers2

3

Without doing a lot of search I recommend you to look at the MAGMA library. Freely available code with continuous support. NVIDIA recognized MAGMA as a "A Breakthrough in Solvers for Eigenvalue Problems".

There is also CULA library, which is generally commercial product, although recently it has been made free for academic usage (see details here).

Alexander
  • 1,111
  • 1
  • 8
  • 15
  • Thank you for your reply Alexander. I've looked into both libraries before, and as far as I know, the functions are called from the host and the memory needs to be in global memory. I believe the overhead would be too much to justify the use. All of these matrices are generated in shared memory, used in the kernel and then discarded. I'd like to keep them there without having to put them back into global memory. Even if I did push them there, there would still be the issue of calling many kernel functions from the host (albeit in multiple streams). – Kantoku Aug 30 '12 at 13:49
  • 1
    @Kantoku, yes, those libraries are more general and they store the whole matrix in the global memory. If your matrices are in the shared memory only one SM can work on them, doesn't it? The implementation of EVD thus should be quite straightforward. – Alexander Aug 30 '12 at 15:33
  • Yes I would imagine so, which is why I was fishing for algorithms that would be appropriate for the situation. I'm not overly familiar with non symmetric evd, so I was looking for suggestions. – Kantoku Aug 30 '12 at 15:48
  • @Kantoku (and Alexander). Nonsymmetric EVD's are far from straightforward, even in the sequential case. It is still an active area of research. – Jack Poulson Aug 30 '12 at 20:49
  • @JackPoulson Ah yes, you are right, but I (and I assume Alexander as well) meant that it would be straightforward to apply an established algorithm to the problem, considering there are many simplifications that can be made when we take the size and nature of the matrix into consideration. The problem is: which algorithm. – Kantoku Aug 30 '12 at 22:19
2

Use the functions in LAPACK, it's unlikely that you can beat them in your own implementation.

Wolfgang Bangerth
  • 55,373
  • 59
  • 119
  • Hi Wolfgang. Thanks for the answer, but I intend to implement this on a GPU using CUDA and for several thousand of these tiny matrices (where each block handles the decomposition of a single matrix), and the matrices aren't necessarily the same size, so implementing something myself that uses shared memory seems to be my only choice. Any idea what algorithm would be best suited for these types of matrices? P.S. Thanks for the deal.II lectures you gave at KAUST last semester. I enjoyed them :) – Kantoku Aug 29 '12 at 17:30
  • 2
    @Kantoku You should add these details in your question, otherwise it is misleading. – Alexander Aug 30 '12 at 08:18
  • @Alexander I've updated the question with more details. Thanks for the suggestion! – Kantoku Aug 30 '12 at 09:27
  • 1
    @Kantoku: GPUs are a bit beyond my realm but I'm sure there are libraries out there already that do what you want (and in fact I see that other answers already link to them). Glad to hear you liked my classes! – Wolfgang Bangerth Aug 30 '12 at 13:26