Actual source code: sbaijspooles.c

      1: #define PETSCMAT_DLL

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

     7:  #include src/mat/impls/aij/seq/spooles/spooles.h

     11: PetscErrorCode MatDestroy_SeqSBAIJSpooles(Mat A)
     12: {
     14:
     16:   /* SeqSBAIJ_Spooles isn't really the matrix that USES spooles, */
     17:   /* rather it is a factory class for creating a symmetric matrix that can */
     18:   /* invoke Spooles' sequential cholesky solver. */
     19:   /* As a result, we don't have to clean up the stuff set by spooles */
     20:   /* as in MatDestroy_SeqAIJ_Spooles. */
     21:   MatConvert_Spooles_Base(A,MATSEQSBAIJ,MAT_REUSE_MATRIX,&A);
     22:   (*A->ops->destroy)(A);
     23:   return(0);
     24: }

     28: PetscErrorCode MatAssemblyEnd_SeqSBAIJSpooles(Mat A,MatAssemblyType mode) \
     29: {
     31:   PetscInt       bs;
     32:   Mat_Spooles    *lu=(Mat_Spooles *)(A->spptr);

     35:   (*lu->MatAssemblyEnd)(A,mode);
     36:   MatGetBlockSize(A,&bs);
     37:   if (bs > 1) SETERRQ1(PETSC_ERR_SUP,"Block size %D not supported by Spooles",bs);
     38:   lu->MatCholeskyFactorSymbolic  = A->ops->choleskyfactorsymbolic;
     39:   A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJSpooles;
     40:   return(0);
     41: }

     43: /*
     44:   input:
     45:    F:                 numeric factor
     46:   output:
     47:    nneg, nzero, npos: matrix inertia
     48: */

     52: PetscErrorCode MatGetInertia_SeqSBAIJSpooles(Mat F,int *nneg,int *nzero,int *npos)
     53: {
     54:   Mat_Spooles *lu = (Mat_Spooles*)F->spptr;
     55:   int         neg,zero,pos;

     58:   FrontMtx_inertia(lu->frontmtx, &neg, &zero, &pos);
     59:   if(nneg)  *nneg  = neg;
     60:   if(nzero) *nzero = zero;
     61:   if(npos)  *npos  = pos;
     62:   return(0);
     63: }

     65: /* Note the Petsc r permutation is ignored */
     68: PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJSpooles(Mat A,IS r,MatFactorInfo *info,Mat *F)
     69: {
     70:   Mat         B;
     71:   Mat_Spooles *lu;
     73:   int m=A->m,n=A->n;

     76:   /* Create the factorization matrix */
     77:   MatCreate(A->comm,&B);
     78:   MatSetSizes(B,m,n,m,n);
     79:   MatSetType(B,A->type_name);
     80:   MatSeqSBAIJSetPreallocation(B,1,PETSC_NULL,PETSC_NULL);

     82:   B->ops->choleskyfactornumeric  = MatFactorNumeric_SeqAIJSpooles;
     83:   B->ops->getinertia             = MatGetInertia_SeqSBAIJSpooles;
     84:   B->factor                      = FACTOR_CHOLESKY;

     86:   lu                        = (Mat_Spooles *)(B->spptr);
     87:   lu->options.pivotingflag  = SPOOLES_NO_PIVOTING;
     88:   lu->options.symflag       = SPOOLES_SYMMETRIC;   /* default */
     89:   lu->flg                   = DIFFERENT_NONZERO_PATTERN;
     90:   lu->options.useQR         = PETSC_FALSE;

     92:   *F = B;
     93:   return(0);
     94: }

     99: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqSBAIJSpooles(Mat A,MatType type,MatReuse reuse,Mat *newmat) {
    100:   /* This routine is only called to convert a MATSEQSBAIJ matrix */
    101:   /* to a MATSEQSBAIJSPOOLES matrix, so we will ignore 'MatType type'. */
    103:   Mat         B=*newmat;
    104:   Mat_Spooles *lu;

    107:   if (reuse == MAT_INITIAL_MATRIX) {
    108:     /* This routine is inherited, so we know the type is correct. */
    109:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    110:   }

    112:   PetscNew(Mat_Spooles,&lu);
    113:   B->spptr                       = (void*)lu;

    115:   lu->basetype                   = MATSEQSBAIJ;
    116:   lu->CleanUpSpooles             = PETSC_FALSE;
    117:   lu->MatDuplicate               = A->ops->duplicate;
    118:   lu->MatCholeskyFactorSymbolic  = A->ops->choleskyfactorsymbolic;
    119:   lu->MatLUFactorSymbolic        = A->ops->lufactorsymbolic;
    120:   lu->MatView                    = A->ops->view;
    121:   lu->MatAssemblyEnd             = A->ops->assemblyend;
    122:   lu->MatDestroy                 = A->ops->destroy;

    124:   B->ops->duplicate              = MatDuplicate_Spooles;
    125:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJSpooles;
    126:   B->ops->assemblyend            = MatAssemblyEnd_SeqSBAIJSpooles;
    127:   B->ops->destroy                = MatDestroy_SeqSBAIJSpooles;
    128:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaijspooles_seqsbaij_C",
    129:                                            "MatConvert_Spooles_Base",MatConvert_Spooles_Base);
    130:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqsbaijspooles_C",
    131:                                            "MatConvert_SeqSBAIJ_SeqSBAIJSpooles",MatConvert_SeqSBAIJ_SeqSBAIJSpooles);
    132:   PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJSPOOLES);
    133:   *newmat = B;
    134:   return(0);
    135: }

    138: /*MC
    139:   MATSEQSBAIJSPOOLES - MATSEQSBAIJSPOOLES = "seqsbaijspooles" - A matrix type providing direct solvers (Cholesky) for sequential symmetric
    140:   matrices via the external package Spooles.

    142:   If Spooles is installed (see the manual for
    143:   instructions on how to declare the existence of external packages),
    144:   a matrix type can be constructed which invokes Spooles solvers.
    145:   After calling MatCreate(...,A), simply call MatSetType(A,MATSEQSBAIJSPOOLES).

    147:   This matrix inherits from MATSEQSBAIJ.  As a result, MatSeqSBAIJSetPreallocation is
    148:   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from
    149:   the MATSEQSBAIJ type without data copy.

    151:   Options Database Keys:
    152: + -mat_type seqsbaijspooles - sets the matrix type to seqsbaijspooles during calls to MatSetFromOptions()
    153: . -mat_spooles_tau <tau> - upper bound on the magnitude of the largest element in L or U
    154: . -mat_spooles_seed <seed> - random number seed used for ordering
    155: . -mat_spooles_msglvl <msglvl> - message output level
    156: . -mat_spooles_ordering <BestOfNDandMS,MMD,MS,ND> - ordering used
    157: . -mat_spooles_maxdomainsize <n> - maximum subgraph size used by Spooles orderings
    158: . -mat_spooles_maxzeros <n> - maximum number of zeros inside a supernode
    159: . -mat_spooles_maxsize <n> - maximum size of a supernode
    160: . -mat_spooles_FrontMtxInfo <true,fase> - print Spooles information about the computed factorization
    161: . -mat_spooles_symmetryflag <0,1,2> - 0: SPOOLES_SYMMETRIC, 1: SPOOLES_HERMITIAN, 2: SPOOLES_NONSYMMETRIC
    162: . -mat_spooles_patchAndGoFlag <0,1,2> - 0: no patch, 1: use PatchAndGo strategy 1, 2: use PatchAndGo strategy 2
    163: . -mat_spooles_toosmall <dt> - drop tolerance for PatchAndGo strategy 1
    164: . -mat_spooles_storeids <bool integer> - if nonzero, stores row and col numbers where patches were applied in an IV object
    165: . -mat_spooles_fudge <delta> - fudge factor for rescaling diagonals with PatchAndGo strategy 2
    166: - -mat_spooles_storevalues <bool integer> - if nonzero and PatchAndGo strategy 2 is used, store change in diagonal value in a DV object

    168:    Level: beginner

    170: .seealso: MATMPISBAIJSPOOLES, MATSEQAIJSPOOLES, MATMPIAIJSPOOLES, PCCHOLESKY
    171: M*/

    176: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJSpooles(Mat A)
    177: {

    181:   /* Change type name before calling MatSetType to force proper construction of SeqSBAIJ */
    182:   /*   and SeqSBAIJSpooles types */
    183:   PetscObjectChangeTypeName((PetscObject)A,MATSEQSBAIJSPOOLES);
    184:   MatSetType(A,MATSEQSBAIJ);
    185:   MatConvert_SeqSBAIJ_SeqSBAIJSpooles(A,MATSEQSBAIJSPOOLES,MAT_REUSE_MATRIX,&A);
    186:   return(0);
    187: }

    190: /*MC
    191:   MATSBAIJSPOOLES - MATSBAIJSPOOLES = "sbaijspooles" - A matrix type providing direct solvers (Cholesky) for sequential and parallel symmetric matrices via the external package Spooles.

    193:   If Spooles is installed (see the manual for
    194:   instructions on how to declare the existence of external packages),
    195:   a matrix type can be constructed which invokes Spooles solvers.
    196:   After calling MatCreate(...,A), simply call MatSetType(A,MATSBAIJSPOOLES).

    198:   This matrix inherits from MATSBAIJ.  As a result, MatSeqSBAIJSetPreallocation and MatMPISBAIJSetPreallocation are
    199:   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from
    200:   the MATSBAIJ type without data copy.

    202:   Options Database Keys:
    203: + -mat_type sbaijspooles - sets the matrix type to sbaijspooles during calls to MatSetFromOptions()
    204: . -mat_spooles_tau <tau> - upper bound on the magnitude of the largest element in L or U
    205: . -mat_spooles_seed <seed> - random number seed used for ordering
    206: . -mat_spooles_msglvl <msglvl> - message output level
    207: . -mat_spooles_ordering <BestOfNDandMS,MMD,MS,ND> - ordering used
    208: . -mat_spooles_maxdomainsize <n> - maximum subgraph size used by Spooles orderings
    209: . -mat_spooles_maxzeros <n> - maximum number of zeros inside a supernode
    210: . -mat_spooles_maxsize <n> - maximum size of a supernode
    211: . -mat_spooles_FrontMtxInfo <true,fase> - print Spooles information about the computed factorization
    212: . -mat_spooles_symmetryflag <0,1,2> - 0: SPOOLES_SYMMETRIC, 1: SPOOLES_HERMITIAN, 2: SPOOLES_NONSYMMETRIC
    213: . -mat_spooles_patchAndGoFlag <0,1,2> - 0: no patch, 1: use PatchAndGo strategy 1, 2: use PatchAndGo strategy 2
    214: . -mat_spooles_toosmall <dt> - drop tolerance for PatchAndGo strategy 1
    215: . -mat_spooles_storeids <bool integer> - if nonzero, stores row and col numbers where patches were applied in an IV object
    216: . -mat_spooles_fudge <delta> - fudge factor for rescaling diagonals with PatchAndGo strategy 2
    217: - -mat_spooles_storevalues <bool integer> - if nonzero and PatchAndGo strategy 2 is used, store change in diagonal value in a DV object

    219:    Level: beginner

    221: .seealso: MATMPISBAIJSPOOLES, MATSEQAIJSPOOLES, MATMPIAIJSPOOLES, PCCHOLESKY
    222: M*/

    227: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SBAIJSpooles(Mat A)
    228: {
    230:   int size;

    233:   /* Change type name before calling MatSetType to force proper construction of SeqSBAIJSpooles or MPISBAIJSpooles */
    234:   PetscObjectChangeTypeName((PetscObject)A,MATSBAIJSPOOLES);
    235:   MPI_Comm_size(A->comm,&size);
    236:   if (size == 1) {
    237:     MatSetType(A,MATSEQSBAIJ);
    238:     MatConvert_SeqSBAIJ_SeqSBAIJSpooles(A,MATSEQSBAIJSPOOLES,MAT_REUSE_MATRIX,&A);
    239:   } else {
    240:     MatSetType(A,MATMPISBAIJ);
    241:     MatConvert_MPISBAIJ_MPISBAIJSpooles(A,MATMPISBAIJSPOOLES,MAT_REUSE_MATRIX,&A);
    242:   }
    243:
    244:   return(0);
    245: }