Actual source code: spooles.c

      1: #define PETCSMAT_DLL

      3: /*
      4:    Provides an interface to the Spooles serial sparse solver
      5: */
     6:  #include src/mat/impls/aij/seq/aij.h
     7:  #include src/mat/impls/sbaij/seq/sbaij.h
     8:  #include src/mat/impls/aij/seq/spooles/spooles.h

     13: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_Spooles_Base(Mat A,MatType type,MatReuse reuse,Mat *newmat) {
     14:   /* This routine is only called to convert an unfactored PETSc-Spooles matrix */
     15:   /* to its base PETSc type, so we will ignore 'MatType type'. */
     17:   Mat            B=*newmat;
     18:   Mat_Spooles    *lu=(Mat_Spooles*)A->spptr;
     19:   void           (*f)(void);

     22:   if (reuse == MAT_INITIAL_MATRIX) {
     23:     MatDuplicate(A,MAT_COPY_VALUES,&B);
     24:   }
     25:   /* Reset the stashed function pointers set by inherited routines */
     26:   B->ops->duplicate              = lu->MatDuplicate;
     27:   B->ops->choleskyfactorsymbolic = lu->MatCholeskyFactorSymbolic;
     28:   B->ops->lufactorsymbolic       = lu->MatLUFactorSymbolic;
     29:   B->ops->view                   = lu->MatView;
     30:   B->ops->assemblyend            = lu->MatAssemblyEnd;
     31:   B->ops->destroy                = lu->MatDestroy;

     33:   PetscObjectQueryFunction((PetscObject)B,"MatMPISBAIJSetPreallocation_C",&f);
     34:   if (f) {
     35:     PetscObjectComposeFunction((PetscObject)B,"MatMPISBAIJSetPreallocation_C","",(FCNVOID)lu->MatPreallocate);
     36:   }
     37:   PetscFree(lu);

     39:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijspooles_seqaij_C","",PETSC_NULL);
     40:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijspooles_C","",PETSC_NULL);
     41:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpiaijspooles_mpiaij_C","",PETSC_NULL);
     42:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpiaij_mpiaijspooles_C","",PETSC_NULL);
     43:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaijspooles_seqsbaij_C","",PETSC_NULL);
     44:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbaijspooles_C","",PETSC_NULL);
     45:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisbaijspooles_mpisbaij_C","",PETSC_NULL);
     46:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisbaij_mpisbaijspooles_C","",PETSC_NULL);

     48:   PetscObjectChangeTypeName((PetscObject)B,type);
     49:   *newmat = B;
     50:   return(0);
     51: }

     56: PetscErrorCode MatDestroy_SeqAIJSpooles(Mat A)
     57: {
     58:   Mat_Spooles    *lu = (Mat_Spooles*)A->spptr;
     60:
     62:   if (lu->CleanUpSpooles) {
     63:     FrontMtx_free(lu->frontmtx);
     64:     IV_free(lu->newToOldIV);
     65:     IV_free(lu->oldToNewIV);
     66:     InpMtx_free(lu->mtxA);
     67:     ETree_free(lu->frontETree);
     68:     IVL_free(lu->symbfacIVL);
     69:     SubMtxManager_free(lu->mtxmanager);
     70:     Graph_free(lu->graph);
     71:   }
     72:   MatConvert_Spooles_Base(A,lu->basetype,MAT_REUSE_MATRIX,&A);
     73:   (*A->ops->destroy)(A);
     74:   return(0);
     75: }

     79: PetscErrorCode MatSolve_SeqAIJSpooles(Mat A,Vec b,Vec x)
     80: {
     81:   Mat_Spooles      *lu = (Mat_Spooles*)A->spptr;
     82:   PetscScalar      *array;
     83:   DenseMtx         *mtxY, *mtxX ;
     84:   PetscErrorCode   ierr;
     85:   PetscInt         irow,neqns=A->n,nrow=A->m,*iv;
     86: #if defined(PETSC_USE_COMPLEX)
     87:   double           x_real,x_imag;
     88: #else
     89:   double           *entX;
     90: #endif

     93:   mtxY = DenseMtx_new();
     94:   DenseMtx_init(mtxY, lu->options.typeflag, 0, 0, nrow, 1, 1, nrow); /* column major */
     95:   VecGetArray(b,&array);

     97:   if (lu->options.useQR) {   /* copy b to mtxY */
     98:     for ( irow = 0 ; irow < nrow; irow++ )
     99: #if !defined(PETSC_USE_COMPLEX)
    100:       DenseMtx_setRealEntry(mtxY, irow, 0, *array++);
    101: #else
    102:       DenseMtx_setComplexEntry(mtxY, irow, 0, PetscRealPart(array[irow]), PetscImaginaryPart(array[irow]));
    103: #endif
    104:   } else {                   /* copy permuted b to mtxY */
    105:     iv = IV_entries(lu->oldToNewIV);
    106:     for ( irow = 0 ; irow < nrow; irow++ )
    107: #if !defined(PETSC_USE_COMPLEX)
    108:       DenseMtx_setRealEntry(mtxY, *iv++, 0, *array++);
    109: #else
    110:       DenseMtx_setComplexEntry(mtxY,*iv++,0,PetscRealPart(array[irow]),PetscImaginaryPart(array[irow]));
    111: #endif
    112:   }
    113:   VecRestoreArray(b,&array);

    115:   mtxX = DenseMtx_new();
    116:   DenseMtx_init(mtxX, lu->options.typeflag, 0, 0, neqns, 1, 1, neqns);
    117:   if (lu->options.useQR) {
    118:     FrontMtx_QR_solve(lu->frontmtx, lu->mtxA, mtxX, mtxY, lu->mtxmanager,
    119:                   lu->cpus, lu->options.msglvl, lu->options.msgFile);
    120:   } else {
    121:     FrontMtx_solve(lu->frontmtx, mtxX, mtxY, lu->mtxmanager,
    122:                  lu->cpus, lu->options.msglvl, lu->options.msgFile);
    123:   }
    124:   if ( lu->options.msglvl > 2 ) {
    125:     PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n right hand side matrix after permutation");
    126:     DenseMtx_writeForHumanEye(mtxY, lu->options.msgFile);
    127:     PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n solution matrix in new ordering");
    128:     DenseMtx_writeForHumanEye(mtxX, lu->options.msgFile);
    129:     fflush(lu->options.msgFile);
    130:   }

    132:   /* permute solution into original ordering, then copy to x */
    133:   DenseMtx_permuteRows(mtxX, lu->newToOldIV);
    134:   VecGetArray(x,&array);

    136: #if !defined(PETSC_USE_COMPLEX)
    137:   entX = DenseMtx_entries(mtxX);
    138:   DVcopy(neqns, array, entX);
    139: #else
    140:   for (irow=0; irow<nrow; irow++){
    141:     DenseMtx_complexEntry(mtxX,irow,0,&x_real,&x_imag);
    142:     array[irow] = x_real+x_imag*PETSC_i;
    143:   }
    144: #endif

    146:   VecRestoreArray(x,&array);
    147:
    148:   /* free memory */
    149:   DenseMtx_free(mtxX);
    150:   DenseMtx_free(mtxY);
    151:   return(0);
    152: }

    156: PetscErrorCode MatFactorNumeric_SeqAIJSpooles(Mat A,MatFactorInfo *info,Mat *F)
    157: {
    158:   Mat_Spooles        *lu = (Mat_Spooles*)(*F)->spptr;
    159:   ChvManager         *chvmanager ;
    160:   Chv                *rootchv ;
    161:   IVL                *adjIVL;
    162:   PetscErrorCode     ierr;
    163:   PetscInt           nz,nrow=A->m,irow,nedges,neqns=A->n,*ai,*aj,i,*diag=0,fierr;
    164:   PetscScalar        *av;
    165:   double             cputotal,facops;
    166: #if defined(PETSC_USE_COMPLEX)
    167:   PetscInt           nz_row,*aj_tmp;
    168:   PetscScalar        *av_tmp;
    169: #else
    170:   PetscInt           *ivec1,*ivec2,j;
    171:   double             *dvec;
    172: #endif
    173:   PetscTruth         isAIJ,isSeqAIJ;
    174:
    176:   if (lu->flg == DIFFERENT_NONZERO_PATTERN) { /* first numeric factorization */
    177:     (*F)->ops->solve   = MatSolve_SeqAIJSpooles;
    178:     (*F)->ops->destroy = MatDestroy_SeqAIJSpooles;
    179:     (*F)->assembled    = PETSC_TRUE;
    180:
    181:     /* set Spooles options */
    182:     SetSpoolesOptions(A, &lu->options);

    184:     lu->mtxA = InpMtx_new();
    185:   }

    187:   /* copy A to Spooles' InpMtx object */
    188:   PetscTypeCompare((PetscObject)A,MATSEQAIJSPOOLES,&isSeqAIJ);
    189:   PetscTypeCompare((PetscObject)A,MATAIJSPOOLES,&isAIJ);
    190:   if (isSeqAIJ || isAIJ){
    191:     Mat_SeqAIJ   *mat = (Mat_SeqAIJ*)A->data;
    192:     ai=mat->i; aj=mat->j; av=mat->a;
    193:     if (lu->options.symflag == SPOOLES_NONSYMMETRIC) {
    194:       nz=mat->nz;
    195:     } else { /* SPOOLES_SYMMETRIC || SPOOLES_HERMITIAN */
    196:       nz=(mat->nz + A->m)/2;
    197:       if (!mat->diag){
    198:         MatMarkDiagonal_SeqAIJ(A);
    199:       }
    200:       diag=mat->diag;
    201:     }
    202:   } else { /* A is SBAIJ */
    203:       Mat_SeqSBAIJ *mat = (Mat_SeqSBAIJ*)A->data;
    204:       ai=mat->i; aj=mat->j; av=mat->a;
    205:       nz=mat->nz;
    206:   }
    207:   InpMtx_init(lu->mtxA, INPMTX_BY_ROWS, lu->options.typeflag, nz, 0);
    208:
    209: #if defined(PETSC_USE_COMPLEX)
    210:     for (irow=0; irow<nrow; irow++) {
    211:       if ( lu->options.symflag == SPOOLES_NONSYMMETRIC || !isAIJ){
    212:         nz_row = ai[irow+1] - ai[irow];
    213:         aj_tmp = aj + ai[irow];
    214:         av_tmp = av + ai[irow];
    215:       } else {
    216:         nz_row = ai[irow+1] - diag[irow];
    217:         aj_tmp = aj + diag[irow];
    218:         av_tmp = av + diag[irow];
    219:       }
    220:       for (i=0; i<nz_row; i++){
    221:         InpMtx_inputComplexEntry(lu->mtxA, irow, *aj_tmp++,PetscRealPart(*av_tmp),PetscImaginaryPart(*av_tmp));
    222:         av_tmp++;
    223:       }
    224:     }
    225: #else
    226:     ivec1 = InpMtx_ivec1(lu->mtxA);
    227:     ivec2 = InpMtx_ivec2(lu->mtxA);
    228:     dvec  = InpMtx_dvec(lu->mtxA);
    229:     if ( lu->options.symflag == SPOOLES_NONSYMMETRIC || !isAIJ){
    230:       for (irow = 0; irow < nrow; irow++){
    231:         for (i = ai[irow]; i<ai[irow+1]; i++) ivec1[i] = irow;
    232:       }
    233:       IVcopy(nz, ivec2, aj);
    234:       DVcopy(nz, dvec, av);
    235:     } else {
    236:       nz = 0;
    237:       for (irow = 0; irow < nrow; irow++){
    238:         for (j = diag[irow]; j<ai[irow+1]; j++) {
    239:           ivec1[nz] = irow;
    240:           ivec2[nz] = aj[j];
    241:           dvec[nz]  = av[j];
    242:           nz++;
    243:         }
    244:       }
    245:     }
    246:     InpMtx_inputRealTriples(lu->mtxA, nz, ivec1, ivec2, dvec);
    247: #endif

    249:   InpMtx_changeStorageMode(lu->mtxA, INPMTX_BY_VECTORS);
    250:   if ( lu->options.msglvl > 0 ) {
    251:     printf("\n\n input matrix");
    252:     PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n input matrix");
    253:     InpMtx_writeForHumanEye(lu->mtxA, lu->options.msgFile);
    254:     fflush(lu->options.msgFile);
    255:   }

    257:   if ( lu->flg == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization */
    258:     /*---------------------------------------------------
    259:     find a low-fill ordering
    260:          (1) create the Graph object
    261:          (2) order the graph
    262:     -------------------------------------------------------*/
    263:     if (lu->options.useQR){
    264:       adjIVL = InpMtx_adjForATA(lu->mtxA);
    265:     } else {
    266:       adjIVL = InpMtx_fullAdjacency(lu->mtxA);
    267:     }
    268:     nedges = IVL_tsize(adjIVL);

    270:     lu->graph = Graph_new();
    271:     Graph_init2(lu->graph, 0, neqns, 0, nedges, neqns, nedges, adjIVL, NULL, NULL);
    272:     if ( lu->options.msglvl > 2 ) {
    273:       if (lu->options.useQR){
    274:         PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n graph of A^T A");
    275:       } else {
    276:         PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n graph of the input matrix");
    277:       }
    278:       Graph_writeForHumanEye(lu->graph, lu->options.msgFile);
    279:       fflush(lu->options.msgFile);
    280:     }

    282:     switch (lu->options.ordering) {
    283:     case 0:
    284:       lu->frontETree = orderViaBestOfNDandMS(lu->graph,
    285:                      lu->options.maxdomainsize, lu->options.maxzeros, lu->options.maxsize,
    286:                      lu->options.seed, lu->options.msglvl, lu->options.msgFile); break;
    287:     case 1:
    288:       lu->frontETree = orderViaMMD(lu->graph,lu->options.seed,lu->options.msglvl,lu->options.msgFile); break;
    289:     case 2:
    290:       lu->frontETree = orderViaMS(lu->graph, lu->options.maxdomainsize,
    291:                      lu->options.seed,lu->options.msglvl,lu->options.msgFile); break;
    292:     case 3:
    293:       lu->frontETree = orderViaND(lu->graph, lu->options.maxdomainsize,
    294:                      lu->options.seed,lu->options.msglvl,lu->options.msgFile); break;
    295:     default:
    296:       SETERRQ(PETSC_ERR_ARG_WRONG,"Unknown Spooles's ordering");
    297:     }

    299:     if ( lu->options.msglvl > 0 ) {
    300:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n front tree from ordering");
    301:       ETree_writeForHumanEye(lu->frontETree, lu->options.msgFile);
    302:       fflush(lu->options.msgFile);
    303:     }
    304:
    305:     /* get the permutation, permute the front tree */
    306:     lu->oldToNewIV = ETree_oldToNewVtxPerm(lu->frontETree);
    307:     lu->oldToNew   = IV_entries(lu->oldToNewIV);
    308:     lu->newToOldIV = ETree_newToOldVtxPerm(lu->frontETree);
    309:     if (!lu->options.useQR) ETree_permuteVertices(lu->frontETree, lu->oldToNewIV);

    311:     /* permute the matrix */
    312:     if (lu->options.useQR){
    313:       InpMtx_permute(lu->mtxA, NULL, lu->oldToNew);
    314:     } else {
    315:       InpMtx_permute(lu->mtxA, lu->oldToNew, lu->oldToNew);
    316:       if ( lu->options.symflag == SPOOLES_SYMMETRIC) {
    317:         InpMtx_mapToUpperTriangle(lu->mtxA);
    318:       }
    319: #if defined(PETSC_USE_COMPLEX)
    320:       if ( lu->options.symflag == SPOOLES_HERMITIAN ) {
    321:         InpMtx_mapToUpperTriangleH(lu->mtxA);
    322:       }
    323: #endif
    324:       InpMtx_changeCoordType(lu->mtxA, INPMTX_BY_CHEVRONS);
    325:     }
    326:     InpMtx_changeStorageMode(lu->mtxA, INPMTX_BY_VECTORS);

    328:     /* get symbolic factorization */
    329:     if (lu->options.useQR){
    330:       lu->symbfacIVL = SymbFac_initFromGraph(lu->frontETree, lu->graph);
    331:       IVL_overwrite(lu->symbfacIVL, lu->oldToNewIV);
    332:       IVL_sortUp(lu->symbfacIVL);
    333:       ETree_permuteVertices(lu->frontETree, lu->oldToNewIV);
    334:     } else {
    335:       lu->symbfacIVL = SymbFac_initFromInpMtx(lu->frontETree, lu->mtxA);
    336:     }
    337:     if ( lu->options.msglvl > 2 ) {
    338:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n old-to-new permutation vector");
    339:       IV_writeForHumanEye(lu->oldToNewIV, lu->options.msgFile);
    340:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n new-to-old permutation vector");
    341:       IV_writeForHumanEye(lu->newToOldIV, lu->options.msgFile);
    342:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n front tree after permutation");
    343:       ETree_writeForHumanEye(lu->frontETree, lu->options.msgFile);
    344:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n input matrix after permutation");
    345:       InpMtx_writeForHumanEye(lu->mtxA, lu->options.msgFile);
    346:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n symbolic factorization");
    347:       IVL_writeForHumanEye(lu->symbfacIVL, lu->options.msgFile);
    348:       fflush(lu->options.msgFile);
    349:     }

    351:     lu->frontmtx   = FrontMtx_new();
    352:     lu->mtxmanager = SubMtxManager_new();
    353:     SubMtxManager_init(lu->mtxmanager, NO_LOCK, 0);

    355:   } else { /* new num factorization using previously computed symbolic factor */

    357:     if (lu->options.pivotingflag) { /* different FrontMtx is required */
    358:       FrontMtx_free(lu->frontmtx);
    359:       lu->frontmtx   = FrontMtx_new();
    360:     } else {
    361:       FrontMtx_clearData (lu->frontmtx);
    362:     }

    364:     SubMtxManager_free(lu->mtxmanager);
    365:     lu->mtxmanager = SubMtxManager_new();
    366:     SubMtxManager_init(lu->mtxmanager, NO_LOCK, 0);

    368:     /* permute mtxA */
    369:     if (lu->options.useQR){
    370:       InpMtx_permute(lu->mtxA, NULL, lu->oldToNew);
    371:     } else {
    372:       InpMtx_permute(lu->mtxA, lu->oldToNew, lu->oldToNew);
    373:       if ( lu->options.symflag == SPOOLES_SYMMETRIC ) {
    374:         InpMtx_mapToUpperTriangle(lu->mtxA);
    375:       }
    376:       InpMtx_changeCoordType(lu->mtxA, INPMTX_BY_CHEVRONS);
    377:     }
    378:     InpMtx_changeStorageMode(lu->mtxA, INPMTX_BY_VECTORS);
    379:     if ( lu->options.msglvl > 2 ) {
    380:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n input matrix after permutation");
    381:       InpMtx_writeForHumanEye(lu->mtxA, lu->options.msgFile);
    382:     }
    383:   } /* end of if( lu->flg == DIFFERENT_NONZERO_PATTERN) */
    384:
    385:   if (lu->options.useQR){
    386:     FrontMtx_init(lu->frontmtx, lu->frontETree, lu->symbfacIVL, lu->options.typeflag,
    387:                  SPOOLES_SYMMETRIC, FRONTMTX_DENSE_FRONTS,
    388:                  SPOOLES_NO_PIVOTING, NO_LOCK, 0, NULL,
    389:                  lu->mtxmanager, lu->options.msglvl, lu->options.msgFile);
    390:   } else {
    391:     FrontMtx_init(lu->frontmtx, lu->frontETree, lu->symbfacIVL, lu->options.typeflag, lu->options.symflag,
    392:                 FRONTMTX_DENSE_FRONTS, lu->options.pivotingflag, NO_LOCK, 0, NULL,
    393:                 lu->mtxmanager, lu->options.msglvl, lu->options.msgFile);
    394:   }

    396:   if ( lu->options.symflag == SPOOLES_SYMMETRIC ) {  /* || SPOOLES_HERMITIAN ? */
    397:     if ( lu->options.patchAndGoFlag == 1 ) {
    398:       lu->frontmtx->patchinfo = PatchAndGoInfo_new();
    399:       PatchAndGoInfo_init(lu->frontmtx->patchinfo, 1, lu->options.toosmall, lu->options.fudge,
    400:                        lu->options.storeids, lu->options.storevalues);
    401:     } else if ( lu->options.patchAndGoFlag == 2 ) {
    402:       lu->frontmtx->patchinfo = PatchAndGoInfo_new();
    403:       PatchAndGoInfo_init(lu->frontmtx->patchinfo, 2, lu->options.toosmall, lu->options.fudge,
    404:                        lu->options.storeids, lu->options.storevalues);
    405:     }
    406:   }

    408:   /* numerical factorization */
    409:   chvmanager = ChvManager_new();
    410:   ChvManager_init(chvmanager, NO_LOCK, 1);
    411:   DVfill(10, lu->cpus, 0.0);
    412:   if (lu->options.useQR){
    413:     facops = 0.0 ;
    414:     FrontMtx_QR_factor(lu->frontmtx, lu->mtxA, chvmanager,
    415:                    lu->cpus, &facops, lu->options.msglvl, lu->options.msgFile);
    416:     if ( lu->options.msglvl > 1 ) {
    417:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n factor matrix");
    418:       PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n facops = %9.2f", facops);
    419:     }
    420:   } else {
    421:     IVfill(20, lu->stats, 0);
    422:     rootchv = FrontMtx_factorInpMtx(lu->frontmtx, lu->mtxA, lu->options.tau, 0.0,
    423:             chvmanager, &fierr, lu->cpus,lu->stats,lu->options.msglvl,lu->options.msgFile);
    424:     if (rootchv) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"\n matrix found to be singular");
    425:     if (fierr >= 0) SETERRQ1(PETSC_ERR_LIB,"\n error encountered at front %D", fierr);
    426:
    427:     if(lu->options.FrontMtxInfo){
    428:       PetscPrintf(PETSC_COMM_SELF,"\n %8d pivots, %8d pivot tests, %8d delayed rows and columns\n",lu->stats[0], lu->stats[1], lu->stats[2]);
    429:       cputotal = lu->cpus[8] ;
    430:       if ( cputotal > 0.0 ) {
    431:         PetscPrintf(PETSC_COMM_SELF,
    432:            "\n                               cpus   cpus/totaltime"
    433:            "\n    initialize fronts       %8.3f %6.2f"
    434:            "\n    load original entries   %8.3f %6.2f"
    435:            "\n    update fronts           %8.3f %6.2f"
    436:            "\n    assemble postponed data %8.3f %6.2f"
    437:            "\n    factor fronts           %8.3f %6.2f"
    438:            "\n    extract postponed data  %8.3f %6.2f"
    439:            "\n    store factor entries    %8.3f %6.2f"
    440:            "\n    miscellaneous           %8.3f %6.2f"
    441:            "\n    total time              %8.3f \n",
    442:            lu->cpus[0], 100.*lu->cpus[0]/cputotal,
    443:            lu->cpus[1], 100.*lu->cpus[1]/cputotal,
    444:            lu->cpus[2], 100.*lu->cpus[2]/cputotal,
    445:            lu->cpus[3], 100.*lu->cpus[3]/cputotal,
    446:            lu->cpus[4], 100.*lu->cpus[4]/cputotal,
    447:            lu->cpus[5], 100.*lu->cpus[5]/cputotal,
    448:            lu->cpus[6], 100.*lu->cpus[6]/cputotal,
    449:            lu->cpus[7], 100.*lu->cpus[7]/cputotal, cputotal);
    450:       }
    451:     }
    452:   }
    453:   ChvManager_free(chvmanager);

    455:   if ( lu->options.msglvl > 0 ) {
    456:     PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n factor matrix");
    457:     FrontMtx_writeForHumanEye(lu->frontmtx, lu->options.msgFile);
    458:     fflush(lu->options.msgFile);
    459:   }

    461:   if ( lu->options.symflag == SPOOLES_SYMMETRIC ) { /* || SPOOLES_HERMITIAN ? */
    462:     if ( lu->options.patchAndGoFlag == 1 ) {
    463:       if ( lu->frontmtx->patchinfo->fudgeIV != NULL ) {
    464:         if (lu->options.msglvl > 0 ){
    465:           PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n small pivots found at these locations");
    466:           IV_writeForHumanEye(lu->frontmtx->patchinfo->fudgeIV, lu->options.msgFile);
    467:         }
    468:       }
    469:       PatchAndGoInfo_free(lu->frontmtx->patchinfo);
    470:     } else if ( lu->options.patchAndGoFlag == 2 ) {
    471:       if (lu->options.msglvl > 0 ){
    472:         if ( lu->frontmtx->patchinfo->fudgeIV != NULL ) {
    473:           PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n small pivots found at these locations");
    474:           IV_writeForHumanEye(lu->frontmtx->patchinfo->fudgeIV, lu->options.msgFile);
    475:         }
    476:         if ( lu->frontmtx->patchinfo->fudgeDV != NULL ) {
    477:           PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n perturbations");
    478:           DV_writeForHumanEye(lu->frontmtx->patchinfo->fudgeDV, lu->options.msgFile);
    479:         }
    480:       }
    481:       PatchAndGoInfo_free(lu->frontmtx->patchinfo);
    482:     }
    483:   }

    485:   /* post-process the factorization */
    486:   FrontMtx_postProcess(lu->frontmtx, lu->options.msglvl, lu->options.msgFile);
    487:   if ( lu->options.msglvl > 2 ) {
    488:     PetscFPrintf(PETSC_COMM_SELF,lu->options.msgFile, "\n\n factor matrix after post-processing");
    489:     FrontMtx_writeForHumanEye(lu->frontmtx, lu->options.msgFile);
    490:     fflush(lu->options.msgFile);
    491:   }

    493:   lu->flg = SAME_NONZERO_PATTERN;
    494:   lu->CleanUpSpooles = PETSC_TRUE;
    495:   return(0);
    496: }

    501: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqAIJSpooles(Mat A,MatType type,MatReuse reuse,Mat *newmat) {
    502:   /* This routine is only called to convert a MATSEQAIJ matrix */
    503:   /* to a MATSEQAIJSPOOLES matrix, so we will ignore 'MatType type'. */
    505:   Mat            B=*newmat;
    506:   Mat_Spooles    *lu;

    509:   if (reuse == MAT_INITIAL_MATRIX) {
    510:     /* This routine is inherited, so we know the type is correct. */
    511:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    512:   }
    513:   PetscNew(Mat_Spooles,&lu);
    514:   B->spptr = (void*)lu;

    516:   lu->basetype                   = MATSEQAIJ;
    517:   lu->useQR                      = PETSC_FALSE;
    518:   lu->CleanUpSpooles             = PETSC_FALSE;
    519:   lu->MatDuplicate               = A->ops->duplicate;
    520:   lu->MatCholeskyFactorSymbolic  = A->ops->choleskyfactorsymbolic;
    521:   lu->MatLUFactorSymbolic        = A->ops->lufactorsymbolic;
    522:   lu->MatView                    = A->ops->view;
    523:   lu->MatAssemblyEnd             = A->ops->assemblyend;
    524:   lu->MatDestroy                 = A->ops->destroy;
    525:   B->ops->duplicate              = MatDuplicate_Spooles;
    526:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqAIJSpooles;
    527:   B->ops->lufactorsymbolic       = MatLUFactorSymbolic_SeqAIJSpooles;
    528:   B->ops->view                   = MatView_SeqAIJSpooles;
    529:   B->ops->assemblyend            = MatAssemblyEnd_SeqAIJSpooles;
    530:   B->ops->destroy                = MatDestroy_SeqAIJSpooles;

    532:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaijspooles_seqaij_C",
    533:                                            "MatConvert_Spooles_Base",MatConvert_Spooles_Base);
    534:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijspooles_C",
    535:                                            "MatConvert_SeqAIJ_SeqAIJSpooles",MatConvert_SeqAIJ_SeqAIJSpooles);
    536:   /* PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJSPOOLES); */
    537:   PetscObjectChangeTypeName((PetscObject)B,type);
    538:   *newmat = B;
    539:   return(0);
    540: }

    545: PetscErrorCode MatDuplicate_Spooles(Mat A, MatDuplicateOption op, Mat *M) {
    547:   Mat_Spooles    *lu=(Mat_Spooles *)A->spptr;

    550:   (*lu->MatDuplicate)(A,op,M);
    551:   PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Spooles));
    552:   return(0);
    553: }

    555: /*MC
    556:   MATSEQAIJSPOOLES - MATSEQAIJSPOOLES = "seqaijspooles" - A matrix type providing direct solvers (LU or Cholesky) for sequential matrices
    557:   via the external package SPOOLES.

    559:   If SPOOLES is installed (see the manual for
    560:   instructions on how to declare the existence of external packages),
    561:   a matrix type can be constructed which invokes SPOOLES solvers.
    562:   After calling MatCreate(...,A), simply call MatSetType(A,MATSEQAIJSPOOLES).

    564:   This matrix inherits from MATSEQAIJ.  As a result, MatSeqAIJSetPreallocation is
    565:   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from
    566:   the MATSEQAIJ type without data copy.

    568:   Options Database Keys:
    569: + -mat_type seqaijspooles - sets the matrix type to "seqaijspooles" during a call to MatSetFromOptions()
    570: . -mat_spooles_tau <tau> - upper bound on the magnitude of the largest element in L or U
    571: . -mat_spooles_seed <seed> - random number seed used for ordering
    572: . -mat_spooles_msglvl <msglvl> - message output level
    573: . -mat_spooles_ordering <BestOfNDandMS,MMD,MS,ND> - ordering used
    574: . -mat_spooles_maxdomainsize <n> - maximum subgraph size used by Spooles orderings
    575: . -mat_spooles_maxzeros <n> - maximum number of zeros inside a supernode
    576: . -mat_spooles_maxsize <n> - maximum size of a supernode
    577: . -mat_spooles_FrontMtxInfo <true,fase> - print Spooles information about the computed factorization
    578: . -mat_spooles_symmetryflag <0,1,2> - 0: SPOOLES_SYMMETRIC, 1: SPOOLES_HERMITIAN, 2: SPOOLES_NONSYMMETRIC
    579: . -mat_spooles_patchAndGoFlag <0,1,2> - 0: no patch, 1: use PatchAndGo strategy 1, 2: use PatchAndGo strategy 2
    580: . -mat_spooles_toosmall <dt> - drop tolerance for PatchAndGo strategy 1
    581: . -mat_spooles_storeids <bool integer> - if nonzero, stores row and col numbers where patches were applied in an IV object
    582: . -mat_spooles_fudge <delta> - fudge factor for rescaling diagonals with PatchAndGo strategy 2
    583: - -mat_spooles_storevalues <bool integer> - if nonzero and PatchAndGo strategy 2 is used, store change in diagonal value in a DV object

    585:    Level: beginner

    587: .seealso: PCLU
    588: M*/

    593: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJSpooles(Mat A)
    594: {

    598:   /* Change type name before calling MatSetType to force proper construction of SeqAIJ and SeqAIJSpooles types */
    599:   PetscObjectChangeTypeName((PetscObject)A,MATSEQAIJSPOOLES);
    600:   MatSetType(A,MATSEQAIJ);
    601:   MatConvert_SeqAIJ_SeqAIJSpooles(A,MATSEQAIJSPOOLES,MAT_REUSE_MATRIX,&A);
    602:   return(0);
    603: }

    606: /*MC
    607:   MATAIJSPOOLES - MATAIJSPOOLES = "aijspooles" - A matrix type providing direct solvers (LU or Cholesky) for sequential and parellel matrices
    608:   via the external package SPOOLES.

    610:   If SPOOLES is installed (see the manual for
    611:   instructions on how to declare the existence of external packages),
    612:   a matrix type can be constructed which invokes SPOOLES solvers.
    613:   After calling MatCreate(...,A), simply call MatSetType(A,MATAIJSPOOLES).
    614:   This matrix type is supported for double precision real and complex.

    616:   This matrix inherits from MATAIJ.  As a result, MatSeqAIJSetPreallocation and MatMPIAIJSetPreallocation are
    617:   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from
    618:   the MATAIJ type without data copy.

    620:   Options Database Keys:
    621: + -mat_type aijspooles - sets the matrix type to "aijspooles" during a call to MatSetFromOptions()
    622: . -mat_spooles_tau <tau> - upper bound on the magnitude of the largest element in L or U
    623: . -mat_spooles_seed <seed> - random number seed used for ordering
    624: . -mat_spooles_msglvl <msglvl> - message output level
    625: . -mat_spooles_ordering <BestOfNDandMS,MMD,MS,ND> - ordering used
    626: . -mat_spooles_maxdomainsize <n> - maximum subgraph size used by Spooles orderings
    627: . -mat_spooles_maxzeros <n> - maximum number of zeros inside a supernode
    628: . -mat_spooles_maxsize <n> - maximum size of a supernode
    629: . -mat_spooles_FrontMtxInfo <true,fase> - print Spooles information about the computed factorization
    630: . -mat_spooles_symmetryflag <0,1,2> - 0: SPOOLES_SYMMETRIC, 1: SPOOLES_HERMITIAN, 2: SPOOLES_NONSYMMETRIC
    631: . -mat_spooles_patchAndGoFlag <0,1,2> - 0: no patch, 1: use PatchAndGo strategy 1, 2: use PatchAndGo strategy 2
    632: . -mat_spooles_toosmall <dt> - drop tolerance for PatchAndGo strategy 1
    633: . -mat_spooles_storeids <bool integer> - if nonzero, stores row and col numbers where patches were applied in an IV object
    634: . -mat_spooles_fudge <delta> - fudge factor for rescaling diagonals with PatchAndGo strategy 2
    635: - -mat_spooles_storevalues <bool integer> - if nonzero and PatchAndGo strategy 2 is used, store change in diagonal value in a DV object

    637:    Level: beginner

    639: .seealso: PCLU
    640: M*/
    644: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_AIJSpooles(Mat A)
    645: {
    647:   PetscMPIInt    size;

    650:   /* Change type name before calling MatSetType to force proper construction of SeqAIJSpooles or MPIAIJSpooles */
    651:   PetscObjectChangeTypeName((PetscObject)A,MATAIJSPOOLES);
    652:   MPI_Comm_size(A->comm,&size);
    653:   if (size == 1) {
    654:     MatSetType(A,MATSEQAIJ);
    655:     MatConvert_SeqAIJ_SeqAIJSpooles(A,MATSEQAIJSPOOLES,MAT_REUSE_MATRIX,&A);
    656:   } else {
    657:     MatSetType(A,MATMPIAIJ);
    658:     MatConvert_MPIAIJ_MPIAIJSpooles(A,MATMPIAIJSPOOLES,MAT_REUSE_MATRIX,&A);
    659:   }
    660:   return(0);
    661: }