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