1 #include <../src/ksp/pc/impls/deflation/deflation.h> /*I "petscksp.h" I*/ 2 3 PetscScalar db2[] = {0.7071067811865476,0.7071067811865476}; 4 5 PetscScalar db4[] = {-0.12940952255092145,0.22414386804185735,0.836516303737469,0.48296291314469025}; 6 7 PetscScalar db8[] = {-0.010597401784997278, 8 0.032883011666982945, 9 0.030841381835986965, 10 -0.18703481171888114, 11 -0.02798376941698385, 12 0.6308807679295904, 13 0.7148465705525415, 14 0.23037781330885523}; 15 16 PetscScalar db16[] = {-0.00011747678400228192, 17 0.0006754494059985568, 18 -0.0003917403729959771, 19 -0.00487035299301066, 20 0.008746094047015655, 21 0.013981027917015516, 22 -0.04408825393106472, 23 -0.01736930100202211, 24 0.128747426620186, 25 0.00047248457399797254, 26 -0.2840155429624281, 27 -0.015829105256023893, 28 0.5853546836548691, 29 0.6756307362980128, 30 0.3128715909144659, 31 0.05441584224308161}; 32 33 PetscScalar biorth22[] = {0.0, 34 -0.1767766952966369, 35 0.3535533905932738, 36 1.0606601717798214, 37 0.3535533905932738, 38 -0.1767766952966369}; 39 40 PetscScalar meyer[] = {0.0,-1.009999956941423e-12,8.519459636796214e-09,-1.111944952595278e-08,-1.0798819539621958e-08,6.066975741351135e-08,-1.0866516536735883e-07,8.200680650386481e-08,1.1783004497663934e-07,-5.506340565252278e-07,1.1307947017916706e-06,-1.489549216497156e-06,7.367572885903746e-07,3.20544191334478e-06,-1.6312699734552807e-05,6.554305930575149e-05,-0.0006011502343516092,-0.002704672124643725,0.002202534100911002,0.006045814097323304,-0.006387718318497156,-0.011061496392513451,0.015270015130934803,0.017423434103729693,-0.03213079399021176,-0.024348745906078023,0.0637390243228016,0.030655091960824263,-0.13284520043622938,-0.035087555656258346,0.44459300275757724,0.7445855923188063,0.44459300275757724,-0.035087555656258346,-0.13284520043622938,0.030655091960824263,0.0637390243228016,-0.024348745906078023,-0.03213079399021176,0.017423434103729693,0.015270015130934803,-0.011061496392513451,-0.006387718318497156,0.006045814097323304,0.002202534100911002,-0.002704672124643725,-0.0006011502343516092,6.554305930575149e-05,-1.6312699734552807e-05,3.20544191334478e-06,7.367572885903746e-07,-1.489549216497156e-06,1.1307947017916706e-06,-5.506340565252278e-07,1.1783004497663934e-07,8.200680650386481e-08,-1.0866516536735883e-07,6.066975741351135e-08,-1.0798819539621958e-08,-1.111944952595278e-08,8.519459636796214e-09,-1.009999956941423e-12}; 41 42 static PetscErrorCode PCDeflationCreateSpaceWave(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc,Mat *H) 43 { 44 PetscErrorCode ierr; 45 Mat defl; 46 PetscInt i,j,k,ilo,ihi,*Iidx; 47 48 PetscFunctionBegin; 49 ierr = PetscMalloc1(ncoeffs,&Iidx);CHKERRQ(ierr); 50 51 ierr = MatCreate(comm,&defl);CHKERRQ(ierr); 52 ierr = MatSetSizes(defl,m,n,M,N);CHKERRQ(ierr); 53 ierr = MatSetUp(defl);CHKERRQ(ierr); 54 ierr = MatSeqAIJSetPreallocation(defl,ncoeffs,NULL);CHKERRQ(ierr); 55 ierr = MatMPIAIJSetPreallocation(defl,ncoeffs,NULL,ncoeffs,NULL);CHKERRQ(ierr); 56 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 57 ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 58 59 /* Alg 735 Taswell: fvecmat */ 60 k = ncoeffs -2; 61 if (trunc) k = k/2; 62 63 ierr = MatGetOwnershipRange(defl,&ilo,&ihi);CHKERRQ(ierr); 64 for (i=0; i<ncoeffs; i++) { 65 Iidx[i] = i+ilo*2 -k; 66 if (Iidx[i] >= N) Iidx[i] = PETSC_MIN_INT; 67 } 68 for (i=ilo; i<ihi; i++) { 69 ierr = MatSetValues(defl,1,&i,ncoeffs,Iidx,coeffs,INSERT_VALUES);CHKERRQ(ierr); 70 for (j=0; j<ncoeffs; j++) { 71 Iidx[j] += 2; 72 if (Iidx[j] >= N) Iidx[j] = PETSC_MIN_INT; 73 } 74 } 75 76 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 77 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 78 79 ierr = PetscFree(Iidx);CHKERRQ(ierr); 80 *H = defl; 81 PetscFunctionReturn(0); 82 } 83 84 PetscErrorCode PCDeflationGetSpaceHaar(PC pc,Mat *W,PetscInt size) 85 { 86 PetscErrorCode ierr; 87 Mat A,defl; 88 PetscInt i,j,len,ilo,ihi,*Iidx,m,M; 89 PetscScalar *col,val; 90 91 PetscFunctionBegin; 92 /* Haar basis wavelet, level=size */ 93 len = pow(2,size); 94 ierr = PetscMalloc1(len,&col);CHKERRQ(ierr); 95 ierr = PetscMalloc1(len,&Iidx);CHKERRQ(ierr); 96 val = 1./pow(2,size/2.); 97 for (i=0; i<len; i++) col[i] = val; 98 99 ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr); 100 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 101 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 102 ierr = MatCreate(PetscObjectComm((PetscObject)A),&defl);CHKERRQ(ierr); 103 ierr = MatSetSizes(defl,m,PETSC_DECIDE,M,(PetscInt)ceil(M/(float)len));CHKERRQ(ierr); 104 ierr = MatSetUp(defl);CHKERRQ(ierr); 105 ierr = MatSeqAIJSetPreallocation(defl,size,NULL);CHKERRQ(ierr); 106 ierr = MatMPIAIJSetPreallocation(defl,size,NULL,size,NULL);CHKERRQ(ierr); 107 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 108 109 ierr = MatGetOwnershipRangeColumn(defl,&ilo,&ihi);CHKERRQ(ierr); 110 for (i=0; i<len; i++) Iidx[i] = i+ilo*len; 111 if (M%len && ihi == (int)ceil(M/(float)len)) ihi -= 1; 112 for (i=ilo; i<ihi; i++) { 113 ierr = MatSetValues(defl,len,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr); 114 for (j=0; j<len; j++) Iidx[j] += len; 115 } 116 if (M%len && ihi+1 == ceil(M/(float)len)) { 117 len = M%len; 118 val = 1./pow(pow(2,len),0.5); 119 for (i=0; i<len; i++) col[i] = val; 120 ierr = MatSetValues(defl,len,Iidx,1,&ihi,col,INSERT_VALUES);CHKERRQ(ierr); 121 } 122 123 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 124 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 125 126 ierr = PetscFree(col);CHKERRQ(ierr); 127 ierr = PetscFree(Iidx);CHKERRQ(ierr); 128 *W = defl; 129 PetscFunctionReturn(0); 130 } 131 132 PetscErrorCode PCDeflationGetSpaceWave(PC pc,Mat *W,PetscInt size,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc) 133 { 134 PetscErrorCode ierr; 135 Mat A,*H,defl; 136 PetscInt i,m,M,Mdefl,Ndefl; 137 MPI_Comm comm; 138 139 PetscFunctionBegin; 140 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 141 ierr = PetscMalloc1(size,&H);CHKERRQ(ierr); 142 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); 143 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 144 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 145 Mdefl = M; 146 Ndefl = M; 147 for (i=0; i<size; i++) { 148 if (Mdefl%2) { 149 if (trunc) { 150 Mdefl = (PetscInt)PetscCeilReal(Mdefl/2.); 151 } else { 152 Mdefl = (PetscInt)PetscFloorReal((ncoeffs+Mdefl-1)/2.); 153 } 154 } else { 155 Mdefl = Mdefl/2; 156 } 157 ierr = PCDeflationCreateSpaceWave(comm,PETSC_DECIDE,m,Mdefl,Ndefl,ncoeffs,coeffs,trunc,&H[i]);CHKERRQ(ierr); 158 ierr = MatGetLocalSize(H[i],&m,NULL);CHKERRQ(ierr); 159 Ndefl = Mdefl; 160 } 161 ierr = MatCreateComposite(comm,size,H,&defl);CHKERRQ(ierr); 162 ierr = MatCompositeSetType(defl,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 163 *W = defl; 164 165 for (i=0; i<size; i++) { 166 ierr = MatDestroy(&H[i]);CHKERRQ(ierr); 167 } 168 ierr = PetscFree(H);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 PetscErrorCode PCDeflationGetSpaceAggregation(PC pc,Mat *W) 173 { 174 PetscErrorCode ierr; 175 Mat A,defl; 176 PetscInt i,ilo,ihi,*Iidx,m,M; 177 PetscScalar *col; 178 MPI_Comm comm; 179 180 PetscFunctionBegin; 181 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); 182 ierr = MatGetOwnershipRangeColumn(A,&ilo,&ihi);CHKERRQ(ierr); 183 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 184 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 185 ierr = MPI_Comm_size(comm,&m);CHKERRQ(ierr); 186 ierr = MatCreate(comm,&defl);CHKERRQ(ierr); 187 ierr = MatSetSizes(defl,ihi-ilo,1,M,m);CHKERRQ(ierr); 188 ierr = MatSetUp(defl);CHKERRQ(ierr); 189 ierr = MatSeqAIJSetPreallocation(defl,1,NULL);CHKERRQ(ierr); 190 ierr = MatMPIAIJSetPreallocation(defl,1,NULL,0,NULL);CHKERRQ(ierr); 191 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 192 ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 193 194 ierr = PetscMalloc1(ihi-ilo,&col);CHKERRQ(ierr); 195 ierr = PetscMalloc1(ihi-ilo,&Iidx);CHKERRQ(ierr); 196 for (i=ilo; i<ihi; i++) { 197 Iidx[i-ilo] = i; 198 col[i-ilo] = 1; 199 } 200 ierr = MPI_Comm_rank(comm,&i);CHKERRQ(ierr); 201 ierr = MatSetValues(defl,ihi-ilo,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr); 202 203 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 204 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 205 206 ierr = PetscFree(col);CHKERRQ(ierr); 207 ierr = PetscFree(Iidx);CHKERRQ(ierr); 208 *W = defl; 209 PetscFunctionReturn(0); 210 } 211 212 PetscErrorCode PCDeflationComputeSpace(PC pc) 213 { 214 PetscErrorCode ierr; 215 Mat defl; 216 PetscBool transp=PETSC_TRUE; 217 PC_Deflation *def = (PC_Deflation*)pc->data; 218 219 PetscFunctionBegin; 220 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 221 if (def->spacesize < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PCDeflation space size specified: %D",def->spacesize); 222 switch (def->spacetype) { 223 case PC_DEFLATION_SPACE_HAAR: 224 transp = PETSC_FALSE; 225 ierr = PCDeflationGetSpaceHaar(pc,&defl,def->spacesize);CHKERRQ(ierr);break; 226 case PC_DEFLATION_SPACE_DB2: 227 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,2,db2,!def->extendsp);CHKERRQ(ierr);break; 228 case PC_DEFLATION_SPACE_DB4: 229 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,4,db4,!def->extendsp);CHKERRQ(ierr);break; 230 case PC_DEFLATION_SPACE_DB8: 231 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,8,db8,!def->extendsp);CHKERRQ(ierr);break; 232 case PC_DEFLATION_SPACE_DB16: 233 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,16,db16,!def->extendsp);CHKERRQ(ierr);break; 234 case PC_DEFLATION_SPACE_BIORTH22: 235 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,6,biorth22,!def->extendsp);CHKERRQ(ierr);break; 236 case PC_DEFLATION_SPACE_MEYER: 237 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,62,meyer,!def->extendsp);CHKERRQ(ierr);break; 238 case PC_DEFLATION_SPACE_AGGREGATION: 239 transp = PETSC_FALSE; 240 ierr = PCDeflationGetSpaceAggregation(pc,&defl);CHKERRQ(ierr);break; 241 default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PCDeflationSpaceType specified"); 242 } 243 244 ierr = PCDeflationSetSpace(pc,defl,transp);CHKERRQ(ierr); 245 ierr = MatDestroy(&defl);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249