Actual source code: pbjacobi.c
petsc-3.7.5 2017-01-01
2: /*
3: Include files needed for the PBJacobi preconditioner:
4: pcimpl.h - private include file intended for use by all preconditioners
5: */
7: #include <petsc/private/matimpl.h>
8: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
10: /*
11: Private context (data structure) for the PBJacobi preconditioner.
12: */
13: typedef struct {
14: const MatScalar *diag;
15: PetscInt bs,mbs;
16: } PC_PBJacobi;
21: static PetscErrorCode PCApply_PBJacobi_1(PC pc,Vec x,Vec y)
22: {
23: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
24: PetscErrorCode ierr;
25: PetscInt i,m = jac->mbs;
26: const MatScalar *diag = jac->diag;
27: const PetscScalar *xx;
28: PetscScalar *yy;
31: VecGetArrayRead(x,&xx);
32: VecGetArray(y,&yy);
33: for (i=0; i<m; i++) yy[i] = diag[i]*xx[i];
34: VecRestoreArrayRead(x,&xx);
35: VecRestoreArray(y,&yy);
36: PetscLogFlops(2.0*m);
37: return(0);
38: }
42: static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y)
43: {
44: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
45: PetscErrorCode ierr;
46: PetscInt i,m = jac->mbs;
47: const MatScalar *diag = jac->diag;
48: PetscScalar x0,x1,*yy;
49: const PetscScalar *xx;
52: VecGetArrayRead(x,&xx);
53: VecGetArray(y,&yy);
54: for (i=0; i<m; i++) {
55: x0 = xx[2*i]; x1 = xx[2*i+1];
56: yy[2*i] = diag[0]*x0 + diag[2]*x1;
57: yy[2*i+1] = diag[1]*x0 + diag[3]*x1;
58: diag += 4;
59: }
60: VecRestoreArrayRead(x,&xx);
61: VecRestoreArray(y,&yy);
62: PetscLogFlops(6.0*m);
63: return(0);
64: }
67: static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y)
68: {
69: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
70: PetscErrorCode ierr;
71: PetscInt i,m = jac->mbs;
72: const MatScalar *diag = jac->diag;
73: PetscScalar x0,x1,x2,*yy;
74: const PetscScalar *xx;
77: VecGetArrayRead(x,&xx);
78: VecGetArray(y,&yy);
79: for (i=0; i<m; i++) {
80: x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2];
82: yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2;
83: yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2;
84: yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2;
85: diag += 9;
86: }
87: VecRestoreArrayRead(x,&xx);
88: VecRestoreArray(y,&yy);
89: PetscLogFlops(15.0*m);
90: return(0);
91: }
94: static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y)
95: {
96: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
97: PetscErrorCode ierr;
98: PetscInt i,m = jac->mbs;
99: const MatScalar *diag = jac->diag;
100: PetscScalar x0,x1,x2,x3,*yy;
101: const PetscScalar *xx;
104: VecGetArrayRead(x,&xx);
105: VecGetArray(y,&yy);
106: for (i=0; i<m; i++) {
107: x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3];
109: yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3;
110: yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3;
111: yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3;
112: yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3;
113: diag += 16;
114: }
115: VecRestoreArrayRead(x,&xx);
116: VecRestoreArray(y,&yy);
117: PetscLogFlops(28.0*m);
118: return(0);
119: }
122: static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y)
123: {
124: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
125: PetscErrorCode ierr;
126: PetscInt i,m = jac->mbs;
127: const MatScalar *diag = jac->diag;
128: PetscScalar x0,x1,x2,x3,x4,*yy;
129: const PetscScalar *xx;
132: VecGetArrayRead(x,&xx);
133: VecGetArray(y,&yy);
134: for (i=0; i<m; i++) {
135: x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4];
137: yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4;
138: yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4;
139: yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4;
140: yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4;
141: yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4;
142: diag += 25;
143: }
144: VecRestoreArrayRead(x,&xx);
145: VecRestoreArray(y,&yy);
146: PetscLogFlops(45.0*m);
147: return(0);
148: }
151: static PetscErrorCode PCApply_PBJacobi_6(PC pc,Vec x,Vec y)
152: {
153: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
154: PetscErrorCode ierr;
155: PetscInt i,m = jac->mbs;
156: const MatScalar *diag = jac->diag;
157: PetscScalar x0,x1,x2,x3,x4,x5,*yy;
158: const PetscScalar *xx;
161: VecGetArrayRead(x,&xx);
162: VecGetArray(y,&yy);
163: for (i=0; i<m; i++) {
164: x0 = xx[6*i]; x1 = xx[6*i+1]; x2 = xx[6*i+2]; x3 = xx[6*i+3]; x4 = xx[6*i+4]; x5 = xx[6*i+5];
166: yy[6*i] = diag[0]*x0 + diag[6]*x1 + diag[12]*x2 + diag[18]*x3 + diag[24]*x4 + diag[30]*x5;
167: yy[6*i+1] = diag[1]*x0 + diag[7]*x1 + diag[13]*x2 + diag[19]*x3 + diag[25]*x4 + diag[31]*x5;
168: yy[6*i+2] = diag[2]*x0 + diag[8]*x1 + diag[14]*x2 + diag[20]*x3 + diag[26]*x4 + diag[32]*x5;
169: yy[6*i+3] = diag[3]*x0 + diag[9]*x1 + diag[15]*x2 + diag[21]*x3 + diag[27]*x4 + diag[33]*x5;
170: yy[6*i+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2 + diag[22]*x3 + diag[28]*x4 + diag[34]*x5;
171: yy[6*i+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2 + diag[23]*x3 + diag[29]*x4 + diag[35]*x5;
172: diag += 36;
173: }
174: VecRestoreArrayRead(x,&xx);
175: VecRestoreArray(y,&yy);
176: PetscLogFlops(66.0*m);
177: return(0);
178: }
181: static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y)
182: {
183: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
184: PetscErrorCode ierr;
185: PetscInt i,m = jac->mbs;
186: const MatScalar *diag = jac->diag;
187: PetscScalar x0,x1,x2,x3,x4,x5,x6,*yy;
188: const PetscScalar *xx;
191: VecGetArrayRead(x,&xx);
192: VecGetArray(y,&yy);
193: for (i=0; i<m; i++) {
194: x0 = xx[7*i]; x1 = xx[7*i+1]; x2 = xx[7*i+2]; x3 = xx[7*i+3]; x4 = xx[7*i+4]; x5 = xx[7*i+5]; x6 = xx[7*i+6];
196: yy[7*i] = diag[0]*x0 + diag[7]*x1 + diag[14]*x2 + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
197: yy[7*i+1] = diag[1]*x0 + diag[8]*x1 + diag[15]*x2 + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6;
198: yy[7*i+2] = diag[2]*x0 + diag[9]*x1 + diag[16]*x2 + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6;
199: yy[7*i+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2 + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6;
200: yy[7*i+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2 + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6;
201: yy[7*i+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2 + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6;
202: yy[7*i+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2 + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6;
203: diag += 49;
204: }
205: VecRestoreArrayRead(x,&xx);
206: VecRestoreArray(y,&yy);
207: PetscLogFlops(80.0*m);
208: return(0);
209: }
210: /* -------------------------------------------------------------------------- */
213: static PetscErrorCode PCSetUp_PBJacobi(PC pc)
214: {
215: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
217: Mat A = pc->pmat;
220: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage");
222: MatInvertBlockDiagonal(A,&jac->diag);
223: if (A->errortype) pc->failedreason = (PCFailedReason)A->errortype;
224:
225: MatGetBlockSize(A,&jac->bs);
226: jac->mbs = A->rmap->n/jac->bs;
227: switch (jac->bs) {
228: case 1:
229: pc->ops->apply = PCApply_PBJacobi_1;
230: break;
231: case 2:
232: pc->ops->apply = PCApply_PBJacobi_2;
233: break;
234: case 3:
235: pc->ops->apply = PCApply_PBJacobi_3;
236: break;
237: case 4:
238: pc->ops->apply = PCApply_PBJacobi_4;
239: break;
240: case 5:
241: pc->ops->apply = PCApply_PBJacobi_5;
242: break;
243: case 6:
244: pc->ops->apply = PCApply_PBJacobi_6;
245: break;
246: case 7:
247: pc->ops->apply = PCApply_PBJacobi_7;
248: break;
249: default:
250: SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"not supported for block size %D",jac->bs);
251: }
252: return(0);
253: }
254: /* -------------------------------------------------------------------------- */
257: static PetscErrorCode PCDestroy_PBJacobi(PC pc)
258: {
262: /*
263: Free the private data structure that was hanging off the PC
264: */
265: PetscFree(pc->data);
266: return(0);
267: }
271: static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer)
272: {
274: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
275: PetscBool iascii;
278: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
279: if (iascii) {
280: PetscViewerASCIIPrintf(viewer," point-block Jacobi: block size %D\n",jac->bs);
281: }
282: return(0);
283: }
285: /* -------------------------------------------------------------------------- */
286: /*MC
287: PCPBJACOBI - Point block Jacobi preconditioner
290: Notes: See PCJACOBI for point Jacobi preconditioning
292: This works for AIJ and BAIJ matrices and uses the blocksize provided to the matrix
294: Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot
295: is detected a PETSc error is generated.
297: Developer Notes: This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow
298: the factorization to continue even after a zero pivot is found resulting in a Nan and hence
299: terminating KSP with a KSP_DIVERGED_NANORIF allowing
300: a nonlinear solver/ODE integrator to recover without stopping the program as currently happens.
302: Developer Note: Perhaps should provide an option that allows generation of a valid preconditioner
303: even if a block is singular as the PCJACOBI does.
305: Level: beginner
307: Concepts: point block Jacobi
310: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI
312: M*/
316: PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc)
317: {
318: PC_PBJacobi *jac;
322: /*
323: Creates the private data structure for this preconditioner and
324: attach it to the PC object.
325: */
326: PetscNewLog(pc,&jac);
327: pc->data = (void*)jac;
329: /*
330: Initialize the pointers to vectors to ZERO; these will be used to store
331: diagonal entries of the matrix for fast preconditioner application.
332: */
333: jac->diag = 0;
335: /*
336: Set the pointers for the functions that are provided above.
337: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
338: are called, they will automatically call these functions. Note we
339: choose not to provide a couple of these functions since they are
340: not needed.
341: */
342: pc->ops->apply = 0; /*set depending on the block size */
343: pc->ops->applytranspose = 0;
344: pc->ops->setup = PCSetUp_PBJacobi;
345: pc->ops->destroy = PCDestroy_PBJacobi;
346: pc->ops->setfromoptions = 0;
347: pc->ops->view = PCView_PBJacobi;
348: pc->ops->applyrichardson = 0;
349: pc->ops->applysymmetricleft = 0;
350: pc->ops->applysymmetricright = 0;
351: return(0);
352: }