Actual source code: cp.c


  2: #include <petsc/private/pcimpl.h>
  3: #include <../src/mat/impls/aij/seq/aij.h>

  5: /*
  6:    Private context (data structure) for the CP preconditioner.
  7: */
  8: typedef struct {
  9:   PetscInt     n, m;
 10:   Vec          work;
 11:   PetscScalar *d;     /* sum of squares of each column */
 12:   PetscScalar *a;     /* non-zeros by column */
 13:   PetscInt    *i, *j; /* offsets of nonzeros by column, non-zero indices by column */
 14: } PC_CP;

 16: static PetscErrorCode PCSetUp_CP(PC pc)
 17: {
 18:   PC_CP      *cp = (PC_CP *)pc->data;
 19:   PetscInt    i, j, *colcnt;
 20:   PetscBool   flg;
 21:   Mat_SeqAIJ *aij = (Mat_SeqAIJ *)pc->pmat->data;

 23:   PetscFunctionBegin;
 24:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg));
 25:   PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only handles SeqAIJ matrices");

 27:   PetscCall(MatGetLocalSize(pc->pmat, &cp->m, &cp->n));
 28:   PetscCheck(cp->m == cp->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Currently only for square matrices");

 30:   if (!cp->work) PetscCall(MatCreateVecs(pc->pmat, &cp->work, NULL));
 31:   if (!cp->d) PetscCall(PetscMalloc1(cp->n, &cp->d));
 32:   if (cp->a && pc->flag != SAME_NONZERO_PATTERN) {
 33:     PetscCall(PetscFree3(cp->a, cp->i, cp->j));
 34:     cp->a = NULL;
 35:   }

 37:   /* convert to column format */
 38:   if (!cp->a) PetscCall(PetscMalloc3(aij->nz, &cp->a, cp->n + 1, &cp->i, aij->nz, &cp->j));
 39:   PetscCall(PetscCalloc1(cp->n, &colcnt));

 41:   for (i = 0; i < aij->nz; i++) colcnt[aij->j[i]]++;
 42:   cp->i[0] = 0;
 43:   for (i = 0; i < cp->n; i++) cp->i[i + 1] = cp->i[i] + colcnt[i];
 44:   PetscCall(PetscArrayzero(colcnt, cp->n));
 45:   for (i = 0; i < cp->m; i++) {                   /* over rows */
 46:     for (j = aij->i[i]; j < aij->i[i + 1]; j++) { /* over columns in row */
 47:       cp->j[cp->i[aij->j[j]] + colcnt[aij->j[j]]]   = i;
 48:       cp->a[cp->i[aij->j[j]] + colcnt[aij->j[j]]++] = aij->a[j];
 49:     }
 50:   }
 51:   PetscCall(PetscFree(colcnt));

 53:   /* compute sum of squares of each column d[] */
 54:   for (i = 0; i < cp->n; i++) { /* over columns */
 55:     cp->d[i] = 0.;
 56:     for (j = cp->i[i]; j < cp->i[i + 1]; j++) cp->d[i] += cp->a[j] * cp->a[j]; /* over rows in column */
 57:     cp->d[i] = 1.0 / cp->d[i];
 58:   }
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

 62: static PetscErrorCode PCApply_CP(PC pc, Vec bb, Vec xx)
 63: {
 64:   PC_CP       *cp = (PC_CP *)pc->data;
 65:   PetscScalar *b, *x, xt;
 66:   PetscInt     i, j;

 68:   PetscFunctionBegin;
 69:   PetscCall(VecCopy(bb, cp->work));
 70:   PetscCall(VecGetArray(cp->work, &b));
 71:   PetscCall(VecGetArray(xx, &x));

 73:   for (i = 0; i < cp->n; i++) { /* over columns */
 74:     xt = 0.;
 75:     for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */
 76:     xt *= cp->d[i];
 77:     x[i] = xt;
 78:     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*/
 79:   }
 80:   for (i = cp->n - 1; i > -1; i--) { /* over columns */
 81:     xt = 0.;
 82:     for (j = cp->i[i]; j < cp->i[i + 1]; j++) xt += cp->a[j] * b[cp->j[j]]; /* over rows in column */
 83:     xt *= cp->d[i];
 84:     x[i] = xt;
 85:     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*/
 86:   }

 88:   PetscCall(VecRestoreArray(cp->work, &b));
 89:   PetscCall(VecRestoreArray(xx, &x));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: static PetscErrorCode PCReset_CP(PC pc)
 94: {
 95:   PC_CP *cp = (PC_CP *)pc->data;

 97:   PetscFunctionBegin;
 98:   PetscCall(PetscFree(cp->d));
 99:   PetscCall(VecDestroy(&cp->work));
100:   PetscCall(PetscFree3(cp->a, cp->i, cp->j));
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: static PetscErrorCode PCDestroy_CP(PC pc)
105: {
106:   PC_CP *cp = (PC_CP *)pc->data;

108:   PetscFunctionBegin;
109:   PetscCall(PCReset_CP(pc));
110:   PetscCall(PetscFree(cp->d));
111:   PetscCall(PetscFree3(cp->a, cp->i, cp->j));
112:   PetscCall(PetscFree(pc->data));
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

116: static PetscErrorCode PCSetFromOptions_CP(PC pc, PetscOptionItems *PetscOptionsObject)
117: {
118:   PetscFunctionBegin;
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: /*MC
123:      PCCP - a "column-projection" preconditioner

125:      This is a terrible preconditioner and is not recommended, ever!

127:      Loops over the entries of x computing dx_i (e_i is the unit vector in the ith direction) to
128: .vb

130:         min || b - A(x + dx_i e_i ||_2
131:         dx_i

133:     That is, it changes a single entry of x to minimize the new residual norm.
134:    Let A_i represent the ith column of A, then the minimization can be written as

136:        min || r - (dx_i) A e_i ||_2
137:        dx_i
138:    or   min || r - (dx_i) A_i ||_2
139:         dx_i

141:     take the derivative with respect to dx_i to obtain
142:         dx_i = (A_i^T A_i)^(-1) A_i^T r

144:     This algorithm can be thought of as Gauss-Seidel on the normal equations
145: .ve

147:     Notes:
148:     This procedure can also be done with block columns or any groups of columns
149:     but this is not coded.

151:     These "projections" can be done simultaneously for all columns (similar to Jacobi)
152:     or sequentially (similar to Gauss-Seidel/SOR). This is only coded for SOR type.

154:     This is related to, but not the same as "row projection" methods.

156:     This is currently coded only for `MATSEQAIJ` matrices

158:   Level: intermediate

160: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCJACOBI`, `PCSOR`
161: M*/

163: PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc)
164: {
165:   PC_CP *cp;

167:   PetscFunctionBegin;
168:   PetscCall(PetscNew(&cp));
169:   pc->data = (void *)cp;

171:   pc->ops->apply           = PCApply_CP;
172:   pc->ops->applytranspose  = PCApply_CP;
173:   pc->ops->setup           = PCSetUp_CP;
174:   pc->ops->reset           = PCReset_CP;
175:   pc->ops->destroy         = PCDestroy_CP;
176:   pc->ops->setfromoptions  = PCSetFromOptions_CP;
177:   pc->ops->view            = NULL;
178:   pc->ops->applyrichardson = NULL;
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }