1 #include <omp.h>
    2 #include <stdio.h>
    3 #include <f2c.h>
    4 #include <string.h>
    5 
    6 #ifdef PARALLEL_SLOW
    7 dosub(s, p, g, cur, sh_s0, sh_s1, nxmax, nzmax, ip, itom, bnd1, sh_rd_glb, sh_array_glb)
    8   sh_ptr_real sh_s0, sh_s1;
    9   integer s, p, g, nxmax, nzmax, ip, itom;
   10   double bnd1;
   11   sh_ptr_rd_glb sh_rd_glb;
   12   sh_ptr_array_glb sh_array_glb;
   13 {
   14   integer i, j, ix, iz, n1, n2, m1, m2, image_dim1, s0_dim1, s1_dim1, count_dim1;
   15   integer t, c, k, l;
   16   rd_ptr_real s0, s1, image, count;
   17   rd_wr_ptr_real sub_s;
   18   real sold, soup, snew, cval, sval, xpr, zpr, b1, b2;
   19   rd_ptr_rd_glb rd_glb = sh_rd_glb;
   20   rd_ptr_array_glb array_glb = sh_array_glb;
   21   extern doublereal slow_();
   22 
   23   s0_dim1 = rd_glb->nv;
   24   s1_dim1 = nxmax;
   25   image_dim1 = nxmax;
   26   count_dim1 = nxmax;
   27   s0 = sh_s0;
   28   s1 = sh_s1;
   29   image = array_glb->cie[0]->image;
   30   count = array_glb->cie[0]->count;
   31   sub_s = array_glb->sub_s[p];
   32 
   33   /* 	calculate new slowness */
   34   b1 = (float)1. - bnd1 / (float)100.;
   35   b2 = bnd1 / (float)100. + (float)1.;
   36   iz = cur / nxmax;
   37   ix = cur % nxmax;
   38   for (c = 0; c < s; c++) {
   39         zpr = rd_glb->zmin + (real) iz * rd_glb->binz;
   40         xpr = rd_glb->xmin + (real) ix * rd_glb->binx;
   41 	if (ip == 1) {
   42 	  soup = slow_(s0, rd_glb->nv, rd_glb->nh,
   43 	      rd_glb->bx0, rd_glb->bz0, rd_glb->x0, rd_glb->z0,
   44 	      rd_glb->x1, rd_glb->z1, xpr, zpr, rd_glb->isp);
   45 	} else {
   46 	  soup = slow_(s1, nxmax, nzmax, rd_glb->binx,
   47 	      rd_glb->binz, rd_glb->xmin, rd_glb->zmin,
   48 	      rd_glb->xmax, rd_glb->zmax, xpr, zpr, rd_glb->isp);
   49 	}
   50 	n1 = max(ix - rd_glb->ipx,0);
   51 	n2 = min(ix + 1 + rd_glb->ipx,nxmax);
   52 	m1 = max(iz - rd_glb->ipz,0);
   53 	m2 = min(iz + 1 + rd_glb->ipz,nzmax);
   54 	cval = (float)0.;
   55 	sval = (float)0.;
   56 	for (j = m1; j < m2; ++j) {
   57 	  for (i = n1; i < n2; ++i) {
   58 	    cval += count[i + j * count_dim1];
   59 	    sval += image[i + j * image_dim1];
   60 	  }
   61 	}
   62 /*
   63 	cover[ix + iz * cover_dim1] = cval;
   64 */
   65 	if (cval == (float)0.) {
   66 	  snew = itom * soup;
   67 	} else {
   68 	  snew = sval / cval + itom * soup;
   69 	}
   70 	/*       apply % bounds if "relative" inversion */
   71 	/* 	and only if this is the last backprojection */
   72 	if (itom == 1 && ip == rd_glb->iback) {
   73 	  sold = slow_(s0, rd_glb->nv, rd_glb->nh,
   74 	      rd_glb->bx0, rd_glb->bz0, rd_glb->x0, rd_glb->z0,
   75 	      rd_glb->x1, rd_glb->z1, xpr, zpr, rd_glb->isp);
   76 	  /* 	precent bounds */
   77 	  if (snew > b2 * sold) {
   78 	    snew = b2 * sold;
   79 	  } else if (snew < b1 * sold) {
   80 	    snew = b1 * sold;
   81 	  }
   82 	  /* 	hard min/max bounds */
   83 	  if (snew > rd_glb->maxslow) {
   84 	    snew = rd_glb->maxslow;
   85 	  }
   86 	  if (snew < rd_glb->minslow) {
   87 	    snew = rd_glb->minslow;
   88 	  }
   89 	}
   90 	/* 	store new slowness */
   91         sub_s[c] = snew;
   92         ix++;
   93         if (ix == nxmax) {
   94           ix = 0;
   95           iz++;
   96         }
   97     }
   98 }
   99 
  100 doslow(sh_s0, sh_s1, nxmax, nzmax, ip, itom, bnd1, sh_rd_glb, sh_array_glb)
  101   sh_ptr_real sh_s0, sh_s1;
  102   integer nxmax, nzmax, ip, itom;
  103   sh_ptr_rd_glb sh_rd_glb;
  104   sh_ptr_array_glb sh_array_glb;
  105   double bnd1;
  106 {
  107   integer start, stop, s1_dim1, t, b, s, p, g, cur, c, iz, ix;
  108   rd_ptr_rd_glb rd_glb = sh_rd_glb;
  109   rd_ptr_array_glb array_glb = sh_array_glb;
  110   rd_wr_ptr_real s1;
  111   rd_ptr_real sub_s;
  112 
  113   get_time(&start);
  114   g = rd_glb->gran;
  115   t = (nzmax * nxmax);
  116   b = t % g;
  117   s = (t / g) + 1;
  118   cur = t;
  119   for (p = g - 1; p >= 0; p--) {
  120     if (b == 0) { s--; }
  121     b--;
  122     cur -= s;
  123     withonly {
  124       rd_wr(array_glb->sub_s[p]);
  125       rd(array_glb->cie[0]);
  126       rd(sh_s0);
  127       rd(sh_s1);
  128       rd(sh_array_glb);
  129       rd(sh_rd_glb);
  130     } do (s, p, g, cur, sh_s0, sh_s1, nxmax, nzmax, ip, itom, bnd1, sh_rd_glb, sh_array_glb) {
  131       int start, stop;
  132       get_time(&start);
  133       dosub(s, p, g, cur, sh_s0, sh_s1, nxmax, nzmax, ip, itom, bnd1, sh_rd_glb, sh_array_glb);
  134       get_time(&stop);
  135       printf("%d : dosub %lf\n", get_this_proc(),
  136         ticks_to_seconds_time(stop-start));
  137     }
  138   }
  139   t = (nzmax * nxmax);
  140   b = t % g;
  141   s = (t / g) + 1;
  142   cur = t;
  143   with { wr(sh_s1); } cont;
  144   s1_dim1 = nxmax;
  145   s1 = sh_s1;
  146   with {
  147     for (p = g - 1; p >= 0; p--) rd(array_glb->sub_s[p]);
  148   } cont;
  149   for (p = g - 1; p >= 0; p--) {
  150     if (b == 0) { s--; }
  151     b--;
  152     cur -= s;
  153     sub_s = array_glb->sub_s[p];
  154     iz = cur / nxmax;
  155     ix = cur % nxmax;
  156     for (c = 0; c < s; c++) {
  157       s1[ix + iz * s1_dim1] = sub_s[c];
  158       ix++;
  159       if (ix == nxmax) {
  160         ix = 0;
  161         iz++;
  162       }
  163     }
  164   }
  165   with {
  166     for (p = g - 1; p >= 0; p--) df_rd(array_glb->sub_s[p]);
  167   } cont;
  168   s1 = NULL;
  169   with { df_wr(sh_s1); } cont;
  170   get_time(&stop);
  171   printf("parallel slowness took %lf\n", ticks_to_seconds_time(stop-start));
  172 }
  173 #else
  174 
  175 doslow(sh_s0, sh_s1, nxmax, nzmax, ip, itom, bnd1, sh_rd_glb, sh_array_glb)
  176   sh_ptr_real sh_s0, sh_s1;
  177   integer nxmax, nzmax, ip, itom;
  178   sh_ptr_rd_glb sh_rd_glb;
  179   sh_ptr_array_glb sh_array_glb;
  180   double bnd1;
  181 {
  182   integer i, j, ix, iz, n1, n2, m1, m2, image_dim1, s0_dim1, s1_dim1, count_dim1;
  183   rd_ptr_real s0, s1, image, count;
  184   rd_wr_ptr_real rd_wr_s1;
  185   real sold, soup, snew, cval, sval, xpr, zpr, b1, b2;
  186   rd_ptr_rd_glb rd_glb = sh_rd_glb;
  187   rd_ptr_array_glb array_glb = sh_array_glb;
  188   extern doublereal slow_();
  189 
  190   s0_dim1 = rd_glb->nv;
  191   s1_dim1 = nxmax;
  192   image_dim1 = nxmax;
  193   count_dim1 = nxmax;
  194   s0 = sh_s0;
  195   s1 = sh_s1;
  196   rd_wr_s1 = sh_s1;
  197   image = array_glb->cie[0]->image;
  198   count = array_glb->cie[0]->count;
  199 
  200   /* 	calculate new slowness */
  201   b1 = (float)1. - bnd1 / (float)100.;
  202   b2 = bnd1 / (float)100. + (float)1.;
  203 
  204 
  205    #pragma omp parallel for private (ix, zpr, i, sold, j, sval, xpr, soup, snew, m1, m2, n1, n2, cval) schedule(static)
  206   for (iz = 0; iz < nzmax; ++iz) {
  207     zpr = rd_glb->zmin + (real) iz * rd_glb->binz;
  208     for (ix = 0; ix < nxmax; ++ix) {
  209       xpr = rd_glb->xmin + (real) ix * rd_glb->binx;
  210       if (ip == 1) {
  211         soup = slow_(s0, rd_glb->nv, rd_glb->nh,
  212           rd_glb->bx0, rd_glb->bz0, rd_glb->x0, rd_glb->z0,
  213           rd_glb->x1, rd_glb->z1, xpr, zpr, rd_glb->isp);
  214       } else {
  215         soup = slow_(s1, nxmax, nzmax, rd_glb->binx,
  216           rd_glb->binz, rd_glb->xmin, rd_glb->zmin,
  217           rd_glb->xmax, rd_glb->zmax, xpr, zpr, rd_glb->isp);
  218       }
  219       n1 = max(ix - rd_glb->ipx,0);
  220       n2 = min(ix + 1 + rd_glb->ipx,nxmax);
  221       m1 = max(iz - rd_glb->ipz,0);
  222       m2 = min(iz + 1 + rd_glb->ipz,nzmax);
  223       cval = (float)0.;
  224       sval = (float)0.;
  225       for (j = m1; j < m2; ++j) {
  226         for (i = n1; i < n2; ++i) {
  227           cval += count[i + j * count_dim1];
  228           sval += image[i + j * image_dim1];
  229         }
  230       }
  231 /*
  232       cover[ix + iz * cover_dim1] = cval;
  233 */
  234       if (cval == (float)0.) {
  235         snew = itom * soup;
  236       } else {
  237         snew = sval / cval + itom * soup;
  238       }
  239       /*       apply % bounds if "relative" inversion */
  240       /*     and only if this is the last backprojection */
  241       if (itom == 1 && ip == rd_glb->iback) {
  242         sold = slow_(s0, rd_glb->nv, rd_glb->nh,
  243             rd_glb->bx0, rd_glb->bz0, rd_glb->x0, rd_glb->z0,
  244             rd_glb->x1, rd_glb->z1, xpr, zpr, rd_glb->isp);
  245       /*     precent bounds */
  246         if (snew > b2 * sold) {
  247           snew = b2 * sold;
  248         } else if (snew < b1 * sold) {
  249           snew = b1 * sold;
  250         }
  251         /*     hard min/max bounds */
  252         if (snew > rd_glb->maxslow) {
  253           snew = rd_glb->maxslow;
  254         }
  255         if (snew < rd_glb->minslow) {
  256           snew = rd_glb->minslow;
  257         }
  258       }
  259     /*     store new slowness */
  260       rd_wr_s1[ix + iz * s1_dim1] = snew;
  261     }
  262   }
  263   rd_wr_s1 = NULL;
  264 }
  265 #endif