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