00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <mpi.h>
00004 #include <math.h>
00005 #include <string.h>
00006
00007 int setpart_gc1(MPI_Comm comm_2d, int nangs, int *psize, int *nbase)
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 {
00019 int ROW = 0, COL = 1;
00020 int dims[2], periods[2], mycoords[2];
00021 int nangsloc, nrem;
00022
00023
00024 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00025
00026 nangsloc = nangs/dims[ROW];
00027 nrem = nangs - dims[ROW]*nangsloc;
00028 if (mycoords[ROW] < nrem) nangsloc++;
00029
00030 for (int i = 0; i < dims[ROW]; i++) {
00031 psize[i] = nangs/dims[ROW];
00032 if (i < nrem) psize[i]++;
00033 }
00034
00035 nbase[0] = 0;
00036 for (int i = 1; i < dims[ROW]; i++) {
00037 nbase[i] = nbase[i-1] + psize[i-1];
00038 }
00039
00040 return nangsloc;
00041 }
00042
00043 int getnnz(int *volsize, int ri, int *origin, int *nrays, int *nnz)
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 {
00058 int ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
00059 int ftm=0, status = 0;
00060
00061 r2 = ri*ri;
00062 *nnz = 0;
00063 *nrays = 0;
00064 int nx = volsize[0];
00065 int ny = volsize[1];
00066 int nz = volsize[2];
00067
00068
00069
00070
00071 int xcent = origin[0];
00072 int ycent = origin[1];
00073 int zcent = origin[2];
00074
00075
00076
00077
00078
00079 for (ix = 1; ix <=nx; ix++) {
00080 xs = ix-xcent;
00081 xx = xs*xs;
00082 for (iy = 1; iy <= ny; iy++) {
00083 ys = iy-ycent;
00084 yy = ys*ys;
00085 ftm = 1;
00086 for (iz = 1; iz <= nz; iz++) {
00087 zs = iz-zcent;
00088 zz = zs*zs;
00089 rs = xx + yy + zz;
00090 if (rs <= r2) {
00091 (*nnz)++;
00092 if (ftm) {
00093 (*nrays)++;
00094 ftm = 0;
00095 }
00096 }
00097 }
00098 }
00099 }
00100 return status;
00101 }
00102
00103
00104 int setpart_gr1(MPI_Comm comm_2d, int nnz, int *nnzpart, int *nnzbase)
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115 {
00116 int ROW = 0, COL = 1;
00117 int dims[2], periods[2], mycoords[2];
00118 int nnzloc, nrem;
00119
00120
00121 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00122
00123 nnzloc = nnz/dims[COL];
00124 nrem = nnz - dims[COL]*nnzloc;
00125 if (mycoords[COL] < nrem) nnzloc++;
00126
00127 for (int i = 0; i < dims[COL]; i++) {
00128 nnzpart[i] = nnz/dims[COL];
00129 if (i < nrem) nnzpart[i]++;
00130 }
00131
00132 nnzbase[0] = 0;
00133 for (int i = 1; i < dims[COL]; i++) {
00134 nnzbase[i] = nnzbase[i-1] + nnzpart[i-1];
00135 }
00136 nnzbase[dims[COL]] = nnz;
00137 return nnzloc;
00138 }
00139
00140 #define cord(i,j) cord[((j)-1)*3 + (i) -1]
00141 #define ptrs(i) ptrs[(i)-1]
00142 #define dm(i) dm[(i)-1]
00143
00144 int getcb2sph(int *volsize, int ri, int *origin, int nnz0, int *ptrs, int *cord)
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158 {
00159 int xs, ys, zs, xx, yy, zz, rs, r2;
00160 int ix, iy, iz, jnz, nnz, nrays;
00161 int ftm = 0, status = 0;
00162
00163 int xcent = (int)origin[0];
00164 int ycent = (int)origin[1];
00165 int zcent = (int)origin[2];
00166
00167 int nx = (int)volsize[0];
00168 int ny = (int)volsize[1];
00169 int nz = (int)volsize[2];
00170
00171 r2 = ri*ri;
00172 nnz = 0;
00173 nrays = 0;
00174 ptrs(1) = 1;
00175
00176 for (ix = 1; ix <= nx; ix++) {
00177 xs = ix-xcent;
00178 xx = xs*xs;
00179 for ( iy = 1; iy <= ny; iy++ ) {
00180 ys = iy-ycent;
00181 yy = ys*ys;
00182 jnz = 0;
00183
00184 ftm = 1;
00185
00186 for (iz = 1; iz <= nz; iz++) {
00187 zs = iz-zcent;
00188 zz = zs*zs;
00189 rs = xx + yy + zz;
00190 if (rs <= r2) {
00191 jnz++;
00192 nnz++;
00193
00194
00195
00196 if (ftm) {
00197 nrays++;
00198 cord(1,nrays) = iz;
00199 cord(2,nrays) = iy;
00200 cord(3,nrays) = ix;
00201 ftm = 0;
00202 }
00203 }
00204 }
00205 if (jnz > 0) {
00206 ptrs(nrays+1) = ptrs(nrays) + jnz;
00207 }
00208 }
00209 }
00210 if (nnz != nnz0) status = -1;
00211 return status;
00212 }
00213
00214 int sphpart(MPI_Comm comm_2d, int nrays, int *ptrs, int *nnzbase, int *ptrstart)
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226 {
00227 int ROW = 0, COL = 1;
00228 int dims[2], periods[2], mycoords[2];
00229 int nraysloc = 0;
00230
00231
00232 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00233
00234 int count = 1;
00235 if (mycoords[COL] == 0){
00236 nraysloc++;
00237 }
00238 ptrstart[0] = 0;
00239
00240 for(int i=1; i<nrays ; i++){
00241 if (ptrs[i]-1 <= nnzbase[count] && ptrs[i+1]-1 >= nnzbase[count]){
00242
00243
00244 if (nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]){
00245 if(mycoords[COL] == count-1){
00246 nraysloc++;
00247 }
00248 ptrstart[count] = i+1;
00249 count++;
00250
00251 } else {
00252 if(mycoords[COL] == count){
00253 nraysloc++;
00254 }
00255 ptrstart[count] = i;
00256 count++;
00257 }
00258
00259 }
00260 else {
00261 if(mycoords[COL] == count-1){
00262 nraysloc++;
00263 }
00264
00265 }
00266 }
00267 ptrstart[dims[COL]] = nrays;
00268 return nraysloc;
00269
00270 }
00271
00272 int make_proj_mat(float phi, float theta, float psi, float * dm)
00273
00274
00275
00276
00277
00278
00279 {
00280
00281
00282
00283
00284
00285
00286
00287 double cphi, sphi, cthe, sthe, cpsi, spsi;
00288 double dphi = phi;
00289 double dthe = theta;
00290 double dpsi = psi;
00291 sincos(dphi, &sphi, &cphi);
00292 sincos(dthe, &sthe, &cthe);
00293 sincos(dpsi, &spsi, &cpsi);
00294
00295 dm[0]=cphi*cthe*cpsi-sphi*spsi;
00296 dm[1]=sphi*cthe*cpsi+cphi*spsi;
00297 dm[2]=-sthe*cpsi;
00298 dm[3]=-cphi*cthe*spsi-sphi*cpsi;
00299 dm[4]=-sphi*cthe*spsi+cphi*cpsi;
00300 dm[5]=sthe*spsi;
00301
00302 return 0;
00303 }
00304
00305 int ifix(float a)
00306
00307
00308 {
00309 int ia = 0;
00310 if (a >= 0.0) {
00311 ia = (int)floor(a);
00312 }
00313 else {
00314 ia = (int)ceil(a);
00315 }
00316 return ia;
00317 }
00318
00319 #define x(i) x[(i)-1]
00320 #define y(i,j) y[((j)-1)*nx + (i) - 1]
00321
00322 int fwdpj3_Cart(int *volsize, int nraysloc, int nnzloc, float *dm, int *origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337 {
00338 int iqx, iqy, i, j, xc, yc, zc, jj;
00339 float ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
00340 int status = 0;
00341
00342
00343 float sx, sy;
00344
00345 sx = dm(7);
00346 sy = dm(8);
00347
00348 int xcent = origin[0];
00349 int ycent = origin[1];
00350 int zcent = origin[2];
00351
00352 int nx = volsize[0];
00353 int ny = volsize[1];
00354
00355 dm1 = dm(1);
00356 dm4 = dm(4);
00357
00358 if ( nx > 2*ri ) {
00359 for (i = 1; i <= nraysloc; i++) {
00360
00361 zc = cord(1,myptrstart+i)-zcent;
00362 yc = cord(2,myptrstart+i)-ycent;
00363 xc = cord(3,myptrstart+i)-xcent;
00364 xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
00365 yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;
00366
00367 for (j = ptrs(myptrstart+i); j< ptrs(myptrstart+i+1); j++) {
00368 jj = j-ptrs[myptrstart]+1;
00369 iqx = ifix(xb);
00370 iqy = ifix(yb);
00371
00372
00373 ct = x(jj);
00374
00375
00376
00377 dipx = xb - iqx;
00378 dipy = (yb - iqy) * ct;
00379
00380 dipy1m = ct - dipy;
00381 dipx1m = 1.0 - dipx;
00382
00383 if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1)
00384
00385 y(iqx ,iqy) += dipx1m*dipy1m;
00386 if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1)
00387
00388 y(iqx+1,iqy) += dipx*dipy1m;
00389 if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0)
00390
00391 y(iqx+1,iqy+1) += dipx*dipy;
00392 if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0)
00393
00394 y(iqx ,iqy+1) += dipx1m*dipy;
00395 xb += dm1;
00396 yb += dm4;
00397 }
00398 }
00399 }
00400 else {
00401 fprintf(stderr, " nx must be greater than 2*ri\n");
00402 exit(1);
00403 }
00404 return status;
00405 }
00406 #undef x
00407 #undef y
00408
00409 #define y(i) y[(i)-1]
00410 #define x(i,j) x[((j)-1)*nx + (i) - 1]
00411
00412 int bckpj3_Cart(int *volsize, int nraysloc, int nnzloc, float *dm, int *origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427 {
00428 int i, j, iqx,iqy, xc, yc, zc, jj;
00429 float xb, yb, dx, dy, dx1m, dy1m, dxdy;
00430 int status = 0;
00431
00432 int xcent = origin[0];
00433 int ycent = origin[1];
00434 int zcent = origin[2];
00435
00436 int nx = volsize[0];
00437 int ny = volsize[1];
00438
00439
00440 float sx, sy;
00441
00442 sx = dm(7);
00443 sy = dm(8);
00444
00445 if ( nx > 2*ri) {
00446 for (i = 1; i <= nraysloc; i++) {
00447 zc = cord(1,myptrstart+i) - zcent;
00448 yc = cord(2,myptrstart+i) - ycent;
00449 xc = cord(3,myptrstart+i) - xcent;
00450
00451 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
00452 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;
00453
00454 for (j = ptrs(myptrstart+i); j <ptrs(myptrstart+i+1); j++) {
00455 jj = j-ptrs[myptrstart]+1;
00456
00457
00458 iqx = ifix(xb);
00459 iqy = ifix(yb);
00460
00461 dx = xb - iqx;
00462 dy = yb - iqy;
00463 dx1m = 1.0 - dx;
00464 dy1m = 1.0 - dy;
00465 dxdy = dx*dy;
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485 if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
00486 y(jj) += x(iqx,iqy);
00487 if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
00488 y(jj) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
00489 }
00490 if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
00491 y(jj) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
00492 }
00493 if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
00494 y(jj) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00495 }
00496 }
00497
00498
00499
00500
00501
00502
00503
00504 xb += dm(1);
00505 yb += dm(4);
00506 }
00507 }
00508 }
00509 else {
00510 fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00511 }
00512
00513 return status;
00514 }
00515
00516 #undef x
00517 #undef y
00518
00519
00520
00521
00522 #define x(i) x[(i)-1]
00523 #define y(i,j) y[((j)-1)*nx + (i) - 1]
00524
00525
00526 int fwdpj3(int *volsize, int nrays, int nnz, float *dm, int *origin, int ri, int *ptrs, int *cord, float *x, float *y)
00527 {
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544 int iqx, iqy, i, j, xc, yc, zc;
00545 float ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
00546 int status = 0;
00547
00548
00549 float sx, sy;
00550
00551 sx = dm(7);
00552 sy = dm(8);
00553
00554 int xcent = origin[0];
00555 int ycent = origin[1];
00556 int zcent = origin[2];
00557
00558 int nx = volsize[0];
00559 int ny = volsize[1];
00560
00561 dm1 = dm(1);
00562 dm4 = dm(4);
00563
00564 if ( nx > 2*ri ) {
00565 for (i = 1; i <= nrays; i++) {
00566
00567 zc = cord(1,i)-zcent;
00568 yc = cord(2,i)-ycent;
00569 xc = cord(3,i)-xcent;
00570 xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
00571 yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;
00572
00573 for (j = ptrs(i); j< ptrs(i+1); j++) {
00574 iqx = ifix(xb);
00575 iqy = ifix(yb);
00576
00577 ct = x(j);
00578
00579
00580
00581 dipx = xb - iqx;
00582 dipy = (yb - iqy) * ct;
00583
00584 dipy1m = ct - dipy;
00585 dipx1m = 1.0 - dipx;
00586
00587 if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1)
00588
00589 y(iqx ,iqy) += dipx1m*dipy1m;
00590 if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1)
00591
00592 y(iqx+1,iqy) += dipx*dipy1m;
00593 if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0)
00594
00595 y(iqx+1,iqy+1) += dipx*dipy;
00596 if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0)
00597
00598 y(iqx ,iqy+1) += dipx1m*dipy;
00599 xb += dm1;
00600 yb += dm4;
00601 }
00602 }
00603 }
00604 else {
00605 fprintf(stderr, " nx must be greater than 2*ri\n");
00606 exit(1);
00607 }
00608 return status;
00609 }
00610 #undef x
00611 #undef y
00612
00613
00614 #define y(i) y[(i)-1]
00615 #define x(i,j) x[((j)-1)*nx + (i) - 1]
00616
00617
00618 int bckpj3(int *volsize, int nrays, int nnz, float *dm, int *origin, int ri, int *ptrs, int *cord, float *x, float *y)
00619 {
00620 int i, j, iqx,iqy, xc, yc, zc;
00621 float xb, yb, dx, dy, dx1m, dy1m, dxdy;
00622 int status = 0;
00623
00624 int xcent = origin[0];
00625 int ycent = origin[1];
00626 int zcent = origin[2];
00627
00628 int nx = volsize[0];
00629 int ny = volsize[1];
00630
00631
00632 float sx, sy;
00633
00634 sx = dm(7);
00635 sy = dm(8);
00636
00637
00638 if ( nx > 2*ri) {
00639 for (i = 1; i <= nrays; i++) {
00640 zc = cord(1,i) - zcent;
00641 yc = cord(2,i) - ycent;
00642 xc = cord(3,i) - xcent;
00643
00644 xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
00645 yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;
00646
00647 for (j = ptrs(i); j <ptrs(i+1); j++) {
00648 iqx = ifix(xb);
00649 iqy = ifix(yb);
00650
00651 dx = xb - iqx;
00652 dy = yb - iqy;
00653 dx1m = 1.0 - dx;
00654 dy1m = 1.0 - dy;
00655 dxdy = dx*dy;
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675 if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
00676 y(j) += x(iqx,iqy);
00677 if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
00678 y(j) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
00679 }
00680 if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
00681 y(j) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
00682 }
00683 if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
00684 y(j) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00685 }
00686 }
00687
00688
00689
00690
00691
00692
00693
00694 xb += dm(1);
00695 yb += dm(4);
00696 }
00697 }
00698 }
00699 else {
00700 fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00701 }
00702
00703 return status;
00704 }
00705
00706 #undef x
00707 #undef y
00708
00709 #undef dm
00710