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