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(PetscOptions *PetscOptionsObject,PC pc) 134 { 135 PetscFunctionBegin; 136 PetscFunctionReturn(0); 137 } 138 139 140 /*MC 141 PCCP - a "column-projection" preconditioner 142 143 This is a terrible preconditioner and is not recommended, ever! 144 145 Loops over the entries of x computing dx_i to 146 $ 147 $ min || b - A(x + dx_i e_i ||_2 148 $ dx_i 149 $ 150 $ That is, it changes a single entry of x to minimize the new residual. 151 $ Let A_i represent the ith column of A, then the minimization can be written as 152 $ 153 $ min || r - (dx_i) A e_i ||_2 154 $ dx_i 155 $ or min || r - (dx_i) A_i ||_2 156 $ dx_i 157 $ 158 $ take the derivative with respect to dx_i to obtain 159 $ dx_i = (A_i^T A_i)^(-1) A_i^T r 160 $ 161 $ This algorithm can be thought of as Gauss-Seidel on the normal equations 162 163 Notes: This proceedure can also be done with block columns or any groups of columns 164 but this is not coded. 165 166 These "projections" can be done simultaneously for all columns (similar to Jacobi) 167 or sequentially (similar to Gauss-Seidel/SOR). This is only coded for SOR type. 168 169 This is related to, but not the same as "row projection" methods. 170 171 This is currently coded only for SeqAIJ matrices in sequential (SOR) form. 172 173 Level: intermediate 174 175 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCJACOBI, PCSOR 176 177 M*/ 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "PCCreate_CP" 181 PETSC_EXTERN PetscErrorCode PCCreate_CP(PC pc) 182 { 183 PC_CP *cp; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 ierr = PetscNewLog(pc,&cp);CHKERRQ(ierr); 188 pc->data = (void*)cp; 189 190 pc->ops->apply = PCApply_CP; 191 pc->ops->applytranspose = PCApply_CP; 192 pc->ops->setup = PCSetUp_CP; 193 pc->ops->reset = PCReset_CP; 194 pc->ops->destroy = PCDestroy_CP; 195 pc->ops->setfromoptions = PCSetFromOptions_CP; 196 pc->ops->view = 0; 197 pc->ops->applyrichardson = 0; 198 PetscFunctionReturn(0); 199 } 200 201 202