diff options
Diffstat (limited to 'src/Core/FiniteDifferenceLibrary.c')
| -rw-r--r-- | src/Core/FiniteDifferenceLibrary.c | 281 |
1 files changed, 81 insertions, 200 deletions
diff --git a/src/Core/FiniteDifferenceLibrary.c b/src/Core/FiniteDifferenceLibrary.c index f5736e2..fbf2646 100644 --- a/src/Core/FiniteDifferenceLibrary.c +++ b/src/Core/FiniteDifferenceLibrary.c @@ -2,161 +2,66 @@ #include "FiniteDifferenceLibrary.h" -DLL_EXPORT int openMPtest(void) +DLL_EXPORT int openMPtest(int nThreads) { - int nThreads; + omp_set_num_threads(nThreads); + + int nThreads_running; #pragma omp parallel { if (omp_get_thread_num() == 0) { - nThreads = omp_get_num_threads(); + nThreads_running = omp_get_num_threads(); } } - return nThreads; + return nThreads_running; } -int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc) +void threads_setup(int nThreads_requested, int *nThreads_current) { - size_t volume = nx * ny * nz; - - const float * inimage = inimagefull; - float * outimageX = outimageXfull; - float * outimageY = outimageYfull; - float * outimageZ = outimageZfull; - - int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice - int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row - - long c; - for (c = 0; c < nc; c++) - { #pragma omp parallel - { - long ind, k, j, i; - float pix0; - //run over all and then fix boundaries -#pragma omp for - for (ind = 0; ind < nx * ny * (nz - 1); ind++) - { - pix0 = -inimage[ind]; - - outimageX[ind] = pix0 + inimage[ind + 1]; - outimageY[ind] = pix0 + inimage[ind + nx]; - outimageZ[ind] = pix0 + inimage[ind + nx * ny]; - } - -#pragma omp for - for (ind = 0; ind < nx * (ny - 1); ind++) - { - pix0 = -inimage[ind + offset1]; - - outimageX[ind + offset1] = pix0 + inimage[ind + offset1 + 1]; - outimageY[ind + offset1] = pix0 + inimage[ind + offset1 + nx]; - } - -#pragma omp for - for (ind = 0; ind < nx - 1; ind++) - { - pix0 = -inimage[ind + offset2]; - - outimageX[ind + offset2] = pix0 + inimage[ind + offset2 + 1]; - } - - //boundaries -#pragma omp for - for (k = 0; k < nz; k++) - { - for (i = 0; i < nx; i++) - { - outimageY[(k * ny * nx) + (ny - 1) * nx + i] = 0; - } - } - -#pragma omp for - for (k = 0; k < nz; k++) - { - for (j = 0; j < ny; j++) - { - outimageX[k * ny * nx + j * nx + nx - 1] = 0; - } - } - - if (nz > 1) - { -#pragma omp for - for (ind = 0; ind < ny * nx; ind++) - { - outimageZ[nx * ny * (nz - 1) + ind] = 0; - } - } - } - - inimage += volume; - outimageX += volume; - outimageY += volume; - outimageZ += volume; - } - - - //now the rest of the channels - if (nc > 1) { - long ind; - - for (c = 0; c < nc - 1; c++) - { - float * outimageC = outimageCfull + c * volume; - const float * inimage = inimagefull + c * volume; - -#pragma omp parallel for - for (ind = 0; ind < volume; ind++) - { - outimageC[ind] = -inimage[ind] + inimage[ind + volume]; - } - } - -#pragma omp parallel for - for (ind = 0; ind < volume; ind++) + if (omp_get_thread_num() == 0) { - outimageCfull[(nc - 1) * volume + ind] = 0; + *nThreads_current = omp_get_num_threads(); } } - - return 0; + omp_set_num_threads(nThreads_requested); } -int fdiff_direct_neumann_L21(const float *inimagefull, float *outimageL21normfull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc) + +int fdiff_direct_neumann(const float *inimagefull, float *outimageXfull, float *outimageYfull, float *outimageZfull, float *outimageCfull, long nx, long ny, long nz, long nc) { size_t volume = nx * ny * nz; - int z_dim = nz - 1; - - const float * inimage = inimagefull; - float * outimageX = outimageXfull; - float * outimageY = outimageYfull; - float * outimageZ = outimageZfull; - float * outimageL21norm = outimageL21normfull; + const float *inimage = inimagefull; + float *outimageX = outimageXfull; + float *outimageY = outimageYfull; + float *outimageZ = outimageZfull; - int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice + int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row long c; + + int z_dim = nz > 1 ? 1: 0; + for (c = 0; c < nc; c++) { #pragma omp parallel { - long ind, k; + long ind, k, j, i; float pix0; //run over all and then fix boundaries -#pragma omp for + +#pragma omp for nowait for (ind = 0; ind < nx * ny * (nz - 1); ind++) { pix0 = -inimage[ind]; - outimageX[ind] = pix0 + inimage[ind + 1]; outimageY[ind] = pix0 + inimage[ind + nx]; outimageZ[ind] = pix0 + inimage[ind + nx * ny]; } -#pragma omp for +#pragma omp for nowait for (ind = 0; ind < nx * (ny - 1); ind++) { pix0 = -inimage[ind + offset1]; @@ -169,29 +74,26 @@ int fdiff_direct_neumann_L21(const float *inimagefull, float *outimageL21normful for (ind = 0; ind < nx - 1; ind++) { pix0 = -inimage[ind + offset2]; - outimageX[ind + offset2] = pix0 + inimage[ind + offset2 + 1]; } //boundaries -#pragma omp for +#pragma omp for nowait for (k = 0; k < nz; k++) { - for (int i = 0; i < nx; i++) + for (i = 0; i < nx; i++) { outimageY[(k * ny * nx) + (ny - 1) * nx + i] = 0; } } - -#pragma omp for +#pragma omp for nowait for (k = 0; k < nz; k++) { - for (int j = 0; j < ny; j++) + for (j = 0; j < ny; j++) { outimageX[k * ny * nx + j * nx + nx - 1] = 0; } } - if (z_dim) { #pragma omp for @@ -199,34 +101,15 @@ int fdiff_direct_neumann_L21(const float *inimagefull, float *outimageL21normful { outimageZ[nx * ny * (nz - 1) + ind] = 0; } - -#pragma omp for - for (ind = 0; ind < volume; ind++) - { - outimageL21norm[ind] = outimageX[ind] * outimageX[ind] + outimageY[ind] * outimageY[ind] + outimageZ[ind] * outimageZ[ind]; - } - - } - else - { - -#pragma omp for - for (ind = 0; ind < volume; ind++) - { - outimageL21norm[ind] = outimageX[ind] * outimageX[ind] + outimageY[ind] * outimageY[ind]; - } } } - inimage += volume; outimageX += volume; outimageY += volume; outimageZ += volume; - outimageL21norm += volume; } - //now the rest of the channels if (nc > 1) { @@ -234,18 +117,13 @@ int fdiff_direct_neumann_L21(const float *inimagefull, float *outimageL21normful for (c = 0; c < nc - 1; c++) { - float * outimageC = outimageCfull + c * volume; - float * outimageL21norm = outimageL21normfull + c * volume; - const float * inimage = inimagefull + c * volume; + float *outimageC = outimageCfull + c * volume; + const float *inimage = inimagefull + c * volume; #pragma omp parallel for for (ind = 0; ind < volume; ind++) { outimageC[ind] = -inimage[ind] + inimage[ind + volume]; - outimageL21norm[ind] += outimageC[ind] * outimageC[ind]; - - //sqrt - outimageL21norm[ind] = (float)sqrt(outimageL21norm[ind]); } } @@ -262,15 +140,14 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float { size_t volume = nx * ny * nz; - const float * inimage = inimagefull; - float * outimageX = outimageXfull; - float * outimageY = outimageYfull; - float * outimageZ = outimageZfull; + const float *inimage = inimagefull; + float *outimageX = outimageXfull; + float *outimageY = outimageYfull; + float *outimageZ = outimageZfull; - int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice + int offset1 = (nz - 1) * nx * ny; //ind to beginning of last slice int offset2 = offset1 + (ny - 1) * nx; //ind to beginning of last row - long c; for (c = 0; c < nc; c++) { @@ -280,7 +157,7 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float long ind, k; float pix0; //run over all and then fix boundaries -#pragma omp for +#pragma omp for nowait for (ind = 0; ind < nx * ny * (nz - 1); ind++) { pix0 = -inimage[ind]; @@ -290,7 +167,7 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float outimageZ[ind] = pix0 + inimage[ind + nx * ny]; } -#pragma omp for +#pragma omp for nowait for (ind = 0; ind < nx * (ny - 1); ind++) { pix0 = -inimage[ind + offset1]; @@ -308,7 +185,7 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float } //boundaries -#pragma omp for +#pragma omp for nowait for (k = 0; k < nz; k++) { for (int i = 0; i < nx; i++) @@ -320,7 +197,7 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float } } -#pragma omp for +#pragma omp for nowait for (k = 0; k < nz; k++) { for (int j = 0; j < ny; j++) @@ -332,10 +209,9 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float } } - if (nz > 1) { -#pragma omp for +#pragma omp for nowait for (ind = 0; ind < ny * nx; ind++) { outimageZ[nx * ny * (nz - 1) + ind] = -inimage[nx * ny * (nz - 1) + ind] + inimage[ind]; @@ -356,8 +232,8 @@ int fdiff_direct_periodic(const float *inimagefull, float *outimageXfull, float for (c = 0; c < nc - 1; c++) { - float * outimageC = outimageCfull + c * volume; - const float * inimage = inimagefull + c * volume; + float *outimageC = outimageCfull + c * volume; + const float *inimage = inimagefull + c * volume; #pragma omp parallel for for (ind = 0; ind < volume; ind++) @@ -383,18 +259,18 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const //assumes nx and ny > 1 int z_dim = nz - 1; - float * outimage = outimagefull; - const float * inimageX = inimageXfull; - const float * inimageY = inimageYfull; - const float * inimageZ = inimageZfull; + float *outimage = outimagefull; + const float *inimageX = inimageXfull; + const float *inimageY = inimageYfull; + const float *inimageZ = inimageZfull; - float * tempX = (float*)malloc(volume * sizeof(float)); - float * tempY = (float*)malloc(volume * sizeof(float)); - float * tempZ; - - if(z_dim) + float *tempX = (float *)malloc(volume * sizeof(float)); + float *tempY = (float *)malloc(volume * sizeof(float)); + float *tempZ; + + if (z_dim) { - tempZ = (float*)malloc(volume * sizeof(float)); + tempZ = (float *)malloc(volume * sizeof(float)); } long c; @@ -413,7 +289,6 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const for (ind = nx; ind < nx * ny * nz; ind++) { tempY[ind] = -inimageY[ind] + inimageY[ind - nx]; - } //boundaries @@ -442,7 +317,7 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const for (ind = nx * ny; ind < nx * ny * nz; ind++) { tempZ[ind] = -inimageZ[ind] + inimageZ[ind - nx * ny]; - } + } #pragma omp for for (ind = 0; ind < ny * nx; ind++) { @@ -454,7 +329,6 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const { outimage[ind] = tempX[ind] + tempY[ind] + tempZ[ind]; } - } else { @@ -464,7 +338,6 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const outimage[ind] = tempX[ind] + tempY[ind]; } } - } outimage += volume; @@ -475,10 +348,9 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const free(tempX); free(tempY); - if(z_dim) + if (z_dim) free(tempZ); - // //now the rest of the channels if (nc > 1) { @@ -500,7 +372,6 @@ int fdiff_adjoint_neumann(float *outimagefull, const float *inimageXfull, const outimagefull[ind] += -inimageCfull[ind]; outimagefull[(nc - 1) * volume + ind] += inimageCfull[(nc - 2) * volume + ind]; } - } return 0; @@ -513,18 +384,18 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const //assumes nx and ny > 1 int z_dim = nz - 1; - float * outimage = outimagefull; - const float * inimageX = inimageXfull; - const float * inimageY = inimageYfull; - const float * inimageZ = inimageZfull; + float *outimage = outimagefull; + const float *inimageX = inimageXfull; + const float *inimageY = inimageYfull; + const float *inimageZ = inimageZfull; - float * tempX = (float*)malloc(volume * sizeof(float)); - float * tempY = (float*)malloc(volume * sizeof(float)); - float * tempZ; + float *tempX = (float *)malloc(volume * sizeof(float)); + float *tempY = (float *)malloc(volume * sizeof(float)); + float *tempZ; if (z_dim) { - tempZ = (float*)malloc(volume * sizeof(float)); + tempZ = (float *)malloc(volume * sizeof(float)); } long c; @@ -583,7 +454,6 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const { outimage[ind] = tempX[ind] + tempY[ind] + tempZ[ind]; } - } else { @@ -593,7 +463,6 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const outimage[ind] = tempX[ind] + tempY[ind]; } } - } outimage += volume; @@ -604,7 +473,7 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const free(tempX); free(tempY); - if(z_dim) + if (z_dim) free(tempZ); //now the rest of the channels @@ -632,9 +501,11 @@ int fdiff_adjoint_periodic(float *outimagefull, const float *inimageXfull, const return 0; } - -DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, float *gradYfull, float *gradXfull, long nc, long nz, long ny, long nx, int boundary, int direction) +DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, float *gradYfull, float *gradXfull, long nc, long nz, long ny, long nx, int boundary, int direction, int nThreads) { + int nThreads_initial; + threads_setup(nThreads, &nThreads_initial); + if (boundary) { if (direction) @@ -643,17 +514,21 @@ DLL_EXPORT int fdiff4D(float *imagefull, float *gradCfull, float *gradZfull, flo fdiff_adjoint_periodic(imagefull, gradXfull, gradYfull, gradZfull, gradCfull, nx, ny, nz, nc); } else - { + { if (direction) fdiff_direct_neumann(imagefull, gradXfull, gradYfull, gradZfull, gradCfull, nx, ny, nz, nc); else fdiff_adjoint_neumann(imagefull, gradXfull, gradYfull, gradZfull, gradCfull, nx, ny, nz, nc); } - + + omp_set_num_threads(nThreads_initial); return 0; } -DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, float *gradXfull, long nz, long ny, long nx, int boundary, int direction) +DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, float *gradXfull, long nz, long ny, long nx, int boundary, int direction, int nThreads) { + int nThreads_initial; + threads_setup(nThreads, &nThreads_initial); + if (boundary) { if (direction) @@ -669,10 +544,14 @@ DLL_EXPORT int fdiff3D(float *imagefull, float *gradZfull, float *gradYfull, flo fdiff_adjoint_neumann(imagefull, gradXfull, gradYfull, gradZfull, NULL, nx, ny, nz, 1); } + omp_set_num_threads(nThreads_initial); return 0; } -DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, long ny, long nx, int boundary, int direction) +DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, long ny, long nx, int boundary, int direction, int nThreads) { + int nThreads_initial; + threads_setup(nThreads, &nThreads_initial); + if (boundary) { if (direction) @@ -688,5 +567,7 @@ DLL_EXPORT int fdiff2D(float *imagefull, float *gradYfull, float *gradXfull, lon fdiff_adjoint_neumann(imagefull, gradXfull, gradYfull, NULL, NULL, nx, ny, 1, 1); } + omp_set_num_threads(nThreads_initial); return 0; } + |