[petsc-dev] patches for ICC/Cholesky preconditioner on GPUs

[Jed Brown]

Paul Mullowney <paulm at txcorp.com> writes:

> Hi,
>
> I've gotten GPU based icc preconditioners working with the following set 
> of patches. One uses ICC by loading an upper triangular matrix into the 
> cusparse class, and then setting the symmetry option, i.e.
>
>      ierr = 
> MatCreateSeqAIJCUSPARSE(comm,n,m,PETSC_NULL,num_entries_per_row,&A);CHKERRQ(ierr);
>      ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
>
> I've organized the implementations into 2 files:
> 1) icc-cholesky-cusparse-organization.patch : This patch has the basic 
> infrastructure for calling the symbolic and numeric factorizations (both 
> ICC and Cholesky). Currently these routines make scoping calls to the 
> implementations in aij/seq/aijfact.c

Can you write context like this into the commit messages? They are
pretty opaque right now, but they are more likely to be read in the
future than this email.

This patch series does not have any tests. Can you add tests so that we
can have a reasonable chance of not breaking it in the future?

We can review on the mailing list, but I think we would prefer to use
pull requests, mostly because comments there sort by context rather than
by author. Barry, Matt, and others, do you prefer the granularity of
email or bitbucket comments?

A few comments now:

> exporting patch:

This first line should not be in the patch file. It should start with the lines below:

> # HG changeset patch
> # User Paul Mullowney <paulm at txcorp.com>
> # Date 1362859983 25200
> # Node ID 6be0d0f06ad18794377015b5ec468e11105b58b1
> # Parent  6d60df3aadfe49f0bf640fa0a41fe3a423b88cf6
> Implementation of ICC in this patch.

Include the component (aijcusparse) in the commit message, otherwise we
don't have a clue from the commit message what is actually being done.

> diff -r 6d60df3aadfe -r 6be0d0f06ad1 src/mat/impls/aij/seq/seqcusparse/aijcusparse.cu
> --- a/src/mat/impls/aij/seq/seqcusparse/aijcusparse.cu	Mon Feb 04 11:42:47 2013 -0700
> +++ b/src/mat/impls/aij/seq/seqcusparse/aijcusparse.cu	Sat Mar 09 13:13:03 2013 -0700
> @@ -1,12 +1,13 @@
>  /*
> -    Defines the basic matrix operations for the AIJ (compressed row)
> +  Defines the basic matrix operations for the AIJ (compressed row)
>    matrix storage format.
>  */
>  
>  #include "petscconf.h"
>  PETSC_CUDA_EXTERN_C_BEGIN
>  #include "../src/mat/impls/aij/seq/aij.h"          /*I "petscmat.h" I*/
> -//#include "petscbt.h"
> +#include <../src/mat/impls/sbaij/seq/sbaij.h>
> +
>  #include "../src/vec/vec/impls/dvecimpl.h"
>  #include "petsc-private/vecimpl.h"
>  PETSC_CUDA_EXTERN_C_END

Not part of this patch, but PETSC_CUDA_EXTERN_C_BEGIN/END should be
replaced with appropriate use of PETSC_EXTERN (as shown below). We are
no longer mangling any public function names, so we can now compile
PETSc using --with-clanguage=C++ and then call it from plain C.

> @@ -53,12 +54,28 @@
>  EXTERN_C_BEGIN
>  extern PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*);
>  EXTERN_C_END

Use this instead:

PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*);

> @@ -404,12 +417,119 @@
>  
>  
>  
> +
> +

Isn't this a lot of blank lines?

> +#undef __FUNCT__
> +#define __FUNCT__ "MatSeqAIJCUSPARSEBuildICCTriMatrices"
> +PetscErrorCode MatSeqAIJCUSPARSEBuildICCTriMatrices(Mat A)

Please use 'static' if the function is really only used in this one
file. This identifier is more than 31 characters, thus too long for a
public symbol.

> +{
> +  Mat_SeqAIJ                   *a                  = (Mat_SeqAIJ*)A->data;
> +  Mat_SeqAIJCUSPARSETriFactors *cusparseTriFactors = (Mat_SeqAIJCUSPARSETriFactors*)A->spptr;
> +  GPU_Matrix_Ifc               * cusparseMatLo       = (GPU_Matrix_Ifc*)cusparseTriFactors->loTriFactorPtr;
> +  GPU_Matrix_Ifc               * cusparseMatUp       = (GPU_Matrix_Ifc*)cusparseTriFactors->upTriFactorPtr;

Either line up consecutive lines or don't, but don't go half way. PETSc
style guide says 'type *ptr', not 'type * ptr'.

> +  cusparseStatus_t             stat;
> +  PetscErrorCode               ierr;
> +  PetscInt                     *AiUp, *AjUp;
> +  PetscScalar                  *AAUp;
> +  PetscScalar                  *AALo;
> +  PetscInt          nzUpper=a->nz,n=A->rmap->n,i,offset,nz,j;
> +  Mat_SeqSBAIJ      *b = (Mat_SeqSBAIJ*)A->data;
> +  const PetscInt    *ai=b->i,*aj=b->j,*vj;
> +  const MatScalar   *aa=b->a,*v;

Again, line up or don't line up.

> +      /* set this flag ... for performance logging */
> +      ((Mat_SeqAIJCUSPARSETriFactors*)A->spptr)->isSymmOrHerm = PETSC_TRUE;

Does this have a MatSetOption() interface?

>  #undef __FUNCT__
>  #define __FUNCT__ "MatSeqAIJCUSPARSEICCAnalysisAndCopyToGPU"
>  PetscErrorCode MatSeqAIJCUSPARSEICCAnalysisAndCopyToGPU(Mat A)

This can have static linkage, right?

> @@ -706,8 +848,28 @@
>        ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);

[Different patch] PETSC_COMM_WORLD cannot be used outside of global
registration.