exporting patch:
# HG changeset patch
# User Paul Mullowney <paulm@txcorp.com>
# Date 1360003367 25200
# Node ID 6bd9465700d3a59249678f17387b3d2f59b08391
# Parent  2ccc7ff4e0471a48a660a80ec2eed85328333712
Adding interfaces to do cholesky and icc factorization from the cusparse class. Matrix construction is currently empty.

diff -r 2ccc7ff4e047 -r 6bd9465700d3 src/mat/impls/aij/seq/seqcusparse/aijcusparse.cu
--- a/src/mat/impls/aij/seq/seqcusparse/aijcusparse.cu	Mon Feb 04 08:00:07 2013 -0500
+++ b/src/mat/impls/aij/seq/seqcusparse/aijcusparse.cu	Mon Feb 04 11:42:47 2013 -0700
@@ -21,15 +21,21 @@
    and reuse it. */
 cudaStream_t theBodyStream=0;
 
+PetscErrorCode MatICCFactorSymbolic_SeqAIJCUSPARSE(Mat,Mat,IS,const MatFactorInfo*);
+PetscErrorCode MatCholeskyFactorSymbolic_SeqAIJCUSPARSE(Mat,Mat,IS,const MatFactorInfo*);
+PetscErrorCode MatCholeskyFactorNumeric_SeqAIJCUSPARSE(Mat,Mat,const MatFactorInfo*);
+PetscErrorCode MatSeqAIJCUSPARSEICCAnalysisAndCopyToGPU(Mat);
+
 PetscErrorCode MatILUFactorSymbolic_SeqAIJCUSPARSE(Mat,Mat,IS,IS,const MatFactorInfo*);
 PetscErrorCode MatLUFactorSymbolic_SeqAIJCUSPARSE(Mat,Mat,IS,IS,const MatFactorInfo*);
 PetscErrorCode MatLUFactorNumeric_SeqAIJCUSPARSE(Mat,Mat,const MatFactorInfo*);
+PetscErrorCode MatSeqAIJCUSPARSEILUAnalysisAndCopyToGPU(Mat);
+
 PetscErrorCode MatSolve_SeqAIJCUSPARSE(Mat,Vec,Vec);
 PetscErrorCode MatSolve_SeqAIJCUSPARSE_NaturalOrdering(Mat,Vec,Vec);
 PetscErrorCode MatSolveTranspose_SeqAIJCUSPARSE(Mat,Vec,Vec);
 PetscErrorCode MatSolveTranspose_SeqAIJCUSPARSE_NaturalOrdering(Mat,Vec,Vec);
 PetscErrorCode MatSetFromOptions_SeqAIJCUSPARSE(Mat);
-PetscErrorCode MatSeqAIJCUSPARSEAnalysisAndCopyToGPU(Mat);
 PetscErrorCode MatMult_SeqAIJCUSPARSE(Mat,Vec,Vec);
 PetscErrorCode MatMultAdd_SeqAIJCUSPARSE(Mat,Vec,Vec,Vec);
 PetscErrorCode MatMultTranspose_SeqAIJCUSPARSE(Mat,Vec,Vec);
@@ -70,6 +76,12 @@
 
     (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqAIJCUSPARSE;
     (*B)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqAIJCUSPARSE;
+  } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
+    ierr = MatSetType(*B,MATSEQAIJCUSPARSE);CHKERRQ(ierr);
+    ierr = MatSetFromOptions_SeqAIJCUSPARSE(*B);CHKERRQ(ierr);
+    ierr = PetscObjectComposeFunctionDynamic((PetscObject)(*B),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_cusparse",MatFactorGetSolverPackage_seqaij_cusparse);CHKERRQ(ierr);    
+    (*B)->ops->iccfactorsymbolic      = MatICCFactorSymbolic_SeqAIJCUSPARSE;
+    (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqAIJCUSPARSE;
   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported for CUSPARSE Matrix Types");
   (*B)->factortype = ftype;
   PetscFunctionReturn(0);
@@ -173,10 +185,8 @@
 PetscErrorCode MatILUFactorSymbolic_SeqAIJCUSPARSE(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
 {
   PetscErrorCode ierr;
-
   PetscFunctionBegin;
   ierr = MatILUFactorSymbolic_SeqAIJ(B,A,isrow,iscol,info);CHKERRQ(ierr);
-
   B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJCUSPARSE;
   PetscFunctionReturn(0);
 }
@@ -186,17 +196,39 @@
 PetscErrorCode MatLUFactorSymbolic_SeqAIJCUSPARSE(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
 {
   PetscErrorCode ierr;
-
   PetscFunctionBegin;
   ierr = MatLUFactorSymbolic_SeqAIJ(B,A,isrow,iscol,info);CHKERRQ(ierr);
-
   B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJCUSPARSE;
   PetscFunctionReturn(0);
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "MatSeqAIJCUSPARSEBuildLowerTriMatrix"
-PetscErrorCode MatSeqAIJCUSPARSEBuildLowerTriMatrix(Mat A)
+#define __FUNCT__ "MatICCFactorSymbolic_SeqAIJCUSPARSE"
+PetscErrorCode MatICCFactorSymbolic_SeqAIJCUSPARSE(Mat B,Mat A,IS perm,const MatFactorInfo *info)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = MatICCFactorSymbolic_SeqAIJ(B,A,perm,info);CHKERRQ(ierr);
+  B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJCUSPARSE;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MatCholeskyFactorSymbolic_SeqAIJCUSPARSE"
+PetscErrorCode MatCholeskyFactorSymbolic_SeqAIJCUSPARSE(Mat B,Mat A,IS perm,const MatFactorInfo *info)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = MatCholeskyFactorSymbolic_SeqAIJ(B,A,perm,info);CHKERRQ(ierr);
+  B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqAIJCUSPARSE;
+  PetscFunctionReturn(0);
+}
+
+
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSeqAIJCUSPARSEBuildILULowerTriMatrix"
+PetscErrorCode MatSeqAIJCUSPARSEBuildILULowerTriMatrix(Mat A)
 {
   Mat_SeqAIJ                   *a                  = (Mat_SeqAIJ*)A->data;
   PetscInt                     n                   = A->rmap->n;
@@ -266,8 +298,8 @@
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "MatSeqAIJCUSPARSEBuildUpperTriMatrix"
-PetscErrorCode MatSeqAIJCUSPARSEBuildUpperTriMatrix(Mat A)
+#define __FUNCT__ "MatSeqAIJCUSPARSEBuildLUUpperTriMatrix"
+PetscErrorCode MatSeqAIJCUSPARSEBuildILUUpperTriMatrix(Mat A)
 {
   Mat_SeqAIJ                   *a                  = (Mat_SeqAIJ*)A->data;
   PetscInt                     n                   = A->rmap->n;
@@ -333,8 +365,8 @@
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "MatSeqAIJCUSPARSEAnalysisAndCopyToGPU"
-PetscErrorCode MatSeqAIJCUSPARSEAnalysisAndCopyToGPU(Mat A)
+#define __FUNCT__ "MatSeqAIJCUSPARSEILUAnalysisAndCopyToGPU"
+PetscErrorCode MatSeqAIJCUSPARSEILUAnalysisAndCopyToGPU(Mat A)
 {
   PetscErrorCode               ierr;
   Mat_SeqAIJ                   *a                  = (Mat_SeqAIJ*)A->data;
@@ -345,8 +377,8 @@
   PetscInt                     n = A->rmap->n;
 
   PetscFunctionBegin;
-  ierr = MatSeqAIJCUSPARSEBuildLowerTriMatrix(A);CHKERRQ(ierr);
-  ierr = MatSeqAIJCUSPARSEBuildUpperTriMatrix(A);CHKERRQ(ierr);
+  ierr = MatSeqAIJCUSPARSEBuildILULowerTriMatrix(A);CHKERRQ(ierr);
+  ierr = MatSeqAIJCUSPARSEBuildILUUpperTriMatrix(A);CHKERRQ(ierr);
 
   cusparseTriFactors->tempvec = new CUSPARRAY;
   cusparseTriFactors->tempvec->resize(n);
@@ -370,6 +402,17 @@
   PetscFunctionReturn(0);
 }
 
+
+
+#undef __FUNCT__
+#define __FUNCT__ "MatSeqAIJCUSPARSEICCAnalysisAndCopyToGPU"
+PetscErrorCode MatSeqAIJCUSPARSEICCAnalysisAndCopyToGPU(Mat A)
+{
+  PetscFunctionBegin;
+  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"ICC and Cholesky Factorization not yet supported for CUSPARSE matrices");
+  PetscFunctionReturn(0);
+}
+
 #undef __FUNCT__
 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJCUSPARSE"
 PetscErrorCode MatLUFactorNumeric_SeqAIJCUSPARSE(Mat B,Mat A,const MatFactorInfo *info)
@@ -393,7 +436,35 @@
   }
 
   /* get the triangular factors */
-  ierr = MatSeqAIJCUSPARSEAnalysisAndCopyToGPU(B);CHKERRQ(ierr);
+  ierr = MatSeqAIJCUSPARSEILUAnalysisAndCopyToGPU(B);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "MatCholeskyFactorNumeric_SeqAIJCUSPARSE"
+PetscErrorCode MatCholeskyFactorNumeric_SeqAIJCUSPARSE(Mat B,Mat A,const MatFactorInfo *info)
+{
+  PetscErrorCode ierr;
+  Mat_SeqAIJ     *b    =(Mat_SeqAIJ*)B->data;
+  IS             isrow = b->row,iscol = b->col;
+  PetscBool      row_identity,col_identity;
+
+  PetscFunctionBegin;
+  ierr = MatCholeskyFactorNumeric_SeqAIJ(B,A,info);CHKERRQ(ierr);
+  /* determine which version of MatSolve needs to be used. */
+  ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
+  ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
+  if (row_identity && col_identity) {
+    B->ops->solve = MatSolve_SeqAIJCUSPARSE_NaturalOrdering;
+    B->ops->solvetranspose = MatSolveTranspose_SeqAIJCUSPARSE_NaturalOrdering;
+  } else {
+    B->ops->solve = MatSolve_SeqAIJCUSPARSE;
+    B->ops->solvetranspose = MatSolveTranspose_SeqAIJCUSPARSE;
+  }
+
+  /* get the triangular factors */
+  ierr = MatSeqAIJCUSPARSEICCAnalysisAndCopyToGPU(B);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }