Actual source code: matusfft.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  2: /*
  3:     Provides an implementation of the Unevenly Sampled FFT algorithm as a Mat.
  4:     Testing examples can be found in ~/src/mat/examples/tests FIX: should these be moved to dm/da/examples/tests?
  5: */

  7:  #include <petsc/private/matimpl.h>
  8: #include <petscdmda.h>                  /*I "petscdmda.h"  I*/ /* Unlike equispaced FFT, USFFT requires geometric information encoded by a DMDA */
  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 on MPI_Comm

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