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