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