1 #include <../src/ksp/pc/impls/deflation/deflation.h> 2 3 #if defined(HAVE_SLEPC) 4 #include <slepceps.h> 5 #include <slepcbv.h> 6 #endif 7 8 PetscScalar db2[] = {0.7071067811865476,0.7071067811865476}; 9 10 //dec low high 11 PetscScalar db4[] = {-0.12940952255092145,0.22414386804185735,0.836516303737469,0.48296291314469025}; 12 //rec low high 13 //PetscScalar db4[] = {0.48296291314469025,0.836516303737469,0.22414386804185735,-0.12940952255092145}; 14 15 PetscScalar db8[] = {-0.010597401784997278, 16 0.032883011666982945, 17 0.030841381835986965, 18 -0.18703481171888114, 19 -0.02798376941698385, 20 0.6308807679295904, 21 0.7148465705525415, 22 0.23037781330885523}; 23 24 //PetscScalar db8[] = {0.23037781330885523, 25 //0.7148465705525415, 26 //0.6308807679295904, 27 //-0.02798376941698385, 28 //-0.18703481171888114, 29 //0.030841381835986965, 30 //0.032883011666982945, 31 //-0.010597401784997278}; 32 33 PetscScalar db16[] = {-0.00011747678400228192, 34 0.0006754494059985568, 35 -0.0003917403729959771, 36 -0.00487035299301066, 37 0.008746094047015655, 38 0.013981027917015516, 39 -0.04408825393106472, 40 -0.01736930100202211, 41 0.128747426620186, 42 0.00047248457399797254, 43 -0.2840155429624281, 44 -0.015829105256023893, 45 0.5853546836548691, 46 0.6756307362980128, 47 0.3128715909144659, 48 0.05441584224308161}; 49 50 //PetscScalar db16[] = {0.05441584224308161, 51 //0.3128715909144659, 52 //0.6756307362980128, 53 //0.5853546836548691, 54 //-0.015829105256023893, 55 //-0.2840155429624281, 56 //0.00047248457399797254, 57 //0.128747426620186, 58 //-0.01736930100202211, 59 //-0.04408825393106472, 60 //0.013981027917015516, 61 //0.008746094047015655, 62 //-0.00487035299301066, 63 //-0.0003917403729959771, 64 //0.0006754494059985568, 65 //-0.00011747678400228192}; 66 67 PetscScalar biorth22[] = {0.0, 68 -0.1767766952966369, 69 0.3535533905932738, 70 1.0606601717798214, 71 0.3535533905932738, 72 -0.1767766952966369}; 73 74 //PetscScalar biorth22[] = {0.0, 75 //0.3535533905932738, 76 //0.7071067811865476, 77 //0.3535533905932738, 78 //0.0, 79 //0.0}; 80 81 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}; 82 83 static PetscErrorCode PCDeflationCreateSpaceWave(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc,Mat *H) 84 { 85 PetscErrorCode ierr; 86 Mat defl; 87 PetscInt i,j,k,ilo,ihi,*Iidx; 88 89 PetscFunctionBegin; 90 ierr = PetscMalloc1(ncoeffs,&Iidx);CHKERRQ(ierr); 91 92 /* TODO pass A instead of PC? */ 93 ierr = MatCreate(comm,&defl);CHKERRQ(ierr); 94 ierr = MatSetSizes(defl,m,n,M,N);CHKERRQ(ierr); 95 ierr = MatSetUp(defl);CHKERRQ(ierr); 96 ierr = MatSeqAIJSetPreallocation(defl,ncoeffs,NULL);CHKERRQ(ierr); 97 ierr = MatMPIAIJSetPreallocation(defl,ncoeffs,NULL,ncoeffs,NULL);CHKERRQ(ierr); 98 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 99 ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 100 101 /* Alg 735 Taswell: fvecmat */ 102 k = ncoeffs -2; 103 if (trunc) k = k/2; 104 105 ierr = MatGetOwnershipRange(defl,&ilo,&ihi);CHKERRQ(ierr); 106 for (i=0; i<ncoeffs; i++) { 107 Iidx[i] = i+ilo*2 -k; 108 if (Iidx[i] >= N) Iidx[i] = PETSC_MIN_INT; 109 } 110 for (i=ilo; i<ihi; i++) { 111 ierr = MatSetValues(defl,1,&i,ncoeffs,Iidx,coeffs,INSERT_VALUES);CHKERRQ(ierr); 112 for (j=0; j<ncoeffs; j++) { 113 Iidx[j] += 2; 114 if (Iidx[j] >= N) Iidx[j] = PETSC_MIN_INT; 115 } 116 } 117 118 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 119 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 120 121 ierr = PetscFree(Iidx);CHKERRQ(ierr); 122 *H = defl; 123 PetscFunctionReturn(0); 124 } 125 126 PetscErrorCode PCDeflationGetSpaceHaar(PC pc,Mat *W,PetscInt size) 127 { 128 PetscErrorCode ierr; 129 Mat A,defl; 130 PetscInt i,j,len,ilo,ihi,*Iidx,m,M; 131 PetscReal *col,val; 132 133 PetscFunctionBegin; 134 /* Haar basis wavelet, level=size */ 135 len = pow(2,size); 136 ierr = PetscMalloc1(len,&col);CHKERRQ(ierr); 137 ierr = PetscMalloc1(len,&Iidx);CHKERRQ(ierr); 138 val = 1./pow(2,size/2.); 139 for (i=0; i<len; i++) col[i] = val; 140 141 /* TODO pass A instead of PC? */ 142 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 143 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 144 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 145 ierr = MatCreate(PetscObjectComm((PetscObject)A),&defl);CHKERRQ(ierr); 146 ierr = MatSetSizes(defl,m,PETSC_DECIDE,M,ceil(M/(float)len));CHKERRQ(ierr); 147 ierr = MatSetUp(defl);CHKERRQ(ierr); 148 ierr = MatSeqAIJSetPreallocation(defl,size,NULL);CHKERRQ(ierr); 149 ierr = MatMPIAIJSetPreallocation(defl,size,NULL,size,NULL);CHKERRQ(ierr); 150 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 151 152 ierr = MatGetOwnershipRangeColumn(defl,&ilo,&ihi);CHKERRQ(ierr); 153 for (i=0; i<len; i++) Iidx[i] = i+ilo*len; 154 if (M%len && ihi == (int)ceil(M/(float)len)) ihi -= 1; 155 for (i=ilo; i<ihi; i++) { 156 ierr = MatSetValues(defl,len,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr); 157 for (j=0; j<len; j++) Iidx[j] += len; 158 } 159 if (M%len && ihi+1 == ceil(M/(float)len)) { 160 len = M%len; 161 val = 1./pow(pow(2,len),0.5); 162 for (i=0; i<len; i++) col[i] = val; 163 ierr = MatSetValues(defl,len,Iidx,1,&ihi,col,INSERT_VALUES);CHKERRQ(ierr); 164 } 165 166 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 167 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 168 169 ierr = PetscFree(col);CHKERRQ(ierr); 170 ierr = PetscFree(Iidx);CHKERRQ(ierr); 171 *W = defl; 172 PetscFunctionReturn(0); 173 } 174 175 PetscErrorCode PCDeflationGetSpaceWave(PC pc,Mat *W,PetscInt size,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc) 176 { 177 PetscErrorCode ierr; 178 Mat A,*H,defl; 179 PetscInt i,m,M,Mdefl,Ndefl; 180 MPI_Comm comm; 181 182 PetscFunctionBegin; 183 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 184 ierr = PetscMalloc1(size,&H);CHKERRQ(ierr); 185 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 186 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 187 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 188 Mdefl = M; 189 Ndefl = M; 190 for (i=0; i<size; i++) { 191 if (Mdefl%2) { 192 if (trunc) { 193 Mdefl = (PetscInt)PetscCeilReal(Mdefl/2.); 194 } else { 195 Mdefl = (PetscInt)PetscFloorReal((ncoeffs+Mdefl-1)/2.); 196 } 197 } else { 198 Mdefl = Mdefl/2; 199 } 200 ierr = PCDeflationCreateSpaceWave(comm,PETSC_DECIDE,m,Mdefl,Ndefl,ncoeffs,coeffs,trunc,&H[i]);CHKERRQ(ierr); 201 ierr = MatGetLocalSize(H[i],&m,NULL);CHKERRQ(ierr); 202 Ndefl = Mdefl; 203 } 204 ierr = MatCreateComposite(comm,size,H,&defl);CHKERRQ(ierr); 205 ierr = MatCompositeSetType(defl,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 206 *W = defl; 207 208 for (i=0; i<size; i++) { 209 ierr = MatDestroy(&H[i]);CHKERRQ(ierr); 210 } 211 ierr = PetscFree(H);CHKERRQ(ierr); 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode PCDeflationGetSpaceAggregation(PC pc,Mat *W) 216 { 217 PetscErrorCode ierr; 218 Mat A,defl; 219 PetscInt i,ilo,ihi,*Iidx,m,M; 220 PetscReal *col; 221 MPI_Comm comm; 222 223 PetscFunctionBegin; 224 /* TODO pass A instead of PC? */ 225 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 226 ierr = MatGetOwnershipRangeColumn(A,&ilo,&ihi);CHKERRQ(ierr); 227 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 228 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 229 ierr = MPI_Comm_size(comm,&m);CHKERRQ(ierr); 230 ierr = MatCreate(comm,&defl);CHKERRQ(ierr); 231 ierr = MatSetSizes(defl,ihi-ilo,1,M,m);CHKERRQ(ierr); 232 ierr = MatSetUp(defl);CHKERRQ(ierr); 233 ierr = MatSeqAIJSetPreallocation(defl,1,NULL);CHKERRQ(ierr); 234 ierr = MatMPIAIJSetPreallocation(defl,1,NULL,0,NULL);CHKERRQ(ierr); 235 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 236 ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 237 238 239 ierr = PetscMalloc1(ihi-ilo,&col);CHKERRQ(ierr); 240 ierr = PetscMalloc1(ihi-ilo,&Iidx);CHKERRQ(ierr); 241 for (i=ilo; i<ihi; i++) { 242 Iidx[i-ilo] = i; 243 col[i-ilo] = 1; 244 } 245 ierr = MPI_Comm_rank(comm,&i);CHKERRQ(ierr); 246 ierr = MatSetValues(defl,ihi-ilo,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr); 247 248 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 249 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 250 251 ierr = PetscFree(col);CHKERRQ(ierr); 252 ierr = PetscFree(Iidx);CHKERRQ(ierr); 253 *W = defl; 254 PetscFunctionReturn(0); 255 } 256 257 PetscErrorCode PCDeflationGetSpaceSLEPc(PC pc,Mat *W,PetscInt size,PetscBool cheapCP) 258 { 259 #if defined(HAVE_SLEPC) 260 PetscErrorCode ierr; 261 PC_Deflation *def = (PC_Deflation*)pc->data; 262 Mat A,defl; 263 Vec vec; 264 EPS eps; 265 PetscScalar *data,*dataScaled,eigval; 266 PetscInt i,nconv,m,M,n=PETSC_DECIDE; 267 PetscBool slepcinit; 268 MPI_Comm comm; 269 270 PetscFunctionBegin; 271 ierr = SlepcInitialized(&slepcinit);CHKERRQ(ierr); 272 if (!slepcinit) { 273 ierr = SlepcInitialize(NULL,NULL,(char*)0,(char*)0);CHKERRQ(ierr); 274 slepcinit = PETSC_TRUE; 275 } 276 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 277 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 278 ierr = EPSCreate(comm,&eps);CHKERRQ(ierr); 279 ierr = EPSSetOperators(eps,A,NULL);CHKERRQ(ierr); 280 ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr); /* Implemented only for def */ 281 ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr); 282 ierr = EPSSetDimensions(eps,size,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 283 ierr = EPSSetFromOptions(eps);CHKERRQ(ierr); 284 285 ierr = EPSSolve(eps);CHKERRQ(ierr); 286 ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr); 287 if (nconv < size) SETERRQ2(comm,PETSC_ERR_CONV_FAILED,"SLEPc: Number of converged eigenpairs (%d) is less than requested (%d)",nconv,size); 288 ierr = MatCreateVecs(A,NULL,&vec);CHKERRQ(ierr); 289 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 290 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 291 ierr = PetscSplitOwnership(comm,&n,&size);CHKERRQ(ierr); 292 ierr = PetscMalloc1(m*size,&data);CHKERRQ(ierr); 293 /* TODO check that eigenvalue is not 0 -> vec is not in Ker A */ 294 for (i=0; i<size; i++) { 295 ierr = VecPlaceArray(vec,&data[i*m]);CHKERRQ(ierr); 296 ierr = EPSGetEigenvector(eps,i,vec,NULL);CHKERRQ(ierr); 297 ierr = VecResetArray(vec);CHKERRQ(ierr); 298 } 299 ierr = MatCreateDense(comm,m,n,M,size,data,&defl);CHKERRQ(ierr); 300 301 if (cheapCP) { 302 ierr = PetscMalloc1(m*size,&dataScaled);CHKERRQ(ierr); 303 for (i=0; i<size; i++) { 304 ierr = VecPlaceArray(vec,&dataScaled[i*m]);CHKERRQ(ierr); 305 ierr = EPSGetEigenpair(eps,i,&eigval,NULL,vec,NULL);CHKERRQ(ierr); 306 ierr = VecScale(vec,eigval);CHKERRQ(ierr); 307 ierr = VecResetArray(vec);CHKERRQ(ierr); 308 } 309 ierr = MatCreateDense(comm,m,n,M,size,dataScaled,&def->AW);CHKERRQ(ierr); 310 } 311 312 //ierr = EPSGetBV(eps,&bv);CHKERRQ(ierr); 313 //ierr = BVCreateMat(bv,&defl);CHKERRQ(ierr); 314 *W = defl; 315 316 ierr = EPSDestroy(&eps);CHKERRQ(ierr); 317 if (slepcinit) ierr = SlepcFinalize();CHKERRQ(ierr); 318 PetscFunctionReturn(0); 319 #else 320 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_CONV_FAILED,"Not compiled with SLEPc support (call make HAVE_SLEPC)"); 321 #endif 322 } 323 324 PetscErrorCode PCDeflationComputeSpace(PC pc) 325 { 326 PetscErrorCode ierr; 327 Mat defl; 328 PetscBool transp=PETSC_TRUE; 329 PC_Deflation *def = (PC_Deflation*)pc->data; 330 331 /* TODO valid header */ 332 PetscFunctionBegin; 333 if (def->spacesize < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space size specified: %d",def->spacesize); 334 switch (def->spacetype) { 335 case PC_DEFLATION_SPACE_HAAR: 336 transp = PETSC_FALSE; 337 ierr = PCDeflationGetSpaceHaar(pc,&defl,def->spacesize);CHKERRQ(ierr);break; 338 case PC_DEFLATION_SPACE_DB2: 339 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,2,db2,!def->extendsp);CHKERRQ(ierr);break; 340 case PC_DEFLATION_SPACE_DB4: 341 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,4,db4,!def->extendsp);CHKERRQ(ierr);break; 342 case PC_DEFLATION_SPACE_DB8: 343 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,8,db8,!def->extendsp);CHKERRQ(ierr);break; 344 case PC_DEFLATION_SPACE_DB16: 345 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,16,db16,!def->extendsp);CHKERRQ(ierr);break; 346 case PC_DEFLATION_SPACE_BIORTH22: 347 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,6,biorth22,!def->extendsp);CHKERRQ(ierr);break; 348 case PC_DEFLATION_SPACE_MEYER: 349 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,62,meyer,!def->extendsp);CHKERRQ(ierr);break; 350 case PC_DEFLATION_SPACE_AGGREGATION: 351 transp = PETSC_FALSE; 352 ierr = PCDeflationGetSpaceAggregation(pc,&defl);CHKERRQ(ierr);break; 353 case PC_DEFLATION_SPACE_SLEPC: 354 transp = PETSC_FALSE; 355 ierr = PCDeflationGetSpaceSLEPc(pc,&defl,def->spacesize,PETSC_FALSE);CHKERRQ(ierr);break; 356 case PC_DEFLATION_SPACE_SLEPC_CHEAP: 357 transp = PETSC_FALSE; 358 ierr = PCDeflationGetSpaceSLEPc(pc,&defl,def->spacesize,PETSC_TRUE);CHKERRQ(ierr);break; 359 default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space Type specified"); 360 } 361 362 ierr = PCDeflationSetSpace(pc,defl,transp);CHKERRQ(ierr); 363 ierr = MatDestroy(&defl);CHKERRQ(ierr); 364 PetscFunctionReturn(0); 365 } 366 367