1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 #include <../src/mat/impls/aij/seq/aij.h> 3 4 /* 5 Private context (data structure) for the CP preconditioner. 6 */ 7 typedef struct { 8 PetscInt n, m; 9 Vec work; 10 PetscScalar *d; /* sum of squares of each column */ 11 PetscScalar *a; /* non-zeros by column */ 12 PetscInt *i, *j; /* offsets of nonzeros by column, non-zero indices by column */ 13 } PC_CP; 14 15 static PetscErrorCode PCSetUp_CP(PC pc) 16 { 17 PC_CP *cp = (PC_CP *)pc->data; 18 PetscInt i, j, *colcnt; 19 PetscBool flg; 20 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)pc->pmat->data; 21 22 PetscFunctionBegin; 23 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg)); 24 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles SeqAIJ matrices"); 25 26 PetscCall(MatGetLocalSize(pc->pmat, &cp->m, &cp->n)); 27 PetscCheck(cp->m == cp->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Currently only for square matrices"); 28 29 if (!cp->work) PetscCall(MatCreateVecs(pc->pmat, &cp->work, NULL)); 30 if (!cp->d) PetscCall(PetscMalloc1(cp->n, &cp->d)); 31 if (cp->a && pc->flag != SAME_NONZERO_PATTERN) { 32 PetscCall(PetscFree3(cp->a, cp->i, cp->j)); 33 cp->a = NULL; 34 } 35 36 /* convert to column format */ 37 if (!cp->a) PetscCall(PetscMalloc3(aij->nz, &cp->a, cp->n + 1, &cp->i, aij->nz, &cp->j)); 38 PetscCall(PetscCalloc1(cp->n, &colcnt)); 39 40 for (i = 0; i < aij->nz; i++) colcnt[aij->j[i]]++; 41 cp->i[0] = 0; 42 for (i = 0; i < cp->n; i++) cp->i[i + 1] = cp->i[i] + colcnt[i]; 43 PetscCall(PetscArrayzero(colcnt, cp->n)); 44 for (i = 0; i < cp->m; i++) { /* over rows */ 45 for (j = aij->i[i]; j < aij->i[i + 1]; j++) { /* over columns in row */ 46 cp->j[cp->i[aij->j[j]] + colcnt[aij->j[j]]] = i; 47 cp->a[cp->i[aij->j[j]] + colcnt[aij->j[j]]++] = aij->a[j]; 48 } 49 } 50 PetscCall(PetscFree(colcnt)); 51 52 /* compute sum of squares of each column d[] */ 53 for (i = 0; i < cp->n; i++) { /* over columns */ 54 cp->d[i] = 0.; 55 for (j = cp->i[i]; j < cp->i[i + 1]; j++) cp->d[i] += cp->a[j] * cp->a[j]; /* over rows in column */ 56 cp->d[i] = 1.0 / cp->d[i]; 57 } 58 PetscFunctionReturn(PETSC_SUCCESS); 59 } 60 61 static PetscErrorCode PCApply_CP(PC pc, Vec bb, Vec xx) 62 { 63 PC_CP *cp = (PC_CP *)pc->data; 64 PetscScalar *b, *x, xt; 65 PetscInt i, j; 66 67 PetscFunctionBegin; 68 PetscCall(VecCopy(bb, cp->work)); 69 PetscCall(VecGetArray(cp->work, &b)); 70 PetscCall(VecGetArray(xx, &x)); 71 72 for (i = 0; i < cp->n; i++) { /* over columns */ 73 xt = 0.; 74 for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */ 75 xt *= cp->d[i]; 76 x[i] = xt; 77 for (j = cp->i[i]; j < cp->i[i + 1]; j++) b[cp->j[j]] -= xt * cp->a[j]; /* over rows in column updating b*/ 78 } 79 for (i = cp->n - 1; i > -1; i--) { /* over columns */ 80 xt = 0.; 81 for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */ 82 xt *= cp->d[i]; 83 x[i] = xt; 84 for (j = cp->i[i]; j < cp->i[i + 1]; j++) b[cp->j[j]] -= xt * cp->a[j]; /* over rows in column updating b*/ 85 } 86 87 PetscCall(VecRestoreArray(cp->work, &b)); 88 PetscCall(VecRestoreArray(xx, &x)); 89 PetscFunctionReturn(PETSC_SUCCESS); 90 } 91 92 static PetscErrorCode PCReset_CP(PC pc) 93 { 94 PC_CP *cp = (PC_CP *)pc->data; 95 96 PetscFunctionBegin; 97 PetscCall(PetscFree(cp->d)); 98 PetscCall(VecDestroy(&cp->work)); 99 PetscCall(PetscFree3(cp->a, cp->i, cp->j)); 100 PetscFunctionReturn(PETSC_SUCCESS); 101 } 102 103 static PetscErrorCode PCDestroy_CP(PC pc) 104 { 105 PC_CP *cp = (PC_CP *)pc->data; 106 107 PetscFunctionBegin; 108 PetscCall(PCReset_CP(pc)); 109 PetscCall(PetscFree(cp->d)); 110 PetscCall(PetscFree3(cp->a, cp->i, cp->j)); 111 PetscCall(PetscFree(pc->data)); 112 PetscFunctionReturn(PETSC_SUCCESS); 113 } 114 115 static PetscErrorCode PCSetFromOptions_CP(PC pc, PetscOptionItems *PetscOptionsObject) 116 { 117 PetscFunctionBegin; 118 PetscFunctionReturn(PETSC_SUCCESS); 119 } 120 121 /*MC 122 PCCP - a "column-projection" preconditioner 123 124 This is a terrible preconditioner and is not recommended, ever! 125 126 Loops over the entries of x computing dx_i (e_i is the unit vector in the ith direction) to 127 .vb 128 129 min || b - A(x + dx_i e_i ||_2 130 dx_i 131 132 That is, it changes a single entry of x to minimize the new residual norm. 133 Let A_i represent the ith column of A, then the minimization can be written as 134 135 min || r - (dx_i) A e_i ||_2 136 dx_i 137 or min || r - (dx_i) A_i ||_2 138 dx_i 139 140 take the derivative with respect to dx_i to obtain 141 dx_i = (A_i^T A_i)^(-1) A_i^T r 142 143 This algorithm can be thought of as Gauss-Seidel on the normal equations 144 .ve 145 146 Notes: 147 This procedure can also be done with block columns or any groups of columns 148 but this is not coded. 149 150 These "projections" can be done simultaneously for all columns (similar to Jacobi) 151 or sequentially (similar to Gauss-Seidel/SOR). This is only coded for SOR type. 152 153 This is related to, but not the same as "row projection" methods. 154 155 This is currently coded only for `MATSEQAIJ` matrices 156 157 Level: intermediate 158 159 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCJACOBI`, `PCSOR` 160 M*/ 161 162 PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc) 163 { 164 PC_CP *cp; 165 166 PetscFunctionBegin; 167 PetscCall(PetscNew(&cp)); 168 pc->data = (void *)cp; 169 170 pc->ops->apply = PCApply_CP; 171 pc->ops->applytranspose = PCApply_CP; 172 pc->ops->setup = PCSetUp_CP; 173 pc->ops->reset = PCReset_CP; 174 pc->ops->destroy = PCDestroy_CP; 175 pc->ops->setfromoptions = PCSetFromOptions_CP; 176 pc->ops->view = NULL; 177 pc->ops->applyrichardson = NULL; 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180