Actual source code: tfs.c

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

  5: #include <petsc/private/pcimpl.h>
  6: #include <../src/mat/impls/aij/mpi/mpiaij.h>
  7: #include <../src/ksp/pc/impls/tfs/tfs.h>

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

 16: PetscErrorCode PCDestroy_TFS(PC pc)
 17: {
 18:   PC_TFS         *tfs = (PC_TFS*)pc->data;

 20:   /* free the XXT datastructures */
 21:   if (tfs->xxt) {
 22:     XXT_free(tfs->xxt);
 23:   }
 24:   if (tfs->xyt) {
 25:     XYT_free(tfs->xyt);
 26:   }
 27:   VecDestroy(&tfs->b);
 28:   VecDestroy(&tfs->xd);
 29:   VecDestroy(&tfs->xo);
 30:   PetscFree(pc->data);
 31:   return 0;
 32: }

 34: static PetscErrorCode PCApply_TFS_XXT(PC pc,Vec x,Vec y)
 35: {
 36:   PC_TFS            *tfs = (PC_TFS*)pc->data;
 37:   PetscScalar       *yy;
 38:   const PetscScalar *xx;

 40:   VecGetArrayRead(x,&xx);
 41:   VecGetArray(y,&yy);
 42:   XXT_solve(tfs->xxt,yy,(PetscScalar*)xx);
 43:   VecRestoreArrayRead(x,&xx);
 44:   VecRestoreArray(y,&yy);
 45:   return 0;
 46: }

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

 54:   VecGetArrayRead(x,&xx);
 55:   VecGetArray(y,&yy);
 56:   XYT_solve(tfs->xyt,yy,(PetscScalar*)xx);
 57:   VecRestoreArrayRead(x,&xx);
 58:   VecRestoreArray(y,&yy);
 59:   return 0;
 60: }

 62: static PetscErrorCode PCTFSLocalMult_TFS(PC pc,PetscScalar *xin,PetscScalar *xout)
 63: {
 64:   PC_TFS         *tfs = (PC_TFS*)pc->data;
 65:   Mat            A    = pc->pmat;
 66:   Mat_MPIAIJ     *a   = (Mat_MPIAIJ*)A->data;

 68:   VecPlaceArray(tfs->b,xout);
 69:   VecPlaceArray(tfs->xd,xin);
 70:   VecPlaceArray(tfs->xo,xin+tfs->nd);
 71:   MatMult(a->A,tfs->xd,tfs->b);
 72:   MatMultAdd(a->B,tfs->xo,tfs->b,tfs->b);
 73:   VecResetArray(tfs->b);
 74:   VecResetArray(tfs->xd);
 75:   VecResetArray(tfs->xo);
 76:   return 0;
 77: }

 79: static PetscErrorCode PCSetUp_TFS(PC pc)
 80: {
 81:   PC_TFS         *tfs = (PC_TFS*)pc->data;
 82:   Mat            A    = pc->pmat;
 83:   Mat_MPIAIJ     *a   = (Mat_MPIAIJ*)A->data;
 84:   PetscInt       *localtoglobal,ncol,i;
 85:   PetscBool      ismpiaij;

 87:   /*
 88:   PetscBool      issymmetric;
 89:   Petsc Real tol = 0.0;
 90:   */

 93:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&ismpiaij);

 96:   /* generate the local to global mapping */
 97:   ncol = a->A->cmap->n + a->B->cmap->n;
 98:   PetscMalloc1(ncol,&localtoglobal);
 99:   for (i=0; i<a->A->cmap->n; i++) localtoglobal[i] = A->cmap->rstart + i + 1;
100:   for (i=0; i<a->B->cmap->n; i++) localtoglobal[i+a->A->cmap->n] = a->garray[i] + 1;

102:   /* generate the vectors needed for the local solves */
103:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->A->rmap->n,NULL,&tfs->b);
104:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->A->cmap->n,NULL,&tfs->xd);
105:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->B->cmap->n,NULL,&tfs->xo);
106:   tfs->nd = a->A->cmap->n;

108:   /*  MatIsSymmetric(A,tol,&issymmetric); */
109:   /*  if (issymmetric) { */
110:   PetscBarrier((PetscObject)pc);
111:   if (A->symmetric) {
112:     tfs->xxt       = XXT_new();
113:     XXT_factor(tfs->xxt,localtoglobal,A->rmap->n,ncol,(PetscErrorCode (*)(void*,PetscScalar*,PetscScalar*))PCTFSLocalMult_TFS,pc);
114:     pc->ops->apply = PCApply_TFS_XXT;
115:   } else {
116:     tfs->xyt       = XYT_new();
117:     XYT_factor(tfs->xyt,localtoglobal,A->rmap->n,ncol,(PetscErrorCode (*)(void*,PetscScalar*,PetscScalar*))PCTFSLocalMult_TFS,pc);
118:     pc->ops->apply = PCApply_TFS_XYT;
119:   }

121:   PetscFree(localtoglobal);
122:   return 0;
123: }

125: static PetscErrorCode PCSetFromOptions_TFS(PetscOptionItems *PetscOptionsObject,PC pc)
126: {
127:   return 0;
128: }
129: static PetscErrorCode PCView_TFS(PC pc,PetscViewer viewer)
130: {
131:   return 0;
132: }

134: /*MC
135:      PCTFS - A parallel direct solver intended for problems with very few unknowns (like the
136:          coarse grid in multigrid). Performs a Cholesky or LU factorization of a matrix defined by
137:          its local matrix vector product.

139:    Implemented by  Henry M. Tufo III and Paul Fischer originally for Nek5000 and called XXT or XYT

141:    Level: beginner

143:    Notes:
144:     Only implemented for the MPIAIJ matrices

146:     Only works on a solver object that lives on all of PETSC_COMM_WORLD!

148:     Only works for real numbers (is not built if PetscScalar is complex)

150: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC
151: M*/
152: PETSC_EXTERN PetscErrorCode PCCreate_TFS(PC pc)
153: {
154:   PC_TFS         *tfs;
155:   PetscMPIInt    cmp;

157:   MPI_Comm_compare(PETSC_COMM_WORLD,PetscObjectComm((PetscObject)pc),&cmp);
159:   PetscNewLog(pc,&tfs);

161:   tfs->xxt = NULL;
162:   tfs->xyt = NULL;
163:   tfs->b   = NULL;
164:   tfs->xd  = NULL;
165:   tfs->xo  = NULL;
166:   tfs->nd  = 0;

168:   pc->ops->apply               = NULL;
169:   pc->ops->applytranspose      = NULL;
170:   pc->ops->setup               = PCSetUp_TFS;
171:   pc->ops->destroy             = PCDestroy_TFS;
172:   pc->ops->setfromoptions      = PCSetFromOptions_TFS;
173:   pc->ops->view                = PCView_TFS;
174:   pc->ops->applyrichardson     = NULL;
175:   pc->ops->applysymmetricleft  = NULL;
176:   pc->ops->applysymmetricright = NULL;
177:   pc->data                     = (void*)tfs;
178:   return 0;
179: }