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 ptrdiff_t ndim_fftw,*dim_fftw; 14 fftw_plan p_forward,p_backward; 15 unsigned p_flag; /* planner flags, FFTW_ESTIMATE,FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE */ 16 PetscScalar *finarray,*foutarray,*binarray,*boutarray; /* keep track of arrays becaue fftw plan should be 17 executed for the arrays with which the plan was created */ 18 } Mat_FFTW; 19 20 extern PetscErrorCode MatMult_SeqFFTW(Mat,Vec,Vec); 21 extern PetscErrorCode MatMultTranspose_SeqFFTW(Mat,Vec,Vec); 22 extern PetscErrorCode MatMult_MPIFFTW(Mat,Vec,Vec); 23 extern PetscErrorCode MatMultTranspose_MPIFFTW(Mat,Vec,Vec); 24 extern PetscErrorCode MatDestroy_FFTW(Mat); 25 extern PetscErrorCode VecDestroy_MPIFFTW(Vec); 26 extern PetscErrorCode MatGetVecs_FFTW(Mat,Vec*,Vec*); 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "MatMult_SeqFFTW" 30 PetscErrorCode MatMult_SeqFFTW(Mat A,Vec x,Vec y) 31 { 32 PetscErrorCode ierr; 33 Mat_FFT *fft = (Mat_FFT*)A->data; 34 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 35 PetscScalar *x_array,*y_array; 36 PetscInt ndim=fft->ndim,*dim=fft->dim; 37 38 PetscFunctionBegin; 39 #if !defined(PETSC_USE_COMPLEX) 40 41 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 42 #endif 43 ierr = VecGetArray(x,&x_array);CHKERRQ(ierr); 44 ierr = VecGetArray(y,&y_array);CHKERRQ(ierr); 45 if (!fftw->p_forward){ /* create a plan, then excute it */ 46 switch (ndim){ 47 case 1: 48 fftw->p_forward = fftw_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_FORWARD,fftw->p_flag); 49 break; 50 case 2: 51 fftw->p_forward = fftw_plan_dft_2d(dim[0],dim[1],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_FORWARD,fftw->p_flag); 52 break; 53 case 3: 54 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); 55 break; 56 default: 57 fftw->p_forward = fftw_plan_dft(ndim,dim,(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_FORWARD,fftw->p_flag); 58 break; 59 } 60 fftw->finarray = x_array; 61 fftw->foutarray = y_array; 62 /* Warning: if (fftw->p_flag!==FFTW_ESTIMATE) The data in the in/out arrays is overwritten! 63 planning should be done before x is initialized! See FFTW manual sec2.1 or sec4 */ 64 fftw_execute(fftw->p_forward); 65 } else { /* use existing plan */ 66 if (fftw->finarray != x_array || fftw->foutarray != y_array){ /* use existing plan on new arrays */ 67 fftw_execute_dft(fftw->p_forward,(fftw_complex*)x_array,(fftw_complex*)y_array); 68 } else { 69 fftw_execute(fftw->p_forward); 70 } 71 } 72 ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr); 73 ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatMultTranspose_SeqFFTW" 79 PetscErrorCode MatMultTranspose_SeqFFTW(Mat A,Vec x,Vec y) 80 { 81 PetscErrorCode ierr; 82 Mat_FFT *fft = (Mat_FFT*)A->data; 83 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 84 PetscScalar *x_array,*y_array; 85 PetscInt ndim=fft->ndim,*dim=fft->dim; 86 87 PetscFunctionBegin; 88 #if !defined(PETSC_USE_COMPLEX) 89 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 90 #endif 91 ierr = VecGetArray(x,&x_array);CHKERRQ(ierr); 92 ierr = VecGetArray(y,&y_array);CHKERRQ(ierr); 93 if (!fftw->p_backward){ /* create a plan, then excute it */ 94 switch (ndim){ 95 case 1: 96 fftw->p_backward = fftw_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_BACKWARD,fftw->p_flag); 97 break; 98 case 2: 99 fftw->p_backward = fftw_plan_dft_2d(dim[0],dim[1],(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_BACKWARD,fftw->p_flag); 100 break; 101 case 3: 102 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); 103 break; 104 default: 105 fftw->p_backward = fftw_plan_dft(ndim,dim,(fftw_complex*)x_array,(fftw_complex*)y_array,FFTW_BACKWARD,fftw->p_flag); 106 break; 107 } 108 fftw->binarray = x_array; 109 fftw->boutarray = y_array; 110 fftw_execute(fftw->p_backward);CHKERRQ(ierr); 111 } else { /* use existing plan */ 112 if (fftw->binarray != x_array || fftw->boutarray != y_array){ /* use existing plan on new arrays */ 113 fftw_execute_dft(fftw->p_backward,(fftw_complex*)x_array,(fftw_complex*)y_array); 114 } else { 115 fftw_execute(fftw->p_backward);CHKERRQ(ierr); 116 } 117 } 118 ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr); 119 ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 #undef __FUNCT__ 124 #define __FUNCT__ "MatMult_MPIFFTW" 125 PetscErrorCode MatMult_MPIFFTW(Mat A,Vec x,Vec y) 126 { 127 PetscErrorCode ierr; 128 Mat_FFT *fft = (Mat_FFT*)A->data; 129 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 130 PetscScalar *x_array,*y_array; 131 PetscInt ndim=fft->ndim,*dim=fft->dim; 132 MPI_Comm comm=((PetscObject)A)->comm; 133 // PetscInt ctr; 134 // ptrdiff_t ndim1=(ptrdiff_t) ndim,*pdim; 135 // ndim1=(ptrdiff_t) ndim; 136 // pdim = (ptrdiff_t *)calloc(ndim,sizeof(ptrdiff_t)); 137 138 // for(ctr=0;ctr<ndim;ctr++) 139 // { 140 // pdim[ctr] = dim[ctr]; 141 // } 142 143 PetscFunctionBegin; 144 //#if !defined(PETSC_USE_COMPLEX) 145 // SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 146 //#endif 147 // pdim = (ptrdiff_t *)calloc(ndim,sizeof(ptrdiff_t)); 148 // for (ctr=0; ctr<ndim; ctr++) pdim[ctr] = dim[ctr]; 149 150 ierr = VecGetArray(x,&x_array);CHKERRQ(ierr); 151 ierr = VecGetArray(y,&y_array);CHKERRQ(ierr); 152 if (!fftw->p_forward){ /* create a plan, then excute it */ 153 switch (ndim){ 154 case 1: 155 #if defined(PETSC_USE_COMPLEX) 156 fftw->p_forward = fftw_mpi_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_FORWARD,fftw->p_flag); 157 #endif 158 break; 159 case 2: 160 #if defined(PETSC_USE_COMPLEX) 161 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); 162 #else 163 printf("The code comes here \n"); 164 fftw->p_forward = fftw_mpi_plan_dft_r2c_2d(dim[0],dim[1],(double *)x_array,(fftw_complex*)y_array,comm,FFTW_ESTIMATE); 165 #endif 166 break; 167 case 3: 168 #if defined(PETSC_USE_COMPLEX) 169 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); 170 #else 171 fftw->p_forward = fftw_mpi_plan_dft_r2c_3d(dim[0],dim[1],dim[2],(double *)x_array,(fftw_complex*)y_array,comm,FFTW_ESTIMATE); 172 #endif 173 break; 174 default: 175 #if defined(PETSC_USE_COMPLEX) 176 fftw->p_forward = fftw_mpi_plan_dft(fftw->ndim_fftw,fftw->dim_fftw,(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_FORWARD,fftw->p_flag); 177 #else 178 fftw->p_forward = fftw_mpi_plan_dft_r2c(fftw->ndim_fftw,fftw->dim_fftw,(double *)x_array,(fftw_complex*)y_array,comm,FFTW_ESTIMATE); 179 #endif 180 // fftw->p_forward = fftw_mpi_plan_dft(ndim,dim,(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_FORWARD,fftw->p_flag); 181 break; 182 } 183 fftw->finarray = x_array; 184 fftw->foutarray = y_array; 185 /* Warning: if (fftw->p_flag!==FFTW_ESTIMATE) The data in the in/out arrays is overwritten! 186 planning should be done before x is initialized! See FFTW manual sec2.1 or sec4 */ 187 fftw_execute(fftw->p_forward); 188 } else { /* use existing plan */ 189 if (fftw->finarray != x_array || fftw->foutarray != y_array){ /* use existing plan on new arrays */ 190 fftw_execute_dft(fftw->p_forward,(fftw_complex*)x_array,(fftw_complex*)y_array); 191 } else { 192 fftw_execute(fftw->p_forward); 193 } 194 } 195 ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr); 196 ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "MatMultTranspose_MPIFFTW" 202 PetscErrorCode MatMultTranspose_MPIFFTW(Mat A,Vec x,Vec y) 203 { 204 PetscErrorCode ierr; 205 Mat_FFT *fft = (Mat_FFT*)A->data; 206 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 207 PetscScalar *x_array,*y_array; 208 PetscInt ndim=fft->ndim,*dim=fft->dim; 209 MPI_Comm comm=((PetscObject)A)->comm; 210 // PetscInt ctr; 211 // ptrdiff_t ndim1=(ptrdiff_t)ndim,*pdim; 212 // ndim1=(ptrdiff_t) ndim; 213 // pdim = (ptrdiff_t *)calloc(ndim,sizeof(ptrdiff_t)); 214 215 // for(ctr=0;ctr<ndim;ctr++) 216 // { 217 // pdim[ctr] = dim[ctr]; 218 // } 219 220 PetscFunctionBegin; 221 //#if !defined(PETSC_USE_COMPLEX) 222 // SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 223 //#endif 224 // ierr = PetscMalloc(ndim*sizeof(ptrdiff_t), (ptrdiff_t *)&pdim);CHKERRQ(ierr); 225 // should pdim be a member of Mat_FFTW? 226 // for (ctr=0; ctr<ndim; ctr++) pdim[ctr] = dim[ctr]; 227 228 ierr = VecGetArray(x,&x_array);CHKERRQ(ierr); 229 ierr = VecGetArray(y,&y_array);CHKERRQ(ierr); 230 if (!fftw->p_backward){ /* create a plan, then excute it */ 231 switch (ndim){ 232 case 1: 233 #if defined(PETSC_USE_COMPLEX) 234 fftw->p_backward = fftw_mpi_plan_dft_1d(dim[0],(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_BACKWARD,fftw->p_flag); 235 #endif 236 break; 237 case 2: 238 #if defined(PETSC_USE_COMPLEX) 239 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); 240 #else 241 fftw->p_backward = fftw_mpi_plan_dft_c2r_2d(dim[0],dim[1],(fftw_complex*)x_array,(double *)y_array,comm,FFTW_ESTIMATE); 242 #endif 243 break; 244 case 3: 245 #if defined(PETSC_USE_COMPLEX) 246 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); 247 #else 248 fftw->p_backward = fftw_mpi_plan_dft_c2r_3d(dim[0],dim[1],dim[2],(fftw_complex*)x_array,(double *)y_array,comm,FFTW_ESTIMATE); 249 #endif 250 break; 251 default: 252 #if defined(PETSC_USE_COMPLEX) 253 fftw->p_backward = fftw_mpi_plan_dft(fftw->ndim_fftw,fftw->dim_fftw,(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_BACKWARD,fftw->p_flag); 254 #else 255 fftw->p_backward = fftw_mpi_plan_dft_c2r(fftw->ndim_fftw,fftw->dim_fftw,(fftw_complex*)x_array,(double *)y_array,comm,FFTW_ESTIMATE); 256 #endif 257 // fftw->p_backward = fftw_mpi_plan_dft(ndim,dim,(fftw_complex*)x_array,(fftw_complex*)y_array,comm,FFTW_BACKWARD,fftw->p_flag); 258 break; 259 } 260 fftw->binarray = x_array; 261 fftw->boutarray = y_array; 262 fftw_execute(fftw->p_backward);CHKERRQ(ierr); 263 } else { /* use existing plan */ 264 if (fftw->binarray != x_array || fftw->boutarray != y_array){ /* use existing plan on new arrays */ 265 fftw_execute_dft(fftw->p_backward,(fftw_complex*)x_array,(fftw_complex*)y_array); 266 } else { 267 fftw_execute(fftw->p_backward);CHKERRQ(ierr); 268 } 269 } 270 ierr = VecRestoreArray(y,&y_array);CHKERRQ(ierr); 271 ierr = VecRestoreArray(x,&x_array);CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274 275 #undef __FUNCT__ 276 #define __FUNCT__ "MatDestroy_FFTW" 277 PetscErrorCode MatDestroy_FFTW(Mat A) 278 { 279 Mat_FFT *fft = (Mat_FFT*)A->data; 280 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 281 PetscErrorCode ierr; 282 283 PetscFunctionBegin; 284 #if !defined(PETSC_USE_COMPLEX) 285 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 286 #endif 287 fftw_destroy_plan(fftw->p_forward); 288 fftw_destroy_plan(fftw->p_backward); 289 ierr = PetscFree(fftw->dim_fftw);CHKERRQ(ierr); 290 ierr = PetscFree(fft->data);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 #include <../src/vec/vec/impls/mpi/pvecimpl.h> /*I "petscvec.h" I*/ 295 #undef __FUNCT__ 296 #define __FUNCT__ "VecDestroy_MPIFFTW" 297 PetscErrorCode VecDestroy_MPIFFTW(Vec v) 298 { 299 PetscErrorCode ierr; 300 PetscScalar *array; 301 302 PetscFunctionBegin; 303 #if !defined(PETSC_USE_COMPLEX) 304 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 305 #endif 306 ierr = VecGetArray(v,&array);CHKERRQ(ierr); 307 fftw_free((fftw_complex*)array);CHKERRQ(ierr); 308 ierr = VecRestoreArray(v,&array);CHKERRQ(ierr); 309 ierr = VecDestroy_MPI(v);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312 313 #undef __FUNCT__ 314 #define __FUNCT__ "MatGetVecs1DC_FFTW" 315 /* 316 MatGetVecs_FFTW1D - Get Vectors(s) compatible with matrix, i.e. with the 317 parallel layout determined by FFTW-1D 318 319 */ 320 PetscErrorCode MatGetVecs_FFTW1D(Mat A,Vec *fin,Vec *fout,Vec *bin,Vec *bout) 321 { 322 PetscErrorCode ierr; 323 PetscMPIInt size,rank; 324 MPI_Comm comm=((PetscObject)A)->comm; 325 Mat_FFT *fft = (Mat_FFT*)A->data; 326 // Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 327 PetscInt N=fft->N; 328 PetscInt ndim=fft->ndim,*dim=fft->dim; 329 ptrdiff_t f_alloc_local,f_local_n0,f_local_0_start; 330 ptrdiff_t f_local_n1,f_local_1_end; 331 ptrdiff_t b_alloc_local,b_local_n0,b_local_0_start; 332 ptrdiff_t b_local_n1,b_local_1_end; 333 fftw_complex *data_fin,*data_fout,*data_bin,*data_bout; 334 335 PetscFunctionBegin; 336 #if !defined(PETSC_USE_COMPLEX) 337 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 338 #endif 339 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 340 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 341 if (size == 1){ 342 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Works only for parallel 1D"); 343 } 344 else { 345 if (ndim>1){ 346 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Works only for parallel 1D");} 347 else { 348 f_alloc_local = fftw_mpi_local_size_1d(dim[0],comm,FFTW_FORWARD,FFTW_ESTIMATE,&f_local_n0,&f_local_0_start,&f_local_n1,&f_local_1_end); 349 b_alloc_local = fftw_mpi_local_size_1d(dim[0],comm,FFTW_BACKWARD,FFTW_ESTIMATE,&b_local_n0,&b_local_0_start,&b_local_n1,&b_local_1_end); 350 if (fin) { 351 data_fin = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*f_alloc_local); 352 ierr = VecCreateMPIWithArray(comm,f_local_n0,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr); 353 (*fin)->ops->destroy = VecDestroy_MPIFFTW; 354 } 355 if (fout) { 356 data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*f_alloc_local); 357 ierr = VecCreateMPIWithArray(comm,f_local_n1,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr); 358 (*fout)->ops->destroy = VecDestroy_MPIFFTW; 359 } 360 if (bin) { 361 data_bin = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*b_alloc_local); 362 ierr = VecCreateMPIWithArray(comm,b_local_n0,N,(const PetscScalar*)data_bin,bin);CHKERRQ(ierr); 363 (*bin)->ops->destroy = VecDestroy_MPIFFTW; 364 } 365 if (bout) { 366 data_bout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*b_alloc_local); 367 ierr = VecCreateMPIWithArray(comm,b_local_n1,N,(const PetscScalar*)data_bout,bout);CHKERRQ(ierr); 368 (*bout)->ops->destroy = VecDestroy_MPIFFTW; 369 } 370 } 371 if (fin){ 372 ierr = PetscLayoutReference(A->cmap,&(*fin)->map);CHKERRQ(ierr); 373 } 374 if (fout){ 375 ierr = PetscLayoutReference(A->rmap,&(*fout)->map);CHKERRQ(ierr); 376 } 377 if (bin){ 378 ierr = PetscLayoutReference(A->rmap,&(*bin)->map);CHKERRQ(ierr); 379 } 380 if (bout){ 381 ierr = PetscLayoutReference(A->rmap,&(*bout)->map);CHKERRQ(ierr); 382 } 383 PetscFunctionReturn(0); 384 } 385 386 387 } 388 389 #undef __FUNCT__ 390 #define __FUNCT__ "MatGetVecs_FFTW" 391 /* 392 MatGetVecs_FFTW - Get vector(s) compatible with the matrix, i.e. with the 393 parallel layout determined by FFTW 394 395 Collective on Mat 396 397 Input Parameter: 398 . mat - the matrix 399 400 Output Parameter: 401 + fin - (optional) input vector of forward FFTW 402 - fout - (optional) output vector of forward FFTW 403 404 Level: advanced 405 406 .seealso: MatCreateFFTW() 407 */ 408 PetscErrorCode MatGetVecs_FFTW(Mat A,Vec *fin,Vec *fout) 409 { 410 PetscErrorCode ierr; 411 PetscMPIInt size,rank; 412 MPI_Comm comm=((PetscObject)A)->comm; 413 Mat_FFT *fft = (Mat_FFT*)A->data; 414 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 415 PetscInt N=fft->N, N1, n1,vsize; 416 417 PetscFunctionBegin; 418 //#if !defined(PETSC_USE_COMPLEX) 419 // SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not support for real numbers"); 420 //#endif 421 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 422 PetscValidType(A,1); 423 424 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 425 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 426 if (size == 1){ /* sequential case */ 427 if (fin) {ierr = VecCreateSeq(PETSC_COMM_SELF,N,fin);CHKERRQ(ierr);} 428 if (fout){ierr = VecCreateSeq(PETSC_COMM_SELF,N,fout);CHKERRQ(ierr);} 429 printf("The code successfully comes at the end of the routine with one processor\n"); 430 } else { /* mpi case */ 431 ptrdiff_t alloc_local,local_n0,local_0_start; 432 ptrdiff_t local_n1,local_1_end; 433 PetscInt ndim=fft->ndim,*dim=fft->dim,n=fft->n; 434 fftw_complex *data_fin,*data_fout; 435 double *data_finr ; 436 ptrdiff_t local_1_start; 437 // PetscInt ctr; 438 // ptrdiff_t ndim1,*pdim; 439 // ndim1=(ptrdiff_t) ndim; 440 // pdim = (ptrdiff_t *)calloc(ndim,sizeof(ptrdiff_t)); 441 442 // for(ctr=0;ctr<ndim;ctr++) 443 // {k 444 // pdim[ctr] = dim[ctr]; 445 // } 446 447 448 449 switch (ndim){ 450 case 1: 451 /* Get local size */ 452 /* We need to write an error message here saying that one cannot call this routine when doing paralllel 1D complex FFTW */ 453 // SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Works only for parallel Multi-dimensional FFTW, Dimension>1. Check Documentation for MatGetVecs_FFTW1D routine"); 454 alloc_local = fftw_mpi_local_size_1d(dim[0],comm,FFTW_FORWARD,FFTW_ESTIMATE,&local_n0,&local_0_start,&local_n1,&local_1_end); 455 if (fin) { 456 data_fin = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 457 ierr = VecCreateMPIWithArray(comm,local_n0,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr); 458 (*fin)->ops->destroy = VecDestroy_MPIFFTW; 459 } 460 if (fout) { 461 data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 462 ierr = VecCreateMPIWithArray(comm,local_n1,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr); 463 (*fout)->ops->destroy = VecDestroy_MPIFFTW; 464 } 465 break; 466 case 2: 467 #if !defined(PETSC_USE_COMPLEX) 468 alloc_local = fftw_mpi_local_size_2d_transposed(dim[0],dim[1]/2+1,comm,&local_n0,&local_0_start,&local_n1,&local_1_start); 469 N1 = 2*dim[0]*(dim[1]/2+1); n1 = 2*local_n0*(dim[1]/2+1); 470 if (fin) { 471 data_finr=(double *)fftw_malloc(sizeof(double)*alloc_local*2); 472 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,(PetscInt)n1,N1,(PetscScalar*)data_finr,fin);CHKERRQ(ierr); 473 ierr = VecGetSize(*fin,&vsize);CHKERRQ(ierr); 474 //printf("The code comes here with vector size %d\n",vsize); 475 (*fin)->ops->destroy = VecDestroy_MPIFFTW; 476 } 477 if (fout) { 478 data_fout=(fftw_complex *)fftw_malloc(sizeof(fftw_complex)*alloc_local); 479 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,n1,N1,(PetscScalar*)data_fout,fout);CHKERRQ(ierr); 480 (*fout)->ops->destroy = VecDestroy_MPIFFTW; 481 } 482 printf("Vector size from fftw.c is given by %d, %d\n",n1,N1); 483 484 #else 485 /* Get local size */ 486 printf("Hope this does not come here"); 487 alloc_local = fftw_mpi_local_size_2d(dim[0],dim[1],comm,&local_n0,&local_0_start); 488 if (fin) { 489 data_fin = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 490 ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr); 491 (*fin)->ops->destroy = VecDestroy_MPIFFTW; 492 } 493 if (fout) { 494 data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 495 ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr); 496 (*fout)->ops->destroy = VecDestroy_MPIFFTW; 497 } 498 printf("Hope this does not come here"); 499 #endif 500 break; 501 case 3: 502 /* Get local size */ 503 #if !defined(PETSC_USE_COMPLEX) 504 alloc_local = fftw_mpi_local_size_3d_transposed(dim[0],dim[1],dim[2]/2+1,comm,&local_n0,&local_0_start,&local_n1,&local_1_start); 505 N1 = 2*dim[0]*dim[1]*(dim[2]/2+1); n1 = 2*local_n0*dim[1]*(dim[2]/2+1); 506 if (fin) { 507 data_finr=(double *)fftw_malloc(sizeof(double)*alloc_local*2); 508 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,(PetscInt)n1,N1,(PetscScalar*)data_finr,fin);CHKERRQ(ierr); 509 ierr = VecGetSize(*fin,&vsize);CHKERRQ(ierr); 510 //printf("The code comes here with vector size %d\n",vsize); 511 (*fin)->ops->destroy = VecDestroy_MPIFFTW; 512 } 513 if (fout) { 514 data_fout=(fftw_complex *)fftw_malloc(sizeof(fftw_complex)*alloc_local); 515 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,n1,N1,(PetscScalar*)data_fout,fout);CHKERRQ(ierr); 516 (*fout)->ops->destroy = VecDestroy_MPIFFTW; 517 } 518 printf("Vector size from fftw.c is given by %d, %d\n",n1,N1); 519 520 521 #else 522 alloc_local = fftw_mpi_local_size_3d(dim[0],dim[1],dim[2],comm,&local_n0,&local_0_start); 523 // printf("The quantity n is %d",n); 524 if (fin) { 525 data_fin = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 526 ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr); 527 (*fin)->ops->destroy = VecDestroy_MPIFFTW; 528 } 529 if (fout) { 530 data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 531 ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr); 532 (*fout)->ops->destroy = VecDestroy_MPIFFTW; 533 } 534 #endif 535 break; 536 default: 537 /* Get local size */ 538 #if !defined(PETSC_USE_COMPLEX) 539 SETERRQ(comm,PETSC_ERR_SUP,"Not done yet"); 540 #else 541 alloc_local = fftw_mpi_local_size(fftw->ndim_fftw,fftw->dim_fftw,comm,&local_n0,&local_0_start); 542 // printf("The value of alloc local is %d from process %d\n",alloc_local,rank); 543 // printf("The value of alloc local is %d",alloc_local); 544 // pdim=(ptrdiff_t *)calloc(ndim,sizeof(ptrdiff_t)); 545 // for(i=0;i<ndim;i++) 546 // { 547 // pdim[i]=dim[i];printf("%d",pdim[i]); 548 // } 549 // alloc_local = fftw_mpi_local_size(ndim,pdim,comm,&local_n0,&local_0_start); 550 // printf("The quantity n is %d",n); 551 if (fin) { 552 data_fin = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 553 ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fin,fin);CHKERRQ(ierr); 554 (*fin)->ops->destroy = VecDestroy_MPIFFTW; 555 } 556 if (fout) { 557 data_fout = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 558 ierr = VecCreateMPIWithArray(comm,n,N,(const PetscScalar*)data_fout,fout);CHKERRQ(ierr); 559 (*fout)->ops->destroy = VecDestroy_MPIFFTW; 560 } 561 #endif 562 break; 563 } 564 } 565 // if (fin){ 566 // ierr = PetscLayoutReference(A->cmap,&(*fin)->map);CHKERRQ(ierr); 567 // } 568 // if (fout){ 569 // ierr = PetscLayoutReference(A->rmap,&(*fout)->map);CHKERRQ(ierr); 570 // } 571 PetscFunctionReturn(0); 572 } 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "InputTransformFFT" 576 PetscErrorCode InputTransformFFT(Mat A,Vec x,Vec y) 577 { 578 PetscErrorCode ierr; 579 PetscFunctionBegin; 580 ierr = PetscTryMethod(A,"InputTransformFFT_C",(Mat,Vec,Vec),(A,x,y));CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 584 /* 585 InputTransformFFT_FFTW - Copies the user data to the vector that goes into FFTW block 586 Input A, x, y 587 A - FFTW matrix 588 x - user data 589 Options Database Keys: 590 + -mat_fftw_plannerflags - set FFTW planner flags 591 592 Level: intermediate 593 594 */ 595 596 EXTERN_C_BEGIN 597 #undef __FUNCT__ 598 #define __FUNCT__ "InputTransformFFT_FTTW" 599 PetscErrorCode InputTransformFFT_FFTW(Mat A,Vec x,Vec y) 600 { 601 PetscErrorCode ierr; 602 MPI_Comm comm=((PetscObject)A)->comm; 603 Mat_FFT *fft = (Mat_FFT*)A->data; 604 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 605 PetscInt N=fft->N, N1, n1 ,NM; 606 PetscInt ndim=fft->ndim,*dim=fft->dim,n=fft->n; 607 PetscInt low, *indx1, *indx2, tempindx, tempindx1, *indx3, *indx4; 608 PetscInt i,j,k,rank,size; 609 ptrdiff_t alloc_local,local_n0,local_0_start; 610 ptrdiff_t local_n1,local_1_start; 611 VecScatter vecscat; 612 IS list1,list2; 613 614 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 615 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 616 ierr = VecGetOwnershipRange(y,&low,PETSC_NULL); 617 printf("Local ownership starts at %d\n",low); 618 619 switch (ndim){ 620 case 1: 621 SETERRQ(comm,PETSC_ERR_SUP,"Not Supported by FFTW"); 622 break; 623 case 2: 624 alloc_local = fftw_mpi_local_size_2d_transposed(dim[0],dim[1]/2+1,comm,&local_n0,&local_0_start,&local_n1,&local_1_start); 625 N1 = 2*dim[0]*(dim[1]/2+1); n1 = 2*local_n0*(dim[1]/2+1); 626 627 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx1);CHKERRQ(ierr); 628 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx2);CHKERRQ(ierr); 629 printf("Val local_0_start = %ld\n",local_0_start); 630 631 if (dim[1]%2==0) 632 NM = dim[1]+2; 633 else 634 NM = dim[1]+1; 635 636 for (i=0;i<local_n0;i++){ 637 for (j=0;j<dim[1];j++){ 638 tempindx = i*dim[1] + j; 639 tempindx1 = i*NM + j; 640 indx1[tempindx]=local_0_start*dim[1]+tempindx; 641 indx2[tempindx]=low+tempindx1; 642 // printf("Val tempindx1 = %d\n",tempindx1); 643 // printf("index1 %d from proc %d is \n",indx1[tempindx],rank); 644 // printf("index2 %d from proc %d is \n",indx2[tempindx],rank); 645 // printf("-------------------------\n",indx2[tempindx],rank); 646 } 647 } 648 649 ierr = ISCreateGeneral(comm,local_n0*dim[1],indx1,PETSC_COPY_VALUES,&list1);CHKERRQ(ierr); 650 ierr = ISCreateGeneral(comm,local_n0*dim[1],indx2,PETSC_COPY_VALUES,&list2);CHKERRQ(ierr); 651 652 ierr = VecScatterCreate(x,list1,y,list2,&vecscat);CHKERRQ(ierr); 653 ierr = VecScatterBegin(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 654 ierr = VecScatterEnd(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 655 ierr = VecScatterDestroy(&vecscat);CHKERRQ(ierr); 656 break; 657 658 case 3: 659 alloc_local = fftw_mpi_local_size_3d_transposed(dim[0],dim[1],dim[2]/2+1,comm,&local_n0,&local_0_start,&local_n1,&local_1_start); 660 N1 = 2*dim[0]*dim[1]*(dim[2]/2+1); n1 = 2*local_n0*dim[1]*(dim[2]/2+1); 661 662 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx1);CHKERRQ(ierr); 663 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx2);CHKERRQ(ierr); 664 printf("Val local_0_start = %ld\n",local_0_start); 665 666 if (dim[2]%2==0) 667 NM = dim[2]+2; 668 else 669 NM = dim[2]+1; 670 671 for (i=0;i<local_n0;i++){ 672 for (j=0;j<dim[1];j++){ 673 for (k=0;k<dim[2];k++){ 674 tempindx = i*dim[1]*dim[2] + j*dim[2] + k; 675 tempindx1 = i*dim[1]*NM + j*NM + k; 676 indx1[tempindx]=local_0_start*dim[1]*dim[2]+tempindx; 677 indx2[tempindx]=low+tempindx1; 678 } 679 // printf("Val tempindx1 = %d\n",tempindx1); 680 // printf("index1 %d from proc %d is \n",indx1[tempindx],rank); 681 // printf("index2 %d from proc %d is \n",indx2[tempindx],rank); 682 // printf("-------------------------\n",indx2[tempindx],rank); 683 } 684 } 685 686 ierr = ISCreateGeneral(comm,local_n0*dim[1]*dim[2],indx1,PETSC_COPY_VALUES,&list1);CHKERRQ(ierr); 687 ierr = ISCreateGeneral(comm,local_n0*dim[1]*dim[2],indx2,PETSC_COPY_VALUES,&list2);CHKERRQ(ierr); 688 689 ierr = VecScatterCreate(x,list1,y,list2,&vecscat);CHKERRQ(ierr); 690 ierr = VecScatterBegin(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 691 ierr = VecScatterEnd(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 692 ierr = VecScatterDestroy(&vecscat);CHKERRQ(ierr); 693 break; 694 695 default: 696 SETERRQ(comm,PETSC_ERR_SUP,"Not Done Yet"); 697 break; 698 } 699 700 return 0; 701 } 702 EXTERN_C_END 703 704 #undef __FUNCT__ 705 #define __FUNCT__ "OutputTransformFFT" 706 PetscErrorCode OutputTransformFFT(Mat A,Vec x,Vec y) 707 { 708 PetscErrorCode ierr; 709 PetscFunctionBegin; 710 ierr = PetscTryMethod(A,"OutputTransformFFT_C",(Mat,Vec,Vec),(A,x,y));CHKERRQ(ierr); 711 PetscFunctionReturn(0); 712 } 713 714 /* 715 OutputTransformFFT_FFTW - Copies the FFTW output to the PETSc vector that user can use 716 Input A, x, y 717 A - FFTW matrix 718 x - FFTW vector 719 y - PETSc vector that user can read 720 Options Database Keys: 721 + -mat_fftw_plannerflags - set FFTW planner flags 722 723 Level: intermediate 724 725 */ 726 727 EXTERN_C_BEGIN 728 #undef __FUNCT__ 729 #define __FUNCT__ "OutputTransformFFT_FTTW" 730 PetscErrorCode OutputTransformFFT_FFTW(Mat A,Vec x,Vec y) 731 { 732 PetscErrorCode ierr; 733 MPI_Comm comm=((PetscObject)A)->comm; 734 Mat_FFT *fft = (Mat_FFT*)A->data; 735 Mat_FFTW *fftw = (Mat_FFTW*)fft->data; 736 PetscInt N=fft->N, N1, n1 ,NM; 737 PetscInt ndim=fft->ndim,*dim=fft->dim,n=fft->n; 738 PetscInt low, *indx1, *indx2, tempindx, tempindx1, *indx3, *indx4; 739 PetscInt i,j,k,rank,size; 740 ptrdiff_t alloc_local,local_n0,local_0_start; 741 ptrdiff_t local_n1,local_1_start; 742 VecScatter vecscat; 743 IS list1,list2; 744 745 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 746 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 747 ierr = VecGetOwnershipRange(x,&low,PETSC_NULL); 748 printf("Local ownership starts at %d\n",low); 749 750 switch (ndim){ 751 case 1: 752 SETERRQ(comm,PETSC_ERR_SUP,"Not Supported by FFTW"); 753 break; 754 case 2: 755 alloc_local = fftw_mpi_local_size_2d_transposed(dim[0],dim[1]/2+1,comm,&local_n0,&local_0_start,&local_n1,&local_1_start); 756 N1 = 2*dim[0]*(dim[1]/2+1); n1 = 2*local_n0*(dim[1]/2+1); 757 758 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx1);CHKERRQ(ierr); 759 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx2);CHKERRQ(ierr); 760 printf("Val local_0_start = %ld\n",local_0_start); 761 762 if (dim[1]%2==0) 763 NM = dim[1]+2; 764 else 765 NM = dim[1]+1; 766 767 768 769 for (i=0;i<local_n0;i++){ 770 for (j=0;j<dim[1];j++){ 771 tempindx = i*dim[1] + j; 772 tempindx1 = i*NM + j; 773 indx1[tempindx]=local_0_start*dim[1]+tempindx; 774 indx2[tempindx]=low+tempindx1; 775 // printf("Val tempindx1 = %d\n",tempindx1); 776 // printf("index1 %d from proc %d is \n",indx1[tempindx],rank); 777 // printf("index2 %d from proc %d is \n",indx2[tempindx],rank); 778 // printf("-------------------------\n",indx2[tempindx],rank); 779 } 780 } 781 782 ierr = ISCreateGeneral(comm,local_n0*dim[1],indx1,PETSC_COPY_VALUES,&list1);CHKERRQ(ierr); 783 ierr = ISCreateGeneral(comm,local_n0*dim[1],indx2,PETSC_COPY_VALUES,&list2);CHKERRQ(ierr); 784 785 ierr = VecScatterCreate(x,list2,y,list1,&vecscat);CHKERRQ(ierr); 786 ierr = VecScatterBegin(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 787 ierr = VecScatterEnd(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 788 ierr = VecScatterDestroy(&vecscat);CHKERRQ(ierr); 789 break; 790 791 case 3: 792 alloc_local = fftw_mpi_local_size_3d_transposed(dim[0],dim[1],dim[2]/2+1,comm,&local_n0,&local_0_start,&local_n1,&local_1_start); 793 N1 = 2*dim[0]*dim[1]*(dim[2]/2+1); n1 = 2*local_n0*dim[1]*(dim[2]/2+1); 794 795 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx1);CHKERRQ(ierr); 796 ierr = PetscMalloc(sizeof(PetscInt)*((PetscInt)local_n0)*N1,&indx2);CHKERRQ(ierr); 797 printf("Val local_0_start = %ld\n",local_0_start); 798 799 if (dim[2]%2==0) 800 NM = dim[2]+2; 801 else 802 NM = dim[2]+1; 803 804 for (i=0;i<local_n0;i++){ 805 for (j=0;j<dim[1];j++){ 806 for (k=0;k<dim[2];k++){ 807 tempindx = i*dim[1]*dim[2] + j*dim[2] + k; 808 tempindx1 = i*dim[1]*NM + j*NM + k; 809 indx1[tempindx]=local_0_start*dim[1]*dim[2]+tempindx; 810 indx2[tempindx]=low+tempindx1; 811 } 812 // printf("Val tempindx1 = %d\n",tempindx1); 813 // printf("index1 %d from proc %d is \n",indx1[tempindx],rank); 814 // printf("index2 %d from proc %d is \n",indx2[tempindx],rank); 815 // printf("-------------------------\n",indx2[tempindx],rank); 816 } 817 } 818 819 ierr = ISCreateGeneral(comm,local_n0*dim[1]*dim[2],indx1,PETSC_COPY_VALUES,&list1);CHKERRQ(ierr); 820 ierr = ISCreateGeneral(comm,local_n0*dim[1]*dim[2],indx2,PETSC_COPY_VALUES,&list2);CHKERRQ(ierr); 821 822 ierr = VecScatterCreate(x,list2,y,list1,&vecscat);CHKERRQ(ierr); 823 ierr = VecScatterBegin(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 824 ierr = VecScatterEnd(vecscat,x,y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 825 ierr = VecScatterDestroy(&vecscat);CHKERRQ(ierr); 826 break; 827 828 default: 829 SETERRQ(comm,PETSC_ERR_SUP,"Not Done Yet"); 830 break; 831 } 832 return 0; 833 } 834 EXTERN_C_END 835 836 EXTERN_C_BEGIN 837 #undef __FUNCT__ 838 #define __FUNCT__ "MatCreate_FFTW" 839 /* 840 MatCreate_FFTW - Creates a matrix object that provides FFT 841 via the external package FFTW 842 Options Database Keys: 843 + -mat_fftw_plannerflags - set FFTW planner flags 844 845 Level: intermediate 846 847 */ 848 849 PetscErrorCode MatCreate_FFTW(Mat A) 850 { 851 PetscErrorCode ierr; 852 MPI_Comm comm=((PetscObject)A)->comm; 853 Mat_FFT *fft=(Mat_FFT*)A->data; 854 Mat_FFTW *fftw; 855 PetscInt n=fft->n,N=fft->N,ndim=fft->ndim,*dim = fft->dim; 856 const char *p_flags[]={"FFTW_ESTIMATE","FFTW_MEASURE","FFTW_PATIENT","FFTW_EXHAUSTIVE"}; 857 PetscBool flg; 858 PetscInt p_flag,partial_dim=1,ctr; 859 PetscMPIInt size,rank; 860 ptrdiff_t *pdim; 861 ptrdiff_t local_n1,local_1_start; 862 863 PetscFunctionBegin; 864 //#if !defined(PETSC_USE_COMPLEX) 865 // SETERRQ(comm,PETSC_ERR_SUP,"not support for real numbers"); 866 //#endif 867 868 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 869 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 870 ierr = MPI_Barrier(PETSC_COMM_WORLD); 871 872 pdim = (ptrdiff_t *)calloc(ndim,sizeof(ptrdiff_t)); 873 pdim[0] = dim[0]; 874 for(ctr=1;ctr<ndim;ctr++) 875 { 876 partial_dim *= dim[ctr]; 877 pdim[ctr] = dim[ctr]; 878 } 879 //#if !defined(PETSC_USE_COMPLEX) 880 // SETERRQ(comm,PETSC_ERR_SUP,"not support for real numbers"); 881 //#endif 882 883 // printf("partial dimension is %d",partial_dim); 884 if (size == 1) { 885 ierr = MatSetSizes(A,N,N,N,N);CHKERRQ(ierr); 886 n = N; 887 } else { 888 ptrdiff_t alloc_local,local_n0,local_0_start,local_n1,local_1_end; 889 switch (ndim){ 890 case 1: 891 #if !defined(PETSC_USE_COMPLEX) 892 SETERRQ(comm,PETSC_ERR_SUP,"FFTW does not support parallel 1D real transform"); 893 #endif 894 alloc_local = fftw_mpi_local_size_1d(dim[0],comm,FFTW_FORWARD,FFTW_ESTIMATE,&local_n0,&local_0_start,&local_n1,&local_1_end); 895 n = (PetscInt)local_n0; 896 ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr); 897 // PetscObjectComposeFunctionDynamic((PetscObject)A,"MatGetVecs1DC_C","MatGetVecs1DC_FFTW",MatGetVecs1DC_FFTW); 898 break; 899 case 2: 900 #if defined(PETSC_USE_COMPLEX) 901 alloc_local = fftw_mpi_local_size_2d(dim[0],dim[1],comm,&local_n0,&local_0_start); 902 /* 903 PetscMPIInt rank; 904 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]); 905 PetscSynchronizedFlush(comm); 906 */ 907 n = (PetscInt)local_n0*dim[1]; 908 ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr); 909 #else 910 alloc_local = fftw_mpi_local_size_2d_transposed(dim[0],dim[1]/2+1,PETSC_COMM_WORLD,&local_n0,&local_0_start,&local_n1,&local_1_start); 911 n = 2*(PetscInt)local_n0*(dim[1]/2+1); 912 ierr = MatSetSizes(A,n,n,2*dim[0]*(dim[1]/2+1),2*dim[0]*(dim[1]/2+1)); 913 #endif 914 break; 915 case 3: 916 // printf("The value of alloc local is %d",alloc_local); 917 #if defined(PETSC_USE_COMPLEX) 918 alloc_local = fftw_mpi_local_size_3d(dim[0],dim[1],dim[2],comm,&local_n0,&local_0_start); 919 n = (PetscInt)local_n0*dim[1]*dim[2]; 920 ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr); 921 #else 922 printf("The code comes here\n"); 923 alloc_local = fftw_mpi_local_size_3d_transposed(dim[0],dim[1],dim[2]/2+1,PETSC_COMM_WORLD,&local_n0,&local_0_start,&local_n1,&local_1_start); 924 n = 2*(PetscInt)local_n0*dim[1]*(dim[2]/2+1); 925 ierr = MatSetSizes(A,n,n,2*dim[0]*dim[1]*(dim[2]/2+1),2*dim[0]*dim[1]*(dim[2]/2+1)); 926 #endif 927 break; 928 default: 929 alloc_local = fftw_mpi_local_size(ndim,pdim,comm,&local_n0,&local_0_start); 930 // printf("The value of alloc local is %ld from process %d\n",alloc_local,rank); 931 // alloc_local = fftw_mpi_local_size(ndim,dim,comm,&local_n0,&local_0_start); 932 n = (PetscInt)local_n0*partial_dim; 933 // printf("New partial dimension is %d %d %d",n,N,ndim); 934 ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr); 935 break; 936 } 937 } 938 ierr = PetscObjectChangeTypeName((PetscObject)A,MATFFTW);CHKERRQ(ierr); 939 940 ierr = PetscNewLog(A,Mat_FFTW,&fftw);CHKERRQ(ierr); 941 fft->data = (void*)fftw; 942 943 fft->n = n; 944 fftw->ndim_fftw = (ptrdiff_t)ndim; // This is dimension of fft 945 ierr = PetscMalloc(ndim*sizeof(ptrdiff_t), (ptrdiff_t *)&(fftw->dim_fftw));CHKERRQ(ierr); 946 for(ctr=0;ctr<ndim;ctr++) (fftw->dim_fftw)[ctr]=dim[ctr]; 947 948 fftw->p_forward = 0; 949 fftw->p_backward = 0; 950 fftw->p_flag = FFTW_ESTIMATE; 951 952 if (size == 1){ 953 A->ops->mult = MatMult_SeqFFTW; 954 A->ops->multtranspose = MatMultTranspose_SeqFFTW; 955 } else { 956 A->ops->mult = MatMult_MPIFFTW; 957 A->ops->multtranspose = MatMultTranspose_MPIFFTW; 958 } 959 fft->matdestroy = MatDestroy_FFTW; 960 A->ops->getvecs = MatGetVecs_FFTW; 961 A->assembled = PETSC_TRUE; 962 #if !defined(PETSC_USE_COMPLEX) 963 PetscObjectComposeFunctionDynamic((PetscObject)A,"InputTransformFFT_C","InputTransformFFT_FFTW",InputTransformFFT_FFTW); 964 PetscObjectComposeFunctionDynamic((PetscObject)A,"OutputTransformFFT_C","OutputTransformFFT_FFTW",OutputTransformFFT_FFTW); 965 #endif 966 967 /* get runtime options */ 968 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"FFTW Options","Mat");CHKERRQ(ierr); 969 ierr = PetscOptionsEList("-mat_fftw_plannerflags","Planner Flags","None",p_flags,4,p_flags[0],&p_flag,&flg);CHKERRQ(ierr); 970 if (flg) {fftw->p_flag = (unsigned)p_flag;} 971 PetscOptionsEnd(); 972 PetscFunctionReturn(0); 973 } 974 EXTERN_C_END 975 976 977 978 979