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

  5: #include <petsc-private/pcimpl.h>   /*I "petscpc.h" I*/
  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;

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

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

 40: static PetscErrorCode PCApply_TFS_XXT(PC pc,Vec x,Vec y)
 41: {
 42:   PC_TFS *tfs = (PC_TFS*)pc->data;
 43:   PetscScalar    *xx,*yy;

 47:   VecGetArray(x,&xx);
 48:   VecGetArray(y,&yy);
 49:   XXT_solve(tfs->xxt,yy,xx);
 50:   VecRestoreArray(x,&xx);
 51:   VecRestoreArray(y,&yy);
 52:   return(0);
 53: }

 57: static PetscErrorCode PCApply_TFS_XYT(PC pc,Vec x,Vec y)
 58: {
 59:   PC_TFS *tfs = (PC_TFS*)pc->data;
 60:   PetscScalar    *xx,*yy;

 64:   VecGetArray(x,&xx);
 65:   VecGetArray(y,&yy);
 66:   XYT_solve(tfs->xyt,yy,xx);
 67:   VecRestoreArray(x,&xx);
 68:   VecRestoreArray(y,&yy);
 69:   return(0);
 70: }

 74: static PetscErrorCode PCTFSLocalMult_TFS(PC pc,PetscScalar *xin,PetscScalar *xout)
 75: {
 76:   PC_TFS        *tfs = (PC_TFS*)pc->data;
 77:   Mat           A = pc->pmat;
 78:   Mat_MPIAIJ    *a = (Mat_MPIAIJ*)A->data;
 80: 
 82:   VecPlaceArray(tfs->b,xout);
 83:   VecPlaceArray(tfs->xd,xin);
 84:   VecPlaceArray(tfs->xo,xin+tfs->nd);
 85:   MatMult(a->A,tfs->xd,tfs->b);
 86:   MatMultAdd(a->B,tfs->xo,tfs->b,tfs->b);
 87:   VecResetArray(tfs->b);
 88:   VecResetArray(tfs->xd);
 89:   VecResetArray(tfs->xo);
 90:   return(0);
 91: }

 95: static PetscErrorCode PCSetUp_TFS(PC pc)
 96: {
 97:   PC_TFS        *tfs = (PC_TFS*)pc->data;
 98:   Mat            A = pc->pmat;
 99:   Mat_MPIAIJ     *a = (Mat_MPIAIJ*)A->data;
101:   PetscInt      *localtoglobal,ncol,i;
102:   PetscBool      ismpiaij;

104:   /*
105:   PetscBool      issymmetric;
106:   Petsc Real tol = 0.0;
107:   */

110:   if (A->cmap->N != A->rmap->N) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_SIZ,"matrix must be square");
111:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&ismpiaij);
112:   if (!ismpiaij) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Currently only supports MPIAIJ matrices");

114:   /* generate the local to global mapping */
115:   ncol = a->A->cmap->n + a->B->cmap->n;
116:   PetscMalloc((ncol)*sizeof(PetscInt),&localtoglobal);
117:   for (i=0; i<a->A->cmap->n; i++) {
118:     localtoglobal[i] = A->cmap->rstart + i + 1;
119:   }
120:   for (i=0; i<a->B->cmap->n; i++) {
121:     localtoglobal[i+a->A->cmap->n] = a->garray[i] + 1;
122:   }
123:   /* generate the vectors needed for the local solves */
124:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->A->rmap->n,PETSC_NULL,&tfs->b);
125:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->A->cmap->n,PETSC_NULL,&tfs->xd);
126:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->B->cmap->n,PETSC_NULL,&tfs->xo);
127:   tfs->nd = a->A->cmap->n;


130:   /*   MatIsSymmetric(A,tol,&issymmetric); */
131:   /*  if (issymmetric) { */
132:   PetscBarrier((PetscObject)pc);
133:   if (A->symmetric) {
134:     tfs->xxt       = XXT_new();
135:     XXT_factor(tfs->xxt,localtoglobal,A->rmap->n,ncol,(void*)PCTFSLocalMult_TFS,pc);
136:     pc->ops->apply = PCApply_TFS_XXT;
137:   } else {
138:     tfs->xyt       = XYT_new();
139:     XYT_factor(tfs->xyt,localtoglobal,A->rmap->n,ncol,(void*)PCTFSLocalMult_TFS,pc);
140:     pc->ops->apply = PCApply_TFS_XYT;
141:   }

143:   PetscFree(localtoglobal);
144:   return(0);
145: }

149: static PetscErrorCode PCSetFromOptions_TFS(PC pc)
150: {
152:   return(0);
153: }
156: static PetscErrorCode PCView_TFS(PC pc,PetscViewer viewer)
157: {
159:   return(0);
160: }

162: EXTERN_C_BEGIN
165: /*MC
166:      PCTFS - A parallel direct solver intended for problems with very few unknowns (like the 
167:          coarse grid in multigrid).

169:    Implemented by  Henry M. Tufo III and Paul Fischer

171:    Level: beginner

173:    Notes: Only implemented for the MPIAIJ matrices

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

177: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC
178: M*/
179: PetscErrorCode  PCCreate_TFS(PC pc)
180: {
182:   PC_TFS         *tfs;
183:   PetscMPIInt    cmp;

186:   MPI_Comm_compare(PETSC_COMM_WORLD,((PetscObject)pc)->comm,&cmp);
187:   if (cmp != MPI_IDENT && cmp != MPI_CONGRUENT) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"TFS only works with PETSC_COMM_WORLD objects");
188:   PetscNewLog(pc,PC_TFS,&tfs);

190:   tfs->xxt = 0;
191:   tfs->xyt = 0;
192:   tfs->b   = 0;
193:   tfs->xd  = 0;
194:   tfs->xo  = 0;
195:   tfs->nd  = 0;

197:   pc->ops->apply               = 0;
198:   pc->ops->applytranspose      = 0;
199:   pc->ops->setup               = PCSetUp_TFS;
200:   pc->ops->destroy             = PCDestroy_TFS;
201:   pc->ops->setfromoptions      = PCSetFromOptions_TFS;
202:   pc->ops->view                = PCView_TFS;
203:   pc->ops->applyrichardson     = 0;
204:   pc->ops->applysymmetricleft  = 0;
205:   pc->ops->applysymmetricright = 0;
206:   pc->data                     = (void*)tfs;
207:   return(0);
208: }
209: EXTERN_C_END