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, 0.032883011666982945, 0.030841381835986965, -0.18703481171888114, -0.02798376941698385, 0.6308807679295904, 0.7148465705525415, 0.23037781330885523}; 8 9 PetscScalar db16[] = {-0.00011747678400228192, 0.0006754494059985568, -0.0003917403729959771, -0.00487035299301066, 0.008746094047015655, 0.013981027917015516, -0.04408825393106472, -0.01736930100202211, 10 0.128747426620186, 0.00047248457399797254, -0.2840155429624281, -0.015829105256023893, 0.5853546836548691, 0.6756307362980128, 0.3128715909144659, 0.05441584224308161}; 11 12 PetscScalar biorth22[] = {0.0, -0.1767766952966369, 0.3535533905932738, 1.0606601717798214, 0.3535533905932738, -0.1767766952966369}; 13 14 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}; 15 16 static PetscErrorCode PCDeflationCreateSpaceWave(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt ncoeffs, PetscScalar *coeffs, PetscBool trunc, Mat *H) { 17 Mat defl; 18 PetscInt i, j, k, ilo, ihi, *Iidx; 19 20 PetscFunctionBegin; 21 PetscCall(PetscMalloc1(ncoeffs, &Iidx)); 22 23 PetscCall(MatCreate(comm, &defl)); 24 PetscCall(MatSetSizes(defl, m, n, M, N)); 25 PetscCall(MatSetUp(defl)); 26 PetscCall(MatSeqAIJSetPreallocation(defl, ncoeffs, NULL)); 27 PetscCall(MatMPIAIJSetPreallocation(defl, ncoeffs, NULL, ncoeffs, NULL)); 28 PetscCall(MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 29 PetscCall(MatSetOption(defl, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 30 31 /* Alg 735 Taswell: fvecmat */ 32 k = ncoeffs - 2; 33 if (trunc) k = k / 2; 34 35 PetscCall(MatGetOwnershipRange(defl, &ilo, &ihi)); 36 for (i = 0; i < ncoeffs; i++) { 37 Iidx[i] = i + ilo * 2 - k; 38 if (Iidx[i] >= N) Iidx[i] = PETSC_MIN_INT; 39 } 40 for (i = ilo; i < ihi; i++) { 41 PetscCall(MatSetValues(defl, 1, &i, ncoeffs, Iidx, coeffs, INSERT_VALUES)); 42 for (j = 0; j < ncoeffs; j++) { 43 Iidx[j] += 2; 44 if (Iidx[j] >= N) Iidx[j] = PETSC_MIN_INT; 45 } 46 } 47 48 PetscCall(MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY)); 49 PetscCall(MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY)); 50 51 PetscCall(PetscFree(Iidx)); 52 *H = defl; 53 PetscFunctionReturn(0); 54 } 55 56 PetscErrorCode PCDeflationGetSpaceHaar(PC pc, Mat *W, PetscInt size) { 57 Mat A, defl; 58 PetscInt i, j, len, ilo, ihi, *Iidx, m, M; 59 PetscScalar *col, val; 60 61 PetscFunctionBegin; 62 /* Haar basis wavelet, level=size */ 63 len = pow(2, size); 64 PetscCall(PetscMalloc2(len, &col, len, &Iidx)); 65 val = 1. / pow(2, size / 2.); 66 for (i = 0; i < len; i++) col[i] = val; 67 68 PetscCall(PCGetOperators(pc, NULL, &A)); 69 PetscCall(MatGetLocalSize(A, &m, NULL)); 70 PetscCall(MatGetSize(A, &M, NULL)); 71 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &defl)); 72 PetscCall(MatSetSizes(defl, m, PETSC_DECIDE, M, PetscCeilInt(M, len))); 73 PetscCall(MatSetUp(defl)); 74 PetscCall(MatSeqAIJSetPreallocation(defl, size, NULL)); 75 PetscCall(MatMPIAIJSetPreallocation(defl, size, NULL, size, NULL)); 76 PetscCall(MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 77 78 PetscCall(MatGetOwnershipRangeColumn(defl, &ilo, &ihi)); 79 for (i = 0; i < len; i++) Iidx[i] = i + ilo * len; 80 if (M % len && ihi == PetscCeilInt(M, len)) ihi -= 1; 81 for (i = ilo; i < ihi; i++) { 82 PetscCall(MatSetValues(defl, len, Iidx, 1, &i, col, INSERT_VALUES)); 83 for (j = 0; j < len; j++) Iidx[j] += len; 84 } 85 if (M % len && ihi + 1 == PetscCeilInt(M, len)) { 86 len = M % len; 87 val = 1. / pow(pow(2, len), 0.5); 88 for (i = 0; i < len; i++) col[i] = val; 89 PetscCall(MatSetValues(defl, len, Iidx, 1, &ihi, col, INSERT_VALUES)); 90 } 91 92 PetscCall(MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY)); 93 PetscCall(MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY)); 94 95 PetscCall(PetscFree2(col, Iidx)); 96 *W = defl; 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode PCDeflationGetSpaceWave(PC pc, Mat *W, PetscInt size, PetscInt ncoeffs, PetscScalar *coeffs, PetscBool trunc) { 101 Mat A, *H, defl; 102 PetscInt i, m, M, Mdefl, Ndefl; 103 MPI_Comm comm; 104 105 PetscFunctionBegin; 106 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 107 PetscCall(PetscMalloc1(size, &H)); 108 PetscCall(PCGetOperators(pc, &A, NULL)); 109 PetscCall(MatGetLocalSize(A, &m, NULL)); 110 PetscCall(MatGetSize(A, &M, NULL)); 111 Mdefl = M; 112 Ndefl = M; 113 for (i = 0; i < size; i++) { 114 if (Mdefl % 2) { 115 if (trunc) Mdefl = (PetscInt)PetscCeilReal(Mdefl / 2.); 116 else Mdefl = (PetscInt)PetscFloorReal((ncoeffs + Mdefl - 1) / 2.); 117 } else Mdefl = Mdefl / 2; 118 PetscCall(PCDeflationCreateSpaceWave(comm, PETSC_DECIDE, m, Mdefl, Ndefl, ncoeffs, coeffs, trunc, &H[i])); 119 PetscCall(MatGetLocalSize(H[i], &m, NULL)); 120 Ndefl = Mdefl; 121 } 122 PetscCall(MatCreateComposite(comm, size, H, &defl)); 123 PetscCall(MatCompositeSetType(defl, MAT_COMPOSITE_MULTIPLICATIVE)); 124 *W = defl; 125 126 for (i = 0; i < size; i++) { PetscCall(MatDestroy(&H[i])); } 127 PetscCall(PetscFree(H)); 128 PetscFunctionReturn(0); 129 } 130 131 PetscErrorCode PCDeflationGetSpaceAggregation(PC pc, Mat *W) { 132 Mat A, defl; 133 PetscInt i, ilo, ihi, *Iidx, M; 134 PetscMPIInt m; 135 PetscScalar *col; 136 MPI_Comm comm; 137 138 PetscFunctionBegin; 139 PetscCall(PCGetOperators(pc, &A, NULL)); 140 PetscCall(MatGetOwnershipRangeColumn(A, &ilo, &ihi)); 141 PetscCall(MatGetSize(A, &M, NULL)); 142 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 143 PetscCallMPI(MPI_Comm_size(comm, &m)); 144 PetscCall(MatCreate(comm, &defl)); 145 PetscCall(MatSetSizes(defl, ihi - ilo, 1, M, m)); 146 PetscCall(MatSetUp(defl)); 147 PetscCall(MatSeqAIJSetPreallocation(defl, 1, NULL)); 148 PetscCall(MatMPIAIJSetPreallocation(defl, 1, NULL, 0, NULL)); 149 PetscCall(MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 150 PetscCall(MatSetOption(defl, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 151 152 PetscCall(PetscMalloc2(ihi - ilo, &col, ihi - ilo, &Iidx)); 153 for (i = ilo; i < ihi; i++) { 154 Iidx[i - ilo] = i; 155 col[i - ilo] = 1; 156 } 157 PetscCallMPI(MPI_Comm_rank(comm, &m)); 158 i = m; 159 PetscCall(MatSetValues(defl, ihi - ilo, Iidx, 1, &i, col, INSERT_VALUES)); 160 161 PetscCall(MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY)); 162 PetscCall(MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY)); 163 164 PetscCall(PetscFree2(col, Iidx)); 165 *W = defl; 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode PCDeflationComputeSpace(PC pc) { 170 Mat defl; 171 PetscBool transp = PETSC_TRUE; 172 PC_Deflation *def = (PC_Deflation *)pc->data; 173 174 PetscFunctionBegin; 175 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 176 PetscCheck(def->spacesize >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Wrong PCDeflation space size specified: %" PetscInt_FMT, def->spacesize); 177 switch (def->spacetype) { 178 case PC_DEFLATION_SPACE_HAAR: 179 transp = PETSC_FALSE; 180 PetscCall(PCDeflationGetSpaceHaar(pc, &defl, def->spacesize)); 181 break; 182 case PC_DEFLATION_SPACE_DB2: PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 2, db2, PetscNot(def->extendsp))); break; 183 case PC_DEFLATION_SPACE_DB4: PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 4, db4, PetscNot(def->extendsp))); break; 184 case PC_DEFLATION_SPACE_DB8: PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 8, db8, PetscNot(def->extendsp))); break; 185 case PC_DEFLATION_SPACE_DB16: PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 16, db16, PetscNot(def->extendsp))); break; 186 case PC_DEFLATION_SPACE_BIORTH22: PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 6, biorth22, PetscNot(def->extendsp))); break; 187 case PC_DEFLATION_SPACE_MEYER: PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 62, meyer, PetscNot(def->extendsp))); break; 188 case PC_DEFLATION_SPACE_AGGREGATION: 189 transp = PETSC_FALSE; 190 PetscCall(PCDeflationGetSpaceAggregation(pc, &defl)); 191 break; 192 default: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Wrong PCDeflationSpaceType specified"); 193 } 194 195 PetscCall(PCDeflationSetSpace(pc, defl, transp)); 196 PetscCall(MatDestroy(&defl)); 197 PetscFunctionReturn(0); 198 } 199