xref: /petsc/src/mat/impls/fft/fftw/fftw.c (revision d63f877fec171bed822722979cba2adffb37de2a) 
1 
2 /*
3     Provides an interface to the FFTW package.
4     Testing examples can be found in ~src/mat/examples/tests
5 */
6 
7 #include <../src/mat/impls/fft/fft.h>   /*I "petscmat.h" I*/
8 EXTERN_C_BEGIN
9 #include <fftw3-mpi.h>
10 EXTERN_C_END
11 
12 typedef struct {
13   fftw_plan   p_forward,p_backward;
14   unsigned    p_flag; /* planner flags, FFTW_ESTIMATE,FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE */
15   PetscScalar *finarray,*foutarray,*binarray,*boutarray; /* keep track of arrays becaue fftw plan should be
16                                                             executed for the arrays with which the plan was created */
17 } Mat_FFTW;
18 
19 extern PetscErrorCode MatMult_SeqFFTW(Mat,Vec,Vec);
20 extern PetscErrorCode MatMultTranspose_SeqFFTW(Mat,Vec,Vec);
21 extern PetscErrorCode MatMult_MPIFFTW(Mat,Vec,Vec);
22 extern PetscErrorCode MatMultTranspose_MPIFFTW(Mat,Vec,Vec);
23 extern PetscErrorCode MatDestroy_FFTW(Mat);
24 extern PetscErrorCode VecDestroy_MPIFFTW(Vec);
25 extern PetscErrorCode MatGetVecs_FFTW(Mat,Vec*,Vec*);
26 
27 #undef __FUNCT__
28 #define __FUNCT__ "MatMult_SeqFFTW"
29 PetscErrorCode MatMult_SeqFFTW(Mat A,Vec x,Vec y)
30 {
31   PetscErrorCode ierr;
32   Mat_FFT        *fft  = (Mat_FFT*)A->data;
33   Mat_FFTW       *fftw = (Mat_FFTW*)fft->data;
34   PetscScalar    *x_array,*y_array;
35   PetscInt       ndim=fft->ndim,*dim=fft->dim;
36 
37   PetscFunctionBegin;
38 #if !defined(PETSC_USE_COMPLEX)
39   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers");
40 #endif
41   ierr = VecGetArray(x,&x_array);CHKERRQ(ierr);
42   ierr = VecGetArray(y,&y_array);CHKERRQ(ierr);
43   if (!fftw->p_forward){ /* create a plan, then excute it */
44     switch (ndim){
45     case 1:
46       fftw->p_forward = fftw_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_FORWARD,fftw->p_flag);
47       break;
48     case 2:
49       fftw->p_forward = fftw_plan_dft_2d(dim[0],dim[1],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_FORWARD,fftw->p_flag);
50       break;
51     case 3:
52       fftw->p_forward = fftw_plan_dft_3d(dim[0],dim[1],dim[2],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_FORWARD,fftw->p_flag);
53       break;
54     default:
55       fftw->p_forward = fftw_plan_dft(ndim,dim,(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_FORWARD,fftw->p_flag);
56       break;
57     }
58     fftw->finarray  = x_array;
59     fftw->foutarray = y_array;
60     /* Warning: if (fftw->p_flag!==FFTW_ESTIMATE) The data in the in/out arrays is overwritten!
61                 planning should be done before x is initialized! See FFTW manual sec2.1 or sec4 */
62     fftw_execute(fftw->p_forward);
63   } else { /* use existing plan */
64     if (fftw->finarray != x_array || fftw->foutarray != y_array){ /* use existing plan on new arrays */
65       fftw_execute_dft(fftw->p_forward,(fftw_complex*)x_array,(fftw_complex*)y_array);
66     } else {
67       fftw_execute(fftw->p_forward);
68     }
69   }
70   ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr);
71   ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr);
72   PetscFunctionReturn(0);
73 }
74 
75 #undef __FUNCT__
76 #define __FUNCT__ "MatMultTranspose_SeqFFTW"
77 PetscErrorCode MatMultTranspose_SeqFFTW(Mat A,Vec x,Vec y)
78 {
79   PetscErrorCode ierr;
80   Mat_FFT        *fft = (Mat_FFT*)A->data;
81   Mat_FFTW       *fftw = (Mat_FFTW*)fft->data;
82   PetscScalar    *x_array,*y_array;
83   PetscInt       ndim=fft->ndim,*dim=fft->dim;
84 
85   PetscFunctionBegin;
86 #if !defined(PETSC_USE_COMPLEX)
87   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers");
88 #endif
89   ierr = VecGetArray(x,&x_array);CHKERRQ(ierr);
90   ierr = VecGetArray(y,&y_array);CHKERRQ(ierr);
91   if (!fftw->p_backward){ /* create a plan, then excute it */
92     switch (ndim){
93     case 1:
94       fftw->p_backward = fftw_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_BACKWARD,fftw->p_flag);
95       break;
96     case 2:
97       fftw->p_backward = fftw_plan_dft_2d(dim[0],dim[1],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_BACKWARD,fftw->p_flag);
98       break;
99     case 3:
100       fftw->p_backward = fftw_plan_dft_3d(dim[0],dim[1],dim[2],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_BACKWARD,fftw->p_flag);
101       break;
102     default:
103       fftw->p_backward = fftw_plan_dft(ndim,dim,(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_BACKWARD,fftw->p_flag);
104       break;
105     }
106     fftw->binarray  = x_array;
107     fftw->boutarray = y_array;
108     fftw_execute(fftw->p_backward);CHKERRQ(ierr);
109   } else { /* use existing plan */
110     if (fftw->binarray != x_array || fftw->boutarray != y_array){ /* use existing plan on new arrays */
111       fftw_execute_dft(fftw->p_backward,(fftw_complex*)x_array,(fftw_complex*)y_array);
112     } else {
113       fftw_execute(fftw->p_backward);CHKERRQ(ierr);
114     }
115   }
116   ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr);
117   ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr);
118   PetscFunctionReturn(0);
119 }
120 
121 #undef __FUNCT__
122 #define __FUNCT__ "MatMult_MPIFFTW"
123 PetscErrorCode MatMult_MPIFFTW(Mat A,Vec x,Vec y)
124 {
125   PetscErrorCode ierr;
126   Mat_FFT        *fft  = (Mat_FFT*)A->data;
127   Mat_FFTW       *fftw = (Mat_FFTW*)fft->data;
128   PetscScalar    *x_array,*y_array;
129   PetscInt       ndim=fft->ndim,*dim=fft->dim;
130   MPI_Comm       comm=((PetscObject)A)->comm;
131 
132   PetscFunctionBegin;
133 #if !defined(PETSC_USE_COMPLEX)
134   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers");
135 #endif
136   ierr = VecGetArray(x,&x_array);CHKERRQ(ierr);
137   ierr = VecGetArray(y,&y_array);CHKERRQ(ierr);
138   if (!fftw->p_forward){ /* create a plan, then excute it */
139     switch (ndim){
140     case 1:
141       fftw->p_forward = fftw_mpi_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_FORWARD,fftw->p_flag);
142       break;
143     case 2:
144       fftw->p_forward = fftw_mpi_plan_dft_2d(dim[0],dim[1],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_FORWARD,fftw->p_flag);
145       break;
146     case 3:
147       fftw->p_forward = fftw_mpi_plan_dft_3d(dim[0],dim[1],dim[2],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_FORWARD,fftw->p_flag);
148       break;
149     default:
150       fftw->p_forward = fftw_mpi_plan_dft(ndim,(const ptrdiff_t*)dim,(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_FORWARD,fftw->p_flag);
151       break;
152     }
153     fftw->finarray  = x_array;
154     fftw->foutarray = y_array;
155     /* Warning: if (fftw->p_flag!==FFTW_ESTIMATE) The data in the in/out arrays is overwritten!
156                 planning should be done before x is initialized! See FFTW manual sec2.1 or sec4 */
157     fftw_execute(fftw->p_forward);
158   } else { /* use existing plan */
159     if (fftw->finarray != x_array || fftw->foutarray != y_array){ /* use existing plan on new arrays */
160       fftw_execute_dft(fftw->p_forward,(fftw_complex*)x_array,(fftw_complex*)y_array);
161     } else {
162       fftw_execute(fftw->p_forward);
163     }
164   }
165   ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr);
166   ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr);
167   PetscFunctionReturn(0);
168 }
169 
170 #undef __FUNCT__
171 #define __FUNCT__ "MatMultTranspose_MPIFFTW"
172 PetscErrorCode MatMultTranspose_MPIFFTW(Mat A,Vec x,Vec y)
173 {
174   PetscErrorCode ierr;
175   Mat_FFT        *fft  = (Mat_FFT*)A->data;
176   Mat_FFTW       *fftw = (Mat_FFTW*)fft->data;
177   PetscScalar    *x_array,*y_array;
178   PetscInt       ndim=fft->ndim,*dim=fft->dim;
179   MPI_Comm       comm=((PetscObject)A)->comm;
180 
181   PetscFunctionBegin;
182 #if !defined(PETSC_USE_COMPLEX)
183   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers");
184 #endif
185   ierr = VecGetArray(x,&x_array);CHKERRQ(ierr);
186   ierr = VecGetArray(y,&y_array);CHKERRQ(ierr);
187   if (!fftw->p_backward){ /* create a plan, then excute it */
188     switch (ndim){
189     case 1:
190       fftw->p_backward = fftw_mpi_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_BACKWARD,fftw->p_flag);
191       break;
192     case 2:
193       fftw->p_backward = fftw_mpi_plan_dft_2d(dim[0],dim[1],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_BACKWARD,fftw->p_flag);
194       break;
195     case 3:
196       fftw->p_backward = fftw_mpi_plan_dft_3d(dim[0],dim[1],dim[2],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_BACKWARD,fftw->p_flag);
197       break;
198     default:
199       fftw->p_backward = fftw_mpi_plan_dft(ndim,(const ptrdiff_t*)dim,(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_BACKWARD,fftw->p_flag);
200       break;
201     }
202     fftw->binarray  = x_array;
203     fftw->boutarray = y_array;
204     fftw_execute(fftw->p_backward);CHKERRQ(ierr);
205   } else { /* use existing plan */
206     if (fftw->binarray != x_array || fftw->boutarray != y_array){ /* use existing plan on new arrays */
207       fftw_execute_dft(fftw->p_backward,(fftw_complex*)x_array,(fftw_complex*)y_array);
208     } else {
209       fftw_execute(fftw->p_backward);CHKERRQ(ierr);
210     }
211   }
212   ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr);
213   ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr);
214   PetscFunctionReturn(0);
215 }
216 
217 #undef __FUNCT__
218 #define __FUNCT__ "MatDestroy_FFTW"
219 PetscErrorCode MatDestroy_FFTW(Mat A)
220 {
221   Mat_FFT        *fft = (Mat_FFT*)A->data;
222   Mat_FFTW       *fftw = (Mat_FFTW*)fft->data;
223   PetscErrorCode ierr;
224 
225   PetscFunctionBegin;
226 #if !defined(PETSC_USE_COMPLEX)
227   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers");
228 #endif
229   fftw_destroy_plan(fftw->p_forward);
230   fftw_destroy_plan(fftw->p_backward);
231   ierr = PetscFree(fft->data);CHKERRQ(ierr);
232   PetscFunctionReturn(0);
233 }
234 
235 #include <../src/vec/vec/impls/mpi/pvecimpl.h>   /*I  "petscvec.h"   I*/
236 #undef __FUNCT__
237 #define __FUNCT__ "VecDestroy_MPIFFTW"
238 PetscErrorCode VecDestroy_MPIFFTW(Vec v)
239 {
240   PetscErrorCode  ierr;
241   PetscScalar     *array;
242 
243   PetscFunctionBegin;
244 #if !defined(PETSC_USE_COMPLEX)
245   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers");
246 #endif
247   ierr = VecGetArray(v,&array);CHKERRQ(ierr);
248   fftw_free((fftw_complex*)array);CHKERRQ(ierr);
249   ierr = VecRestoreArray(v,&array);CHKERRQ(ierr);
250   ierr = VecDestroy_MPI(v);CHKERRQ(ierr);
251   PetscFunctionReturn(0);
252 }
253 
254 #undef __FUNCT__
255 #define __FUNCT__ "MatGetVecs_FFTW"
256 /*
257    MatGetVecs_FFTW - Get vector(s) compatible with the matrix, i.e. with the
258      parallel layout determined by FFTW
259 
260    Collective on Mat
261 
262    Input Parameter:
263 .  mat - the matrix
264 
265    Output Parameter:
266 +   fin - (optional) input vector of forward FFTW
267 -   fout - (optional) output vector of forward FFTW
268 
269   Level: advanced
270 
271 .seealso: MatCreateFFTW()
272 */
273 PetscErrorCode  MatGetVecs_FFTW(Mat A,Vec *fin,Vec *fout)
274 {
275   PetscErrorCode ierr;
276   PetscMPIInt    size,rank;
277   MPI_Comm       comm=((PetscObject)A)->comm;
278   Mat_FFT        *fft = (Mat_FFT*)A->data;
279   PetscInt       N=fft->N;
280 
281   PetscFunctionBegin;
282 #if !defined(PETSC_USE_COMPLEX)
283   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers");
284 #endif
285   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
286   PetscValidType(A,1);
287 
288   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
289   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
290   if (size == 1){ /* sequential case */
291     if (fin) {ierr = VecCreateSeq(PETSC_COMM_SELF,N,fin);CHKERRQ(ierr);}
292     if (fout){ierr = VecCreateSeq(PETSC_COMM_SELF,N,fout);CHKERRQ(ierr);}
293   } else {        /* mpi case */
294     ptrdiff_t      alloc_local,local_n0,local_0_start;
295     ptrdiff_t      local_n1,local_1_end;
296     PetscInt       ndim=fft->ndim,*dim=fft->dim,n=fft->n;
297     fftw_complex   *data_fin,*data_fout;
298 
299     switch (ndim){
300     case 1:
301       /* Get local size */
302 
303       alloc_local = fftw_mpi_local_size_1d(dim[0],comm,FFTW_FORWARD,FFTW_ESTIMATE,&local_n0,&local_0_start,&local_n1,&local_1_end);
304       if (fin) {
305         data_fin  = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
306         ierr = VecCreateMPIWithArray(comm,local_n0,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr);
307         (*fin)->ops->destroy   = VecDestroy_MPIFFTW;
308       }
309       if (fout) {
310         data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
311         ierr = VecCreateMPIWithArray(comm,local_n1,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr);
312         (*fout)->ops->destroy   = VecDestroy_MPIFFTW;
313       }
314       break;
315     case 2:
316       /* Get local size */
317       alloc_local = fftw_mpi_local_size_2d(dim[0],dim[1],comm,&local_n0,&local_0_start);
318       if (fin) {
319         data_fin  = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
320         ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr);
321         (*fin)->ops->destroy   = VecDestroy_MPIFFTW;
322       }
323       if (fout) {
324         data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
325         ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr);
326         (*fout)->ops->destroy   = VecDestroy_MPIFFTW;
327       }
328       break;
329     case 3:
330       /* Get local size */
331       alloc_local = fftw_mpi_local_size_3d(dim[0],dim[1],dim[2],comm,&local_n0,&local_0_start);
332       if (fin) {
333         data_fin  = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
334         ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr);
335         (*fin)->ops->destroy   = VecDestroy_MPIFFTW;
336       }
337       if (fout) {
338         data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
339         ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr);
340         (*fout)->ops->destroy   = VecDestroy_MPIFFTW;
341       }
342       break;
343     default:
344       /* Get local size */
345       alloc_local = fftw_mpi_local_size(ndim,(const ptrdiff_t*)dim,comm,&local_n0,&local_0_start);
346       if (fin) {
347         data_fin  = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
348         ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr);
349         (*fin)->ops->destroy   = VecDestroy_MPIFFTW;
350       }
351       if (fout) {
352         data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
353         ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr);
354         (*fout)->ops->destroy   = VecDestroy_MPIFFTW;
355       }
356       break;
357     }
358   }
359   if (fin){
360     ierr = PetscLayoutReference(A->cmap,&(*fin)->map);CHKERRQ(ierr);
361   }
362   if (fout){
363     ierr = PetscLayoutReference(A->rmap,&(*fout)->map);CHKERRQ(ierr);
364   }
365   PetscFunctionReturn(0);
366 }
367 
368 EXTERN_C_BEGIN
369 #undef __FUNCT__
370 #define __FUNCT__ "MatCreate_FFTW"
371 /*
372       MatCreate_FFTW - Creates a matrix object that provides FFT
373   via the external package FFTW
374 
375   Options Database Keys:
376 + -mat_fftw_plannerflags - set FFTW planner flags
377 
378    Level: intermediate
379 
380 */
381 PetscErrorCode MatCreate_FFTW(Mat A)
382 {
383   PetscErrorCode ierr;
384   MPI_Comm       comm=((PetscObject)A)->comm;
385   Mat_FFT        *fft=(Mat_FFT*)A->data;
386   Mat_FFTW       *fftw;
387   PetscInt       n=fft->n,N=fft->N,ndim=fft->ndim,*dim = fft->dim;
388   const char     *p_flags[]={"FFTW_ESTIMATE","FFTW_MEASURE","FFTW_PATIENT","FFTW_EXHAUSTIVE"};
389   PetscBool      flg;
390   PetscInt       p_flag,partial_dim=1,ctr;
391   PetscMPIInt    size;
392 
393   PetscFunctionBegin;
394 #if !defined(PETSC_USE_COMPLEX)
395   SETERRQ(comm,PETSC_ERR_SUP,"not support for real numbers");
396 #endif
397 
398   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
399   for(ctr=1;ctr<ndim;ctr++){
400           partial_dim*=dim[ctr];
401       }
402   if (size == 1) {
403     ierr = MatSetSizes(A,N,N,N,N);CHKERRQ(ierr);
404     n = N;
405   } else {
406     ptrdiff_t alloc_local,local_n0,local_0_start,local_n1,local_1_end;
407     switch (ndim){
408     case 1:
409       alloc_local = fftw_mpi_local_size_1d(dim[0],comm,FFTW_FORWARD,FFTW_ESTIMATE,&local_n0,&local_0_start,&local_n1,&local_1_end);
410       n = (PetscInt)local_n0;
411       ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr);
412 
413       break;
414     case 2:
415       alloc_local = fftw_mpi_local_size_2d(dim[0],dim[1],comm,&local_n0,&local_0_start);
416       /*
417        PetscMPIInt    rank;
418        PetscSynchronizedPrintf(comm,"[%d] MatCreateSeqFFTW: local_n0, local_0_start %d %d, N %d,dim %d, %d\n",rank,(PetscInt)local_n0*dim[1],(PetscInt)local_0_start,m,dim[0],dim[1]);
419        PetscSynchronizedFlush(comm);
420        */
421       n = (PetscInt)local_n0*dim[1];
422       ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr);
423       break;
424     case 3:
425       alloc_local = fftw_mpi_local_size_3d(dim[0],dim[1],dim[2],comm,&local_n0,&local_0_start);
426       n = (PetscInt)local_n0*dim[1]*dim[2];
427       ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr);
428       break;
429     default:
430       alloc_local = fftw_mpi_local_size(ndim,(const ptrdiff_t*)dim,comm,&local_n0,&local_0_start);
431       n = (PetscInt)local_n0*partial_dim;
432       ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr);
433       break;
434     }
435   }
436   ierr = PetscObjectChangeTypeName((PetscObject)A,MATFFTW);CHKERRQ(ierr);
437 
438   ierr = PetscNewLog(A,Mat_FFTW,&fftw);CHKERRQ(ierr);
439   fft->data = (void*)fftw;
440 
441   fft->n           = n;
442   fftw->p_forward  = 0;
443   fftw->p_backward = 0;
444   fftw->p_flag     = FFTW_ESTIMATE;
445 
446   if (size == 1){
447     A->ops->mult          = MatMult_SeqFFTW;
448     A->ops->multtranspose = MatMultTranspose_SeqFFTW;
449   } else {
450     A->ops->mult          = MatMult_MPIFFTW;
451     A->ops->multtranspose = MatMultTranspose_MPIFFTW;
452   }
453   fft->matdestroy          = MatDestroy_FFTW;
454   A->ops->getvecs       = MatGetVecs_FFTW;
455   A->assembled          = PETSC_TRUE;
456 
457   /* get runtime options */
458   ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"FFTW Options","Mat");CHKERRQ(ierr);
459     ierr = PetscOptionsEList("-mat_fftw_plannerflags","Planner Flags","None",p_flags,4,p_flags[0],&p_flag,&flg);CHKERRQ(ierr);
460     if (flg) {fftw->p_flag = (unsigned)p_flag;}
461   PetscOptionsEnd();
462   PetscFunctionReturn(0);
463 }
464 EXTERN_C_END
465