14
$\begingroup$

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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
7
$\begingroup$

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).

Improve this answer
$\endgroup$
6
  • $\begingroup$ What would be your answer if I was working on a Shared Memory system? $\endgroup$
    – Inquest
    Feb 3, 2012 at 19:11
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jed Brown
    Feb 3, 2012 at 20:19
  • $\begingroup$ 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? $\endgroup$
    – Inquest
    Feb 3, 2012 at 20:31
  • $\begingroup$ 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. $\endgroup$
    – Jed Brown
    Feb 3, 2012 at 20:47
  • $\begingroup$ Flexibility in terms of? Distributed memory? $\endgroup$
    – Inquest
    Feb 3, 2012 at 20:57
5
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Inquest
    Feb 4, 2012 at 15:20
  • $\begingroup$ 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. $\endgroup$
    – Jack Poulson
    Feb 4, 2012 at 15:35
  • $\begingroup$ 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? $\endgroup$
    – Inquest
    Feb 4, 2012 at 15:37
  • $\begingroup$ 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. $\endgroup$
    – Jack Poulson
    Feb 4, 2012 at 15:42
  • $\begingroup$ Also, yes, the fact that they perform roughly the same amount of flops as memops is why they are avoided as much as possible. $\endgroup$
    – Jack Poulson
    Feb 4, 2012 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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