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 const char *const PCDeflationSpaceTypes[] = { 9 "haar", 10 "jacket-haar", 11 "db2", 12 "db4", 13 "db8", 14 "db16", 15 "biorth22", 16 "meyer", 17 "aggregation", 18 "slepc", 19 "slepc-cheap", 20 "user", 21 "DdefSpaceType", 22 "Ddef_SPACE_", 23 0 24 }; 25 26 PetscScalar db2[] = {0.7071067811865476,0.7071067811865476}; 27 28 //dec low high 29 PetscScalar db4[] = {-0.12940952255092145,0.22414386804185735,0.836516303737469,0.48296291314469025}; 30 //rec low high 31 //PetscScalar db4[] = {0.48296291314469025,0.836516303737469,0.22414386804185735,-0.12940952255092145}; 32 33 PetscScalar db8[] = {-0.010597401784997278, 34 0.032883011666982945, 35 0.030841381835986965, 36 -0.18703481171888114, 37 -0.02798376941698385, 38 0.6308807679295904, 39 0.7148465705525415, 40 0.23037781330885523}; 41 42 //PetscScalar db8[] = {0.23037781330885523, 43 //0.7148465705525415, 44 //0.6308807679295904, 45 //-0.02798376941698385, 46 //-0.18703481171888114, 47 //0.030841381835986965, 48 //0.032883011666982945, 49 //-0.010597401784997278}; 50 51 PetscScalar db16[] = {-0.00011747678400228192, 52 0.0006754494059985568, 53 -0.0003917403729959771, 54 -0.00487035299301066, 55 0.008746094047015655, 56 0.013981027917015516, 57 -0.04408825393106472, 58 -0.01736930100202211, 59 0.128747426620186, 60 0.00047248457399797254, 61 -0.2840155429624281, 62 -0.015829105256023893, 63 0.5853546836548691, 64 0.6756307362980128, 65 0.3128715909144659, 66 0.05441584224308161}; 67 68 //PetscScalar db16[] = {0.05441584224308161, 69 //0.3128715909144659, 70 //0.6756307362980128, 71 //0.5853546836548691, 72 //-0.015829105256023893, 73 //-0.2840155429624281, 74 //0.00047248457399797254, 75 //0.128747426620186, 76 //-0.01736930100202211, 77 //-0.04408825393106472, 78 //0.013981027917015516, 79 //0.008746094047015655, 80 //-0.00487035299301066, 81 //-0.0003917403729959771, 82 //0.0006754494059985568, 83 //-0.00011747678400228192}; 84 85 PetscScalar biorth22[] = {0.0, 86 -0.1767766952966369, 87 0.3535533905932738, 88 1.0606601717798214, 89 0.3535533905932738, 90 -0.1767766952966369}; 91 92 //PetscScalar biorth22[] = {0.0, 93 //0.3535533905932738, 94 //0.7071067811865476, 95 //0.3535533905932738, 96 //0.0, 97 //0.0}; 98 99 PetscReal 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}; 100 101 static PetscErrorCode PCDeflationCreateSpaceJacketHaar(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscBool jacket,Mat *H) 102 { 103 PetscErrorCode ierr; 104 Mat defl; 105 PetscInt i,j,ilo,ihi,alloc=2,*Iidx; 106 PetscReal val,*row; 107 108 PetscFunctionBegin; 109 if (jacket) alloc = 3; 110 ierr = PetscMalloc1(alloc,&row);CHKERRQ(ierr); 111 ierr = PetscMalloc1(alloc,&Iidx);CHKERRQ(ierr); 112 113 val = 1./pow(2,0.5); 114 row[0] = val; 115 row[1] = val; 116 117 /* TODO pass A instead of PC? */ 118 ierr = MatCreate(comm,&defl);CHKERRQ(ierr); 119 ierr = MatSetSizes(defl,m,n,M,N);CHKERRQ(ierr); 120 ierr = MatSetUp(defl);CHKERRQ(ierr); 121 ierr = MatSeqAIJSetPreallocation(defl,alloc,NULL);CHKERRQ(ierr); 122 ierr = MatMPIAIJSetPreallocation(defl,alloc,NULL,alloc,NULL);CHKERRQ(ierr); 123 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 124 ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 125 126 ierr = MatGetOwnershipRange(defl,&ilo,&ihi);CHKERRQ(ierr); 127 for (i=0; i<2; i++) Iidx[i] = i+ilo*2; 128 if (jacket && ihi==M) ihi -=2; 129 if (ihi<ilo) SETERRQ1(comm,PETSC_ERR_ARG_WRONG,"To many cores to assemble Jacket Haar matrix with %d rows",M); 130 for (i=ilo; i<ihi; i++) { 131 ierr = MatSetValues(defl,1,&i,2,Iidx,row,INSERT_VALUES);CHKERRQ(ierr); 132 for (j=0; j<2; j++) Iidx[j] += 2; 133 } 134 if (jacket && ihi == M-2) { 135 for (i=0; i<3; i++) Iidx[i] = i+ilo*2; 136 row[0] = 0.5; row[1] = 0.5; row[2] = val; 137 ierr = MatSetValues(defl,1,&ihi,3,Iidx,row,INSERT_VALUES);CHKERRQ(ierr); 138 ihi += 1; 139 row[2] = -row[2]; 140 ierr = MatSetValues(defl,1,&ihi,3,Iidx,row,INSERT_VALUES);CHKERRQ(ierr); 141 } 142 143 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 144 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 145 146 ierr = PetscFree(row);CHKERRQ(ierr); 147 ierr = PetscFree(Iidx);CHKERRQ(ierr); 148 *H = defl; 149 PetscFunctionReturn(0); 150 } 151 152 static PetscErrorCode PCDeflationCreateSpaceWave(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc,Mat *H) 153 { 154 PetscErrorCode ierr; 155 Mat defl; 156 PetscInt i,j,k,ilo,ihi,*Iidx; 157 158 PetscFunctionBegin; 159 ierr = PetscMalloc1(ncoeffs,&Iidx);CHKERRQ(ierr); 160 161 /* TODO pass A instead of PC? */ 162 ierr = MatCreate(comm,&defl);CHKERRQ(ierr); 163 ierr = MatSetSizes(defl,m,n,M,N);CHKERRQ(ierr); 164 ierr = MatSetUp(defl);CHKERRQ(ierr); 165 ierr = MatSeqAIJSetPreallocation(defl,ncoeffs,NULL);CHKERRQ(ierr); 166 ierr = MatMPIAIJSetPreallocation(defl,ncoeffs,NULL,ncoeffs,NULL);CHKERRQ(ierr); 167 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 168 ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 169 170 /* Alg 735 Taswell: fvecmat */ 171 k = ncoeffs -2; 172 if (trunc) k = k/2; 173 174 ierr = MatGetOwnershipRange(defl,&ilo,&ihi);CHKERRQ(ierr); 175 for (i=0; i<ncoeffs; i++) { 176 Iidx[i] = i+ilo*2 -k; 177 if (Iidx[i] >= N) Iidx[i] = PETSC_MIN_INT; 178 } 179 for (i=ilo; i<ihi; i++) { 180 ierr = MatSetValues(defl,1,&i,ncoeffs,Iidx,coeffs,INSERT_VALUES);CHKERRQ(ierr); 181 for (j=0; j<ncoeffs; j++) { 182 Iidx[j] += 2; 183 if (Iidx[j] >= N) Iidx[j] = PETSC_MIN_INT; 184 } 185 } 186 187 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 188 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 189 190 ierr = PetscFree(Iidx);CHKERRQ(ierr); 191 *H = defl; 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode PCDeflationGetSpaceHaar(PC pc,Mat *W,PetscInt size) 196 { 197 PetscErrorCode ierr; 198 Mat A,defl; 199 PetscInt i,j,len,ilo,ihi,*Iidx,m,M; 200 PetscReal *col,val; 201 202 PetscFunctionBegin; 203 /* Haar basis wavelet, level=size */ 204 len = pow(2,size); 205 ierr = PetscMalloc1(len,&col);CHKERRQ(ierr); 206 ierr = PetscMalloc1(len,&Iidx);CHKERRQ(ierr); 207 val = 1./pow(2,size/2.); 208 for (i=0; i<len; i++) col[i] = val; 209 210 /* TODO pass A instead of PC? */ 211 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 212 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 213 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 214 ierr = MatCreate(PetscObjectComm((PetscObject)A),&defl);CHKERRQ(ierr); 215 ierr = MatSetSizes(defl,m,PETSC_DECIDE,M,ceil(M/(float)len));CHKERRQ(ierr); 216 ierr = MatSetUp(defl);CHKERRQ(ierr); 217 ierr = MatSeqAIJSetPreallocation(defl,size,NULL);CHKERRQ(ierr); 218 ierr = MatMPIAIJSetPreallocation(defl,size,NULL,size,NULL);CHKERRQ(ierr); 219 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 220 221 ierr = MatGetOwnershipRangeColumn(defl,&ilo,&ihi);CHKERRQ(ierr); 222 for (i=0; i<len; i++) Iidx[i] = i+ilo*len; 223 if (M%len && ihi == (int)ceil(M/(float)len)) ihi -= 1; 224 for (i=ilo; i<ihi; i++) { 225 ierr = MatSetValues(defl,len,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr); 226 for (j=0; j<len; j++) Iidx[j] += len; 227 } 228 if (M%len && ihi+1 == ceil(M/(float)len)) { 229 len = M%len; 230 val = 1./pow(pow(2,len),0.5); 231 for (i=0; i<len; i++) col[i] = val; 232 ierr = MatSetValues(defl,len,Iidx,1,&ihi,col,INSERT_VALUES);CHKERRQ(ierr); 233 } 234 235 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 236 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 237 238 ierr = PetscFree(col);CHKERRQ(ierr); 239 ierr = PetscFree(Iidx);CHKERRQ(ierr); 240 *W = defl; 241 PetscFunctionReturn(0); 242 } 243 244 PetscErrorCode PCDeflationGetSpaceJacketHaar(PC pc,Mat *W,PetscInt size) 245 { 246 PetscErrorCode ierr; 247 Mat A,*H,defl; 248 PetscInt i,m,M,Mdefl,Ndefl; 249 PetscBool jh; 250 MPI_Comm comm; 251 252 PetscFunctionBegin; 253 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 254 ierr = PetscMalloc1(size,&H);CHKERRQ(ierr); 255 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 256 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 257 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 258 Mdefl = M; 259 Ndefl = M; 260 for (i=0; i<size; i++) { 261 if (Mdefl%2) { 262 jh=PETSC_TRUE; 263 Mdefl = Mdefl/2 +1; 264 } else { 265 jh=PETSC_FALSE; 266 Mdefl = Mdefl/2; 267 } 268 ierr = PCDeflationCreateSpaceJacketHaar(comm,PETSC_DECIDE,m,Mdefl,Ndefl,jh,&H[i]);CHKERRQ(ierr); 269 ierr = MatGetLocalSize(H[i],&m,NULL);CHKERRQ(ierr); 270 Ndefl = Mdefl; 271 } 272 //ierr = MatCreateProd(comm,size,H,&defl);CHKERRQ(ierr); 273 //ierr = MatCreateComposite(comm,size,H,&defl);CHKERRQ(ierr); 274 //ierr = MatCompositeSetType(defl,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 275 /* TODO allow implicit */ 276 //ierr = MatCompositeMerge(defl);CHKERRQ(ierr); 277 Mat newmat; 278 defl = H[0]; 279 for (i=0; i<size-1; i++) { 280 ierr = MatMatMult(H[i+1],defl,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 281 ierr = MatDestroy(&defl);CHKERRQ(ierr); 282 defl = newmat ; 283 } 284 285 ierr = MatTranspose(defl,MAT_INITIAL_MATRIX,W);CHKERRQ(ierr); 286 287 ierr = MatDestroy(&defl);CHKERRQ(ierr); 288 for (i=1; i<size; i++) { 289 ierr = MatDestroy(&H[i]);CHKERRQ(ierr); 290 } 291 ierr = PetscFree(H);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 PetscErrorCode PCDeflationGetSpaceWave(PC pc,Mat *W,PetscInt size,PetscInt ncoeffs,PetscScalar *coeffs,PetscBool trunc) 296 { 297 PetscErrorCode ierr; 298 Mat A,*H,defl; 299 PetscInt i,m,M,Mdefl,Ndefl; 300 MPI_Comm comm; 301 302 PetscFunctionBegin; 303 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 304 ierr = PetscMalloc1(size,&H);CHKERRQ(ierr); 305 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 306 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 307 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 308 Mdefl = M; 309 Ndefl = M; 310 for (i=0; i<size; i++) { 311 if (Mdefl%2) { 312 if (trunc) { 313 Mdefl = (PetscInt)PetscCeilReal(Mdefl/2.); 314 } else { 315 Mdefl = (PetscInt)PetscFloorReal((ncoeffs+Mdefl-1)/2.); 316 } 317 } else { 318 Mdefl = Mdefl/2; 319 } 320 ierr = PCDeflationCreateSpaceWave(comm,PETSC_DECIDE,m,Mdefl,Ndefl,ncoeffs,coeffs,trunc,&H[i]);CHKERRQ(ierr); 321 ierr = MatGetLocalSize(H[i],&m,NULL);CHKERRQ(ierr); 322 Ndefl = Mdefl; 323 } 324 ierr = MatCreateComposite(comm,size,H,&defl);CHKERRQ(ierr); 325 ierr = MatCompositeSetType(defl,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 326 *W = defl; 327 328 for (i=0; i<size; i++) { 329 ierr = MatDestroy(&H[i]);CHKERRQ(ierr); 330 } 331 ierr = PetscFree(H);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 PetscErrorCode PCDeflationGetSpaceAggregation(PC pc,Mat *W) 336 { 337 PetscErrorCode ierr; 338 Mat A,defl; 339 PetscInt i,ilo,ihi,*Iidx,m,M; 340 PetscReal *col; 341 MPI_Comm comm; 342 343 PetscFunctionBegin; 344 /* TODO pass A instead of PC? */ 345 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 346 ierr = MatGetOwnershipRangeColumn(A,&ilo,&ihi);CHKERRQ(ierr); 347 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 348 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 349 ierr = MPI_Comm_size(comm,&m);CHKERRQ(ierr); 350 ierr = MatCreate(comm,&defl);CHKERRQ(ierr); 351 ierr = MatSetSizes(defl,ihi-ilo,1,M,m);CHKERRQ(ierr); 352 ierr = MatSetUp(defl);CHKERRQ(ierr); 353 ierr = MatSeqAIJSetPreallocation(defl,1,NULL);CHKERRQ(ierr); 354 ierr = MatMPIAIJSetPreallocation(defl,1,NULL,0,NULL);CHKERRQ(ierr); 355 ierr = MatSetOption(defl,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 356 ierr = MatSetOption(defl,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 357 358 359 ierr = PetscMalloc1(ihi-ilo,&col);CHKERRQ(ierr); 360 ierr = PetscMalloc1(ihi-ilo,&Iidx);CHKERRQ(ierr); 361 for (i=ilo; i<ihi; i++) { 362 Iidx[i-ilo] = i; 363 col[i-ilo] = 1; 364 } 365 ierr = MPI_Comm_rank(comm,&i);CHKERRQ(ierr); 366 ierr = MatSetValues(defl,ihi-ilo,Iidx,1,&i,col,INSERT_VALUES);CHKERRQ(ierr); 367 368 ierr = MatAssemblyBegin(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 369 ierr = MatAssemblyEnd(defl,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 370 371 ierr = PetscFree(col);CHKERRQ(ierr); 372 ierr = PetscFree(Iidx);CHKERRQ(ierr); 373 *W = defl; 374 PetscFunctionReturn(0); 375 } 376 377 PetscErrorCode PCDeflationGetSpaceSLEPc(PC pc,Mat *W,PetscInt size,PetscBool cheapCP) 378 { 379 #if defined(HAVE_SLEPC) 380 PetscErrorCode ierr; 381 PC_Deflation *def = (PC_Deflation*)pc->data; 382 Mat A,defl; 383 Vec vec; 384 EPS eps; 385 PetscScalar *data,*dataScaled,eigval; 386 PetscInt i,nconv,m,M,n=PETSC_DECIDE; 387 PetscBool slepcinit; 388 MPI_Comm comm; 389 390 PetscFunctionBegin; 391 ierr = SlepcInitialized(&slepcinit);CHKERRQ(ierr); 392 if (!slepcinit) { 393 ierr = SlepcInitialize(NULL,NULL,(char*)0,(char*)0);CHKERRQ(ierr); 394 slepcinit = PETSC_TRUE; 395 } 396 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr); /* NOTE: Get Pmat instead? */ 397 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 398 ierr = EPSCreate(comm,&eps);CHKERRQ(ierr); 399 ierr = EPSSetOperators(eps,A,NULL);CHKERRQ(ierr); 400 ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr); /* Implemented only for def */ 401 ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr); 402 ierr = EPSSetDimensions(eps,size,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 403 ierr = EPSSetFromOptions(eps);CHKERRQ(ierr); 404 405 ierr = EPSSolve(eps);CHKERRQ(ierr); 406 ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr); 407 if (nconv < size) SETERRQ2(comm,PETSC_ERR_CONV_FAILED,"SLEPc: Number of converged eigenpairs (%d) is less than requested (%d)",nconv,size); 408 ierr = MatCreateVecs(A,NULL,&vec);CHKERRQ(ierr); 409 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr); 410 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 411 ierr = PetscSplitOwnership(comm,&n,&size);CHKERRQ(ierr); 412 ierr = PetscMalloc1(m*size,&data);CHKERRQ(ierr); 413 /* TODO check that eigenvalue is not 0 -> vec is not in Ker A */ 414 for (i=0; i<size; i++) { 415 ierr = VecPlaceArray(vec,&data[i*m]);CHKERRQ(ierr); 416 ierr = EPSGetEigenvector(eps,i,vec,NULL);CHKERRQ(ierr); 417 ierr = VecResetArray(vec);CHKERRQ(ierr); 418 } 419 ierr = MatCreateDense(comm,m,n,M,size,data,&defl);CHKERRQ(ierr); 420 421 if (cheapCP) { 422 ierr = PetscMalloc1(m*size,&dataScaled);CHKERRQ(ierr); 423 for (i=0; i<size; i++) { 424 ierr = VecPlaceArray(vec,&dataScaled[i*m]);CHKERRQ(ierr); 425 ierr = EPSGetEigenpair(eps,i,&eigval,NULL,vec,NULL);CHKERRQ(ierr); 426 ierr = VecScale(vec,eigval);CHKERRQ(ierr); 427 ierr = VecResetArray(vec);CHKERRQ(ierr); 428 } 429 ierr = MatCreateDense(comm,m,n,M,size,dataScaled,&def->AW);CHKERRQ(ierr); 430 } 431 432 //ierr = EPSGetBV(eps,&bv);CHKERRQ(ierr); 433 //ierr = BVCreateMat(bv,&defl);CHKERRQ(ierr); 434 *W = defl; 435 436 ierr = EPSDestroy(&eps);CHKERRQ(ierr); 437 if (slepcinit) ierr = SlepcFinalize();CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 #else 440 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_CONV_FAILED,"Not compiled with SLEPc support (call make HAVE_SLEPC)"); 441 #endif 442 } 443 444 PetscErrorCode PCDeflationComputeSpace(PC pc) 445 { 446 PetscErrorCode ierr; 447 Mat defl; 448 PetscBool transp=PETSC_TRUE; 449 PC_Deflation *def = (PC_Deflation*)pc->data; 450 451 /* TODO valid header */ 452 PetscFunctionBegin; 453 if (def->spacesize < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space size specified: %d",def->spacesize); 454 switch (def->spacetype) { 455 case PC_DEFLATION_SPACE_HAAR: 456 transp = PETSC_FALSE; 457 ierr = PCDeflationGetSpaceHaar(pc,&defl,def->spacesize);CHKERRQ(ierr);break; 458 case PC_DEFLATION_SPACE_JACKET_HAAR: 459 transp = PETSC_FALSE; 460 ierr = PCDeflationGetSpaceJacketHaar(pc,&defl,def->spacesize);CHKERRQ(ierr);break; 461 case PC_DEFLATION_SPACE_DB2: 462 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,2,db2,!def->extendsp);CHKERRQ(ierr);break; 463 case PC_DEFLATION_SPACE_DB4: 464 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,4,db4,!def->extendsp);CHKERRQ(ierr);break; 465 case PC_DEFLATION_SPACE_DB8: 466 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,8,db8,!def->extendsp);CHKERRQ(ierr);break; 467 case PC_DEFLATION_SPACE_DB16: 468 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,16,db16,!def->extendsp);CHKERRQ(ierr);break; 469 case PC_DEFLATION_SPACE_BIORTH22: 470 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,6,biorth22,!def->extendsp);CHKERRQ(ierr);break; 471 case PC_DEFLATION_SPACE_MEYER: 472 ierr = PCDeflationGetSpaceWave(pc,&defl,def->spacesize,62,meyer,!def->extendsp);CHKERRQ(ierr);break; 473 case PC_DEFLATION_SPACE_AGGREGATION: 474 transp = PETSC_FALSE; 475 ierr = PCDeflationGetSpaceAggregation(pc,&defl);CHKERRQ(ierr);break; 476 case PC_DEFLATION_SPACE_SLEPC: 477 transp = PETSC_FALSE; 478 ierr = PCDeflationGetSpaceSLEPc(pc,&defl,def->spacesize,PETSC_FALSE);CHKERRQ(ierr);break; 479 case PC_DEFLATION_SPACE_SLEPC_CHEAP: 480 transp = PETSC_FALSE; 481 ierr = PCDeflationGetSpaceSLEPc(pc,&defl,def->spacesize,PETSC_TRUE);CHKERRQ(ierr);break; 482 default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Wrong PC_DEFLATION Space Type specified"); 483 } 484 485 ierr = PCDeflationSetSpace(pc,defl,transp);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 489