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