Actual source code: tfs.c

      1: #define PETSCKSP_DLL
      2: /*
      3:         Provides an interface to the Tufo-Fischer parallel direct solver
      4: */

     6:  #include src/ksp/pc/pcimpl.h
     7:  #include src/mat/impls/aij/mpi/mpiaij.h
     8:  #include src/ksp/pc/impls/tfs/tfs.h

     10: typedef struct {
     11:   xxt_ADT  xxt;
     12:   xyt_ADT  xyt;
     13:   Vec      b,xd,xo;
     14:   PetscInt nd;
     15: } PC_TFS;

     19: PetscErrorCode PCDestroy_TFS(PC pc)
     20: {
     21:   PC_TFS *tfs = (PC_TFS*)pc->data;

     25:   /* free the XXT datastructures */
     26:   if (tfs->xxt) {
     27:     XXT_free(tfs->xxt);
     28:   }
     29:   if (tfs->xyt) {
     30:     XYT_free(tfs->xyt);
     31:   }
     32:   if (tfs->b) {
     33:   VecDestroy(tfs->b);
     34:   }
     35:   if (tfs->xd) {
     36:   VecDestroy(tfs->xd);
     37:   }
     38:   if (tfs->xo) {
     39:   VecDestroy(tfs->xo);
     40:   }
     41:   PetscFree(tfs);
     42:   return(0);
     43: }

     47: static PetscErrorCode PCApply_TFS_XXT(PC pc,Vec x,Vec y)
     48: {
     49:   PC_TFS *tfs = (PC_TFS*)pc->data;
     50:   PetscScalar    *xx,*yy;

     54:   VecGetArray(x,&xx);
     55:   VecGetArray(y,&yy);
     56:   XXT_solve(tfs->xxt,yy,xx);
     57:   VecRestoreArray(x,&xx);
     58:   VecRestoreArray(y,&yy);
     59:   return(0);
     60: }

     64: static PetscErrorCode PCApply_TFS_XYT(PC pc,Vec x,Vec y)
     65: {
     66:   PC_TFS *tfs = (PC_TFS*)pc->data;
     67:   PetscScalar    *xx,*yy;

     71:   VecGetArray(x,&xx);
     72:   VecGetArray(y,&yy);
     73:   XYT_solve(tfs->xyt,yy,xx);
     74:   VecRestoreArray(x,&xx);
     75:   VecRestoreArray(y,&yy);
     76:   return(0);
     77: }

     81: static PetscErrorCode LocalMult_TFS(PC pc,PetscScalar *xin,PetscScalar *xout)
     82: {
     83:   PC_TFS        *tfs = (PC_TFS*)pc->data;
     84:   Mat           A = pc->pmat;
     85:   Mat_MPIAIJ    *a = (Mat_MPIAIJ*)A->data;
     87:
     89:   VecPlaceArray(tfs->b,xout);
     90:   VecPlaceArray(tfs->xd,xin);
     91:   VecPlaceArray(tfs->xo,xin+tfs->nd);
     92:   MatMult(a->A,tfs->xd,tfs->b);
     93:   MatMultAdd(a->B,tfs->xo,tfs->b,tfs->b);
     94:   return(0);
     95: }

     99: static PetscErrorCode PCSetUp_TFS(PC pc)
    100: {
    101:   PC_TFS        *tfs = (PC_TFS*)pc->data;
    102:   Mat            A = pc->pmat;
    103:   Mat_MPIAIJ     *a = (Mat_MPIAIJ*)A->data;
    105:   PetscInt      *localtoglobal,ncol,i;
    106:   PetscTruth     ismpiaij;

    108:   /*
    109:   PetscTruth     issymmetric;
    110:   Petsc Real tol = 0.0;
    111:   */

    114:   if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
    115:   PetscTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&ismpiaij);
    116:   if (!ismpiaij) {
    117:     SETERRQ(PETSC_ERR_SUP,"Currently only supports MPIAIJ matrices");
    118:   }

    120:   /* generate the local to global mapping */
    121:   ncol = a->A->n + a->B->n;
    122:   PetscMalloc((ncol)*sizeof(int),&localtoglobal);
    123:   for (i=0; i<a->A->n; i++) {
    124:     localtoglobal[i] = a->cstart + i + 1;
    125:   }
    126:   for (i=0; i<a->B->n; i++) {
    127:     localtoglobal[i+a->A->n] = a->garray[i] + 1;
    128:   }
    129:   /* generate the vectors needed for the local solves */
    130:   VecCreateSeqWithArray(PETSC_COMM_SELF,a->A->m,PETSC_NULL,&tfs->b);
    131:   VecCreateSeqWithArray(PETSC_COMM_SELF,a->A->n,PETSC_NULL,&tfs->xd);
    132:   VecCreateSeqWithArray(PETSC_COMM_SELF,a->B->n,PETSC_NULL,&tfs->xo);
    133:   tfs->nd = a->A->n;

    136:   /*   MatIsSymmetric(A,tol,&issymmetric); */
    137:   /*  if (issymmetric) { */
    138:   PetscBarrier((PetscObject)pc);
    139:   if (A->symmetric) {
    140:     tfs->xxt       = XXT_new();
    141:     XXT_factor(tfs->xxt,localtoglobal,A->m,ncol,(void*)LocalMult_TFS,pc);
    142:     pc->ops->apply = PCApply_TFS_XXT;
    143:   } else {
    144:     tfs->xyt       = XYT_new();
    145:     XYT_factor(tfs->xyt,localtoglobal,A->m,ncol,(void*)LocalMult_TFS,pc);
    146:     pc->ops->apply = PCApply_TFS_XYT;
    147:   }

    149:   PetscFree(localtoglobal);
    150:   return(0);
    151: }

    155: static PetscErrorCode PCSetFromOptions_TFS(PC pc)
    156: {
    158:   return(0);
    159: }
    162: static PetscErrorCode PCView_TFS(PC pc,PetscViewer viewer)
    163: {
    165:   return(0);
    166: }

    171: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_TFS(PC pc)
    172: {
    174:   PC_TFS         *tfs;

    177:   PetscNew(PC_TFS,&tfs);
    178:   PetscLogObjectMemory(pc,sizeof(PC_TFS));

    180:   tfs->xxt = 0;
    181:   tfs->xyt = 0;
    182:   tfs->b   = 0;
    183:   tfs->xd  = 0;
    184:   tfs->xo  = 0;
    185:   tfs->nd  = 0;

    187:   pc->ops->apply               = 0;
    188:   pc->ops->applytranspose      = 0;
    189:   pc->ops->setup               = PCSetUp_TFS;
    190:   pc->ops->destroy             = PCDestroy_TFS;
    191:   pc->ops->setfromoptions      = PCSetFromOptions_TFS;
    192:   pc->ops->view                = PCView_TFS;
    193:   pc->ops->applyrichardson     = 0;
    194:   pc->ops->applysymmetricleft  = 0;
    195:   pc->ops->applysymmetricright = 0;
    196:   pc->data                     = (void*)tfs;
    197:   return(0);
    198: }