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 for (iz = 0; iz < nzmax; ++iz) { 206 zpr = rd_glb->zmin + (real) iz * rd_glb->binz; 207 #pragma omp parallel for private (i, sold, j, sval, xpr, soup, snew, m1, m2, n1, n2, cval) schedule(static) 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