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:   PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg);

 26:   MatGetLocalSize(pc->pmat, &cp->m, &cp->n);

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

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

 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:   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:   PetscFree(colcnt);

 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:   return 0;
 59: }

 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;

 67:   VecCopy(bb, cp->work);
 68:   VecGetArray(cp->work, &b);
 69:   VecGetArray(xx, &x);

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

 86:   VecRestoreArray(cp->work, &b);
 87:   VecRestoreArray(xx, &x);
 88:   return 0;
 89: }

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

 95:   PetscFree(cp->d);
 96:   VecDestroy(&cp->work);
 97:   PetscFree3(cp->a, cp->i, cp->j);
 98:   return 0;
 99: }

101: static PetscErrorCode PCDestroy_CP(PC pc)
102: {
103:   PC_CP *cp = (PC_CP *)pc->data;

105:   PCReset_CP(pc);
106:   PetscFree(cp->d);
107:   PetscFree3(cp->a, cp->i, cp->j);
108:   PetscFree(pc->data);
109:   return 0;
110: }

112: static PetscErrorCode PCSetFromOptions_CP(PC pc, PetscOptionItems *PetscOptionsObject)
113: {
114:   return 0;
115: }

117: /*MC
118:      PCCP - a "column-projection" preconditioner

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

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

125:         min || b - A(x + dx_i e_i ||_2
126:         dx_i

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

131:        min || r - (dx_i) A e_i ||_2
132:        dx_i
133:    or   min || r - (dx_i) A_i ||_2
134:         dx_i

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

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

142:     Notes:
143:     This procedure can also be done with block columns or any groups of columns
144:     but this is not coded.

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

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

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

153:   Level: intermediate

155: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCJACOBI`, `PCSOR`
156: M*/

158: PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc)
159: {
160:   PC_CP *cp;

162:   PetscNew(&cp);
163:   pc->data = (void *)cp;

165:   pc->ops->apply           = PCApply_CP;
166:   pc->ops->applytranspose  = PCApply_CP;
167:   pc->ops->setup           = PCSetUp_CP;
168:   pc->ops->reset           = PCReset_CP;
169:   pc->ops->destroy         = PCDestroy_CP;
170:   pc->ops->setfromoptions  = PCSetFromOptions_CP;
171:   pc->ops->view            = NULL;
172:   pc->ops->applyrichardson = NULL;
173:   return 0;
174: }