14

Wherever I have seen, PETSc tutorial/documents etc. say that it is useful for linear algebra and usually specifies that sparse systems will benefit. What about dense matrices? I am concerned about solving $Ax=b$ for dense $A$.

I have written my own code for CG and QMR in Fortran. The basic construct is to rip apart the pseudo code and add BLAS routines wherever possible (ddot, dnrm and dgemv) with a little self tuning. How will this compare to PETSc?

I know the best answer would be for me to try it myself but because of reasons of time and others, that is not possible.

Any help is much appreciated.

Geoff Oxberry
  • 30,394
  • 9
  • 64
  • 127
Inquest
  • 3,394
  • 3
  • 26
  • 44

2 Answers2

7

If you have dense matrices with structure (e.g. fast transforms, Schur complements, etc), PETSc could be useful. In these cases, you won't be assembling the full matrix.

For assembled dense systems, PETSc currently uses PLAPACK, but the matrix distribution in PETSc native format is not the best to minimize communication (for most operations). Jack Poulson, Matt Knepley, and I have discussed porting PETSc's dense linear algebra to use Elemental, Jack's more modern dense linear algebra library. It hasn't happened yet, but we will do it when we have time.

If want a full-featured dense linear algebra library, use Elemental. It is likely to always have a more complete API for those operations than PETSc. If you need to interact with sparse or unassembled systems, it likely makes sense to stick with PETSc and extract what you need to use Elemental for the dense parts (or help us write the PETSc interface to Elemental).

Jed Brown
  • 25,650
  • 3
  • 72
  • 130
  • What would be your answer if I was working on a Shared Memory system? – Inquest Feb 03 '12 at 19:11
  • 1
    You can use multiple MPI processes or multiple pthreads. I don't see much point in writing parallel software that only works with shared memory because usually the memory locality that you think about when writing for distributed memory improves the performance relative to all but the best threaded implementations. If you only want a "serial" API that uses threads internally for dense linear algebra, you can use a threaded BLAS. – Jed Brown Feb 03 '12 at 20:19
  • I want my codes to work on my Xeon 12 core workstation (I'm not looking at Clusters as of now). I am using Intel MKL for (threaded) BLAS. Would you still recommend I try PETSc? – Inquest Feb 03 '12 at 20:31
  • If all you want is BLAS, you are happy with it, and you have already written the code, just use it. If you want more flexibility and/or the opportunity to use distributed memory, you might want PETSc. – Jed Brown Feb 03 '12 at 20:47
  • Flexibility in terms of? Distributed memory? – Inquest Feb 03 '12 at 20:57
  • In terms of mixing sparse or unassembled/fast-transform parts with the dense computation, interactions between threads/processes, and general extensibility. – Jed Brown Feb 03 '12 at 21:20
5

It's important to realize that parallel dense linear algebra libraries usually focus on level 3 BLAS routines (routines that perform $O(n^3)$ work with $O(n^2)$ data) and higher-level functionality like factorizations and eigensolvers. They usually don't tune the level 1 and level 2 operations that you're referring to.

Since you mentioned that you are on a shared-memory system, I would have recommended libFLAME and/or PLASMA, but I do not think that either will be significantly faster than vendor threaded BLAS for your level 1 and level 2 operations.

Jed recommended Elemental, which I happen to develop, but I will again stress that level 1 and level 2 operations are not the main focus of parallel dense linear algebra libraries. I have honestly never benchmarked any of those routines.

Jack Poulson
  • 7,599
  • 32
  • 40
  • So, from what I get, every subsequent parallel dense linear algebra library will try to optimize operations such as eigenvector calculation, solution of system rather than BLAS.

    Also, Elemental looks really impressive. I will definitely give it a run when I can.

     – Inquest Feb 04 '12 at 15:20
  • They will typically optimize level 3 BLAS. The reason is that most people using the library will be doing large calculations which can usually be mapped to level 3 BLAS. – Jack Poulson Feb 04 '12 at 15:35
  • Is that because BLAS 1/2 can't get any better? (Maybe because of the surface-to-volume property?) O(N^2) data and O(N^2) computation? – Inquest Feb 04 '12 at 15:37
  • It's because BLAS 1 and 2 are generally lower-order terms in the computation. The entire 'game' of high-performance dense linear algebra is to squint at every operation in just the right way so that you can call xGEMM and friends for as much of the work as possible. – Jack Poulson Feb 04 '12 at 15:42
  • Also, yes, the fact that they perform roughly the same amount of flops as memops is why they are avoided as much as possible. – Jack Poulson Feb 04 '12 at 15:51
  • Can you please come on chat or something (If you are free)? In applications like Krylov Subspace methods, what choice does one have but to use BLAS 1/2 ? – Inquest Feb 04 '12 at 16:03