Actual source code: matusfft.c
petsc-3.14.6 2021-03-30
2: /*
3: Provides an implementation of the Unevenly Sampled FFT algorithm as a Mat.
4: Testing examples can be found in ~/src/mat/tests FIX: should these be moved to dm/da/tests?
5: */
7: #include <petsc/private/matimpl.h>
8: #include <petscdmda.h>
9: #include <fftw3.h>
11: typedef struct {
12: PetscInt dim;
13: Vec sampleCoords;
14: PetscInt dof;
15: DM freqDA; /* frequency DMDA */
16: PetscInt *freqSizes; /* sizes of the frequency DMDA, one per each dim */
17: DM resampleDa; /* the Battle-Lemarie interpolant DMDA */
18: Vec resample; /* Vec of samples, one per dof per sample point */
19: fftw_plan p_forward,p_backward;
20: unsigned p_flag; /* planner flags, FFTW_ESTIMATE,FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE */
21: } Mat_USFFT;
24: PetscErrorCode MatApply_USFFT_Private(Mat A, fftw_plan *plan, int direction, Vec x,Vec y)
25: {
26: #if 0
28: PetscScalar *r_array, *y_array;
29: Mat_USFFT* = (Mat_USFFT*)(A->data);
30: #endif
33: #if 0
34: /* resample x to usfft->resample */
35: MatResample_USFFT_Private(A, x);
37: /* NB: for now we use outdim for both x and y; this will change once a full USFFT is implemented */
38: VecGetArray(usfft->resample,&r_array);
39: VecGetArray(y,&y_array);
40: if (!*plan) { /* create a plan then execute it*/
41: if (usfft->dof == 1) {
42: #if defined(PETSC_DEBUG_USFFT)
43: PetscPrintf(PetscObjectComm((PetscObject)A), "direction = %d, usfft->ndim = %d\n", direction, usfft->ndim);
44: for (int ii = 0; ii < usfft->ndim; ++ii) {
45: PetscPrintf(PetscObjectComm((PetscObject)A), "usfft->outdim[%d] = %d\n", ii, usfft->outdim[ii]);
46: }
47: #endif
49: switch (usfft->dim) {
50: case 1:
51: *plan = fftw_plan_dft_1d(usfft->outdim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,direction,usfft->p_flag);
52: break;
53: case 2:
54: *plan = fftw_plan_dft_2d(usfft->outdim[0],usfft->outdim[1],(fftw_complex*)x_array,(fftw_complex*)y_array,direction,usfft->p_flag);
55: break;
56: case 3:
57: *plan = fftw_plan_dft_3d(usfft->outdim[0],usfft->outdim[1],usfft->outdim[2],(fftw_complex*)x_array,(fftw_complex*)y_array,direction,usfft->p_flag);
58: break;
59: default:
60: *plan = fftw_plan_dft(usfft->ndim,usfft->outdim,(fftw_complex*)x_array,(fftw_complex*)y_array,direction,usfft->p_flag);
61: break;
62: }
63: fftw_execute(*plan);
64: } /* if (dof == 1) */
65: else { /* if (dof > 1) */
66: *plan = fftw_plan_many_dft(/*rank*/usfft->ndim, /*n*/usfft->outdim, /*howmany*/usfft->dof,
67: (fftw_complex*)x_array, /*nembed*/usfft->outdim, /*stride*/usfft->dof, /*dist*/1,
68: (fftw_complex*)y_array, /*nembed*/usfft->outdim, /*stride*/usfft->dof, /*dist*/1,
69: /*sign*/direction, /*flags*/usfft->p_flag);
70: fftw_execute(*plan);
71: } /* if (dof > 1) */
72: } /* if (!*plan) */
73: else { /* if (*plan) */
74: /* use existing plan */
75: fftw_execute_dft(*plan,(fftw_complex*)x_array,(fftw_complex*)y_array);
76: }
77: VecRestoreArray(y,&y_array);
78: VecRestoreArray(x,&x_array);
79: #endif
80: return(0);
81: } /* MatApply_USFFT_Private() */
83: #if 0
84: PetscErrorCode MatUSFFT_ProjectOnBattleLemarie_Private(Vec x,double *r)
85: /* Project onto the Battle-Lemarie function centered around r */
86: {
88: PetscScalar *x_array, *y_array;
91: return(0);
92: } /* MatUSFFT_ProjectOnBattleLemarie_Private() */
94: PetscErrorCode MatInterpolate_USFFT_Private(Vec x,Vec y)
95: {
97: PetscScalar *x_array, *y_array;
100: return(0);
101: } /* MatInterpolate_USFFT_Private() */
104: PetscErrorCode MatMult_SeqUSFFT(Mat A,Vec x,Vec y)
105: {
107: Mat_USFFT *usfft = (Mat_USFFT*)A->data;
110: /* NB: for now we use outdim for both x and y; this will change once a full USFFT is implemented */
111: MatApply_USFFT_Private(A, &usfft->p_forward, FFTW_FORWARD, x,y);
112: return(0);
113: }
115: PetscErrorCode MatMultTranspose_SeqUSFFT(Mat A,Vec x,Vec y)
116: {
118: Mat_USFFT *usfft = (Mat_USFFT*)A->data;
121: /* NB: for now we use outdim for both x and y; this will change once a full USFFT is implemented */
122: MatApply_USFFT_Private(usfft, &usfft->p_backward, FFTW_BACKWARD, x,y);
123: return(0);
124: }
126: PetscErrorCode MatDestroy_SeqUSFFT(Mat A)
127: {
128: Mat_USFFT *usfft = (Mat_USFFT*)A->data;
132: fftw_destroy_plan(usfft->p_forward);
133: fftw_destroy_plan(usfft->p_backward);
134: PetscFree(usfft->indim);
135: PetscFree(usfft->outdim);
136: PetscFree(usfft);
137: PetscObjectChangeTypeName((PetscObject)A,0);
138: return(0);
139: } /* MatDestroy_SeqUSFFT() */
142: /*@C
143: MatCreateSeqUSFFT - Creates a matrix object that provides sequential USFFT
144: via the external package FFTW
146: Collective
148: Input Parameter:
149: . da - geometry of the domain encoded by a DMDA
151: Output Parameter:
152: . A - the matrix
154: Options Database Keys:
155: . -mat_usfft_plannerflags - set the FFTW planner flags
157: Level: intermediate
159: @*/
160: PetscErrorCode MatCreateSeqUSFFT(Vec sampleCoords, DMDA freqDA, Mat *A)
161: {
163: Mat_USFFT *usfft;
164: PetscInt m,n,M,N,i;
165: const char *p_flags[]={"FFTW_ESTIMATE","FFTW_MEASURE","FFTW_PATIENT","FFTW_EXHAUSTIVE"};
166: PetscBool flg;
167: PetscInt p_flag;
168: PetscInt dof, dim, freqSizes[3];
169: MPI_Comm comm;
170: PetscInt size;
173: PetscObjectGetComm((PetscObject)inda, &comm);
174: MPI_Comm_size(comm, &size);
175: if (size > 1) SETERRQ(comm,PETSC_ERR_USER, "Parallel DMDA (in) not yet supported by USFFT");
176: PetscObjectGetComm((PetscObject)outda, &comm);
177: MPI_Comm_size(comm, &size);
178: if (size > 1) SETERRQ(comm,PETSC_ERR_USER, "Parallel DMDA (out) not yet supported by USFFT");
179: MatCreate(comm,A);
180: PetscNewLog(*A,&usfft);
181: (*A)->data = (void*)usfft;
182: usfft->inda = inda;
183: usfft->outda = outda;
184: /* inda */
185: DMDAGetInfo(usfft->inda, &ndim, dim+0, dim+1, dim+2, NULL, NULL, NULL, &dof, NULL, NULL, NULL);
186: if (ndim <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"ndim %d must be > 0",ndim);
187: if (dof <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"dof %d must be > 0",dof);
188: usfft->ndim = ndim;
189: usfft->dof = dof;
190: usfft->freqDA = freqDA;
191: /* NB: we reverse the freq and resample DMDA sizes, since the DMDA ordering (natural on x-y-z, with x varying the fastest)
192: is the order opposite of that assumed by FFTW: z varying the fastest */
193: PetscMalloc1(usfft->ndim+1,&usfft->indim);
194: for (i = usfft->ndim; i > 0; --i) usfft->indim[usfft->ndim-i] = dim[i-1];
196: /* outda */
197: DMDAGetInfo(usfft->outda, &ndim, dim+0, dim+1, dim+2, NULL, NULL, NULL, &dof, NULL, NULL, NULL);
198: if (ndim != usfft->ndim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"in and out DMDA dimensions must match: %d != %d",usfft->ndim, ndim);
199: if (dof != usfft->dof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"in and out DMDA dof must match: %d != %d",usfft->dof, dof);
200: /* Store output dimensions */
201: /* NB: we reverse the DMDA dimensions, since the DMDA ordering (natural on x-y-z, with x varying the fastest)
202: is the order opposite of that assumed by FFTW: z varying the fastest */
203: PetscMalloc1(usfft->ndim+1,&usfft->outdim);
204: for (i = usfft->ndim; i > 0; --i) usfft->outdim[usfft->ndim-i] = dim[i-1];
206: /* TODO: Use the new form of DMDACreate() */
207: #if 0
208: DMDACreate(comm,usfft->dim, DMDA_NONPERIODIC, DMDA_STENCIL_STAR, usfft->freqSizes[0], usfft->freqSizes[1], usfft->freqSizes[2],
209: PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, 0, NULL, NULL, NULL, 0, &(usfft->resampleDA));
210: #endif
211: DMDAGetVec(usfft->resampleDA, usfft->resample);
214: /* CONTINUE: Need to build the connectivity "Sieve" attaching sample points to the resample points they are close to */
216: /* CONTINUE: recalculate matrix sizes based on the connectivity "Sieve" */
217: /* mat sizes */
218: m = 1; n = 1;
219: for (i=0; i<usfft->ndim; i++) {
220: if (usfft->indim[i] <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"indim[%d]=%d must be > 0",i,usfft->indim[i]);
221: if (usfft->outdim[i] <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"outdim[%d]=%d must be > 0",i,usfft->outdim[i]);
222: n *= usfft->indim[i];
223: m *= usfft->outdim[i];
224: }
225: N = n*usfft->dof;
226: M = m*usfft->dof;
227: MatSetSizes(*A,M,N,M,N); /* "in size" is the number of columns, "out size" is the number of rows" */
228: PetscObjectChangeTypeName((PetscObject)*A,MATSEQUSFFT);
229: usfft->m = m; usfft->n = n; usfft->M = M; usfft->N = N;
230: /* FFTW */
231: usfft->p_forward = 0;
232: usfft->p_backward = 0;
233: usfft->p_flag = FFTW_ESTIMATE;
234: /* set Mat ops */
235: (*A)->ops->mult = MatMult_SeqUSFFT;
236: (*A)->ops->multtranspose = MatMultTranspose_SeqUSFFT;
237: (*A)->assembled = PETSC_TRUE;
238: (*A)->ops->destroy = MatDestroy_SeqUSFFT;
239: /* get runtime options */
240: PetscOptionsBegin(((PetscObject)(*A))->comm,((PetscObject)(*A))->prefix,"USFFT Options","Mat");
241: PetscOptionsEList("-mat_usfft_fftw_plannerflags","Planner Flags","None",p_flags,4,p_flags[0],&p_flag,&flg);
242: if (flg) usfft->p_flag = (unsigned)p_flag;
243: PetscOptionsEnd();
244: return(0);
245: } /* MatCreateSeqUSFFT() */
247: #endif