Actual source code: vpbjacobi.c
petsc-3.10.5 2019-03-28
2: /*
3: Include files needed for the variable size block PBJacobi preconditioner:
4: pcimpl.h - private include file intended for use by all preconditioners
5: */
7: #include <petsc/private/pcimpl.h>
9: /*
10: Private context (data structure) for the VPBJacobi preconditioner.
11: */
12: typedef struct {
13: MatScalar *diag;
14: } PC_VPBJacobi;
17: static PetscErrorCode PCApply_VPBJacobi(PC pc,Vec x,Vec y)
18: {
19: PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data;
20: PetscErrorCode ierr;
21: PetscInt i,ncnt = 0;
22: const MatScalar *diag = jac->diag;
23: PetscInt ib,jb,bs;
24: const PetscScalar *xx;
25: PetscScalar *yy,x0,x1,x2,x3,x4,x5,x6;
26: PetscInt nblocks;
27: const PetscInt *bsizes;
30: MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes);
31: VecGetArrayRead(x,&xx);
32: VecGetArray(y,&yy);
33: for (i=0; i<nblocks; i++) {
34: bs = bsizes[i];
35: switch (bs) {
36: case 1:
37: yy[ncnt] = *diag*xx[ncnt];
38: break;
39: case 2:
40: x0 = xx[ncnt]; x1 = xx[ncnt+1];
41: yy[ncnt] = diag[0]*x0 + diag[2]*x1;
42: yy[ncnt+1] = diag[1]*x0 + diag[3]*x1;
43: break;
44: case 3:
45: x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2];
46: yy[ncnt] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2;
47: yy[ncnt+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2;
48: yy[ncnt+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2;
49: break;
50: case 4:
51: x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3];
52: yy[ncnt] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3;
53: yy[ncnt+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3;
54: yy[ncnt+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3;
55: yy[ncnt+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3;
56: break;
57: case 5:
58: x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4];
59: yy[ncnt] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4;
60: yy[ncnt+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4;
61: yy[ncnt+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4;
62: yy[ncnt+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4;
63: yy[ncnt+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4;
64: break;
65: case 6:
66: x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4]; x5 = xx[ncnt+5];
67: yy[ncnt] = diag[0]*x0 + diag[6]*x1 + diag[12]*x2 + diag[18]*x3 + diag[24]*x4 + diag[30]*x5;
68: yy[ncnt+1] = diag[1]*x0 + diag[7]*x1 + diag[13]*x2 + diag[19]*x3 + diag[25]*x4 + diag[31]*x5;
69: yy[ncnt+2] = diag[2]*x0 + diag[8]*x1 + diag[14]*x2 + diag[20]*x3 + diag[26]*x4 + diag[32]*x5;
70: yy[ncnt+3] = diag[3]*x0 + diag[9]*x1 + diag[15]*x2 + diag[21]*x3 + diag[27]*x4 + diag[33]*x5;
71: yy[ncnt+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2 + diag[22]*x3 + diag[28]*x4 + diag[34]*x5;
72: yy[ncnt+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2 + diag[23]*x3 + diag[29]*x4 + diag[35]*x5;
73: break;
74: case 7:
75: x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4]; x5 = xx[ncnt+5]; x6 = xx[ncnt+6];
76: yy[ncnt] = diag[0]*x0 + diag[7]*x1 + diag[14]*x2 + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
77: yy[ncnt+1] = diag[1]*x0 + diag[8]*x1 + diag[15]*x2 + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6;
78: yy[ncnt+2] = diag[2]*x0 + diag[9]*x1 + diag[16]*x2 + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6;
79: yy[ncnt+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2 + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6;
80: yy[ncnt+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2 + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6;
81: yy[ncnt+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2 + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6;
82: yy[ncnt+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2 + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6;
83: break;
84: default:
85: for (ib=0; ib<bs; ib++){
86: PetscScalar rowsum = 0;
87: for (jb=0; jb<bs; jb++){
88: rowsum += diag[ib+jb*bs] * xx[ncnt+jb];
89: }
90: yy[ncnt+ib] = rowsum;
91: }
92: }
93: ncnt += bsizes[i];
94: diag += bsizes[i]*bsizes[i];
95: }
96: VecRestoreArrayRead(x,&xx);
97: VecRestoreArray(y,&yy);
98: return(0);
99: }
103: /* -------------------------------------------------------------------------- */
104: static PetscErrorCode PCSetUp_VPBJacobi(PC pc)
105: {
106: PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data;
108: Mat A = pc->pmat;
109: MatFactorError err;
110: PetscInt i,nsize = 0,nlocal;
111: PetscInt nblocks;
112: const PetscInt *bsizes;
115: MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes);
116: MatGetLocalSize(pc->pmat,&nlocal,NULL);
117: if (nlocal && !nblocks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatSetVariableBlockSizes() before using PCVPBJACOBI");
118: if (!jac->diag) {
119: for (i=0; i<nblocks; i++) nsize += bsizes[i]*bsizes[i];
120: PetscMalloc1(nsize,&jac->diag);
121: }
122: MatInvertVariableBlockDiagonal(A,nblocks,bsizes,jac->diag);
123: MatFactorGetError(A,&err);
124: if (err) pc->failedreason = (PCFailedReason)err;
125: pc->ops->apply = PCApply_VPBJacobi;
126: return(0);
127: }
128: /* -------------------------------------------------------------------------- */
129: static PetscErrorCode PCDestroy_VPBJacobi(PC pc)
130: {
131: PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data;
135: /*
136: Free the private data structure that was hanging off the PC
137: */
138: PetscFree(jac->diag);
139: PetscFree(pc->data);
140: return(0);
141: }
143: /* -------------------------------------------------------------------------- */
144: /*MC
145: PCVPBJACOBI - Variable size point block Jacobi preconditioner
148: Notes:
149: See PCJACOBI for point Jacobi preconditioning, PCPBJACOBI for fixed point block size, and PCBJACOBI for large size blocks
151: This works for AIJ matrices
153: Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot
154: is detected a PETSc error is generated.
156: One must call MatSetVariableBlockSizes() to use this preconditioner
157: Developer Notes:
158: This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow
159: the factorization to continue even after a zero pivot is found resulting in a Nan and hence
160: terminating KSP with a KSP_DIVERGED_NANORIF allowing
161: a nonlinear solver/ODE integrator to recover without stopping the program as currently happens.
163: Perhaps should provide an option that allows generation of a valid preconditioner
164: even if a block is singular as the PCJACOBI does.
166: Level: beginner
168: Concepts: variable point block Jacobi
170: .seealso: MatSetVariableBlockSizes(), PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI
172: M*/
174: PETSC_EXTERN PetscErrorCode PCCreate_VPBJacobi(PC pc)
175: {
176: PC_VPBJacobi *jac;
180: /*
181: Creates the private data structure for this preconditioner and
182: attach it to the PC object.
183: */
184: PetscNewLog(pc,&jac);
185: pc->data = (void*)jac;
187: /*
188: Initialize the pointers to vectors to ZERO; these will be used to store
189: diagonal entries of the matrix for fast preconditioner application.
190: */
191: jac->diag = NULL;
193: /*
194: Set the pointers for the functions that are provided above.
195: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
196: are called, they will automatically call these functions. Note we
197: choose not to provide a couple of these functions since they are
198: not needed.
199: */
200: pc->ops->apply = PCApply_VPBJacobi;
201: pc->ops->applytranspose = 0;
202: pc->ops->setup = PCSetUp_VPBJacobi;
203: pc->ops->destroy = PCDestroy_VPBJacobi;
204: pc->ops->setfromoptions = 0;
205: pc->ops->applyrichardson = 0;
206: pc->ops->applysymmetricleft = 0;
207: pc->ops->applysymmetricright = 0;
208: return(0);
209: }