2: /* 
  3:    Provides an interface to the Spooles serial sparse solver
  4: */

  6: #include <../src/mat/impls/aij/seq/spooles/spooles.h>

  8: extern PetscErrorCode MatDestroy_SeqSBAIJ(Mat);

 12: PetscErrorCode MatDestroy_SeqSBAIJSpooles(Mat A)
 13: {
 14:   Mat_Spooles    *lu = (Mat_Spooles*)A->spptr;

 18:   if (lu && lu->CleanUpSpooles) {
 19:     FrontMtx_free(lu->frontmtx);
 20:     IV_free(lu->newToOldIV);
 21:     IV_free(lu->oldToNewIV);
 22:     InpMtx_free(lu->mtxA);
 23:     ETree_free(lu->frontETree);
 24:     IVL_free(lu->symbfacIVL);
 25:     SubMtxManager_free(lu->mtxmanager);
 26:     Graph_free(lu->graph);
 27:   }
 28:   PetscFree(A->spptr);
 29:   MatDestroy_SeqSBAIJ(A);
 30:   return(0);
 31: }

 33: #if !defined(PETSC_USE_COMPLEX)
 34: /* 
 35:   input:
 36:    F:                 numeric factor
 37:   output:
 38:    nneg, nzero, npos: matrix inertia 
 39: */

 43: PetscErrorCode MatGetInertia_SeqSBAIJSpooles(Mat F,int *nneg,int *nzero,int *npos)
 44: {
 45:   Mat_Spooles *lu = (Mat_Spooles*)F->spptr;
 46:   int         neg,zero,pos;

 49:   FrontMtx_inertia(lu->frontmtx, &neg, &zero, &pos);
 50:   if(nneg)  *nneg  = neg;
 51:   if(nzero) *nzero = zero;
 52:   if(npos)  *npos  = pos;
 53:   return(0);
 54: }
 55: #endif /* !defined(PETSC_USE_COMPLEX) */

 57: /* Note the Petsc r permutation is ignored */
 60: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJSpooles(Mat B,Mat A,IS r,const MatFactorInfo *info)
 61: {
 63:   B->ops->choleskyfactornumeric  = MatFactorNumeric_SeqSpooles;
 64: #if !defined(PETSC_USE_COMPLEX)
 65:   B->ops->getinertia             = MatGetInertia_SeqSBAIJSpooles;
 66: #endif

 68:   return(0);
 69: }

 71: EXTERN_C_BEGIN
 74: PetscErrorCode MatFactorGetSolverPackage_seqsbaij_spooles(Mat A,const MatSolverPackage *type)
 75: {
 77:   *type = MATSOLVERSPOOLES;
 78:   return(0);
 79: }
 80: EXTERN_C_END

 82: EXTERN_C_BEGIN
 85: PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat A,MatFactorType ftype,Mat *F)
 86: {
 87:   Mat            B;
 89:   Mat_Spooles    *lu;

 92:   if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only Cholesky factorization is support for Spooles from SBAIJ matrix");
 93:   MatCreate(((PetscObject)A)->comm,&B);
 94:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
 95:   MatSetType(B,((PetscObject)A)->type_name);
 96:   MatSeqSBAIJSetPreallocation(B,1,PETSC_NULL,PETSC_NULL);
 97:   PetscNewLog(B,Mat_Spooles,&lu);
 98:   B->spptr = lu;
 99:   lu->options.pivotingflag  = SPOOLES_NO_PIVOTING;
100:   lu->options.symflag       = SPOOLES_SYMMETRIC;   /* default */
101:   lu->flg                   = DIFFERENT_NONZERO_PATTERN;
102:   lu->options.useQR         = PETSC_FALSE;

104:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJSpooles;
105:   B->ops->destroy                = MatDestroy_SeqSBAIJSpooles;
106:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqsbaij_spooles",MatFactorGetSolverPackage_seqsbaij_spooles);
107:   B->factortype                  = ftype;
108:   *F = B;
109:   return(0);
110: }
111: EXTERN_C_END