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) {
 38:     PetscMalloc3(aij->nz,&cp->a,cp->n+1,&cp->i,aij->nz,&cp->j);
 39:   }
 40:   PetscCalloc1(cp->n,&colcnt);

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

 54:   /* compute sum of squares of each column d[] */
 55:   for (i=0; i<cp->n; i++) {  /* over columns */
 56:     cp->d[i] = 0.;
 57:     for (j=cp->i[i]; j<cp->i[i+1]; j++) cp->d[i] += cp->a[j]*cp->a[j]; /* over rows in column */
 58:     cp->d[i] = 1.0/cp->d[i];
 59:   }
 60:   return 0;
 61: }
 62: /* -------------------------------------------------------------------------- */
 63: static PetscErrorCode PCApply_CP(PC pc,Vec bb,Vec xx)
 64: {
 65:   PC_CP          *cp = (PC_CP*)pc->data;
 66:   PetscScalar    *b,*x,xt;
 67:   PetscInt       i,j;

 69:   VecCopy(bb,cp->work);
 70:   VecGetArray(cp->work,&b);
 71:   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:   VecRestoreArray(cp->work,&b);
 89:   VecRestoreArray(xx,&x);
 90:   return 0;
 91: }
 92: /* -------------------------------------------------------------------------- */
 93: static PetscErrorCode PCReset_CP(PC pc)
 94: {
 95:   PC_CP          *cp = (PC_CP*)pc->data;

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

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

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

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

119: /*MC
120:      PCCP - a "column-projection" preconditioner

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

124:      Loops over the entries of x computing dx_i (e_i is the unit vector in the ith direction) to
125: $
126: $        min || b - A(x + dx_i e_i ||_2
127: $        dx_i
128: $
129: $    That is, it changes a single entry of x to minimize the new residual norm.
130: $   Let A_i represent the ith column of A, then the minimization can be written as
131: $
132: $       min || r - (dx_i) A e_i ||_2
133: $       dx_i
134: $   or   min || r - (dx_i) A_i ||_2
135: $        dx_i
136: $
137: $    take the derivative with respect to dx_i to obtain
138: $        dx_i = (A_i^T A_i)^(-1) A_i^T r
139: $
140: $    This algorithm can be thought of as Gauss-Seidel on the normal equations

142:     Notes:
143:     This proceedure 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 SeqAIJ matrices in sequential (SOR) form.

153:   Level: intermediate

155: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCJACOBI, PCSOR

157: M*/

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

163:   PetscNewLog(pc,&cp);
164:   pc->data = (void*)cp;

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