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: }