arrays.c (6768B)
1 #include <stdlib.h> 2 #include <stdio.h> 3 #include <math.h> 4 #include <err.h> 5 6 #define DEG2RAD(x) (x*M_PI/180.0) 7 8 void 9 check_magnitude(const char *func_name, int limit, int value) 10 { 11 if (value < limit) 12 errx(1, "%s: input size %d is less than %d\n", 13 func_name, value, limit); 14 } 15 16 /* Translate a i,j,k index in grid with dimensions nx, ny, nz into a 17 * linear index */ 18 unsigned int 19 idx3(const unsigned int i, 20 const unsigned int j, 21 const unsigned int k, 22 const unsigned int nx, 23 const unsigned int ny) 24 { 25 return i + nx * j + nx * ny * k; 26 } 27 28 /* Translate a i,j,k index in grid with dimensions nx, ny, nz and a 29 * padding of single ghost nodes into a linear index */ 30 unsigned int 31 idx3g(const unsigned int i, 32 const unsigned int j, 33 const unsigned int k, 34 const unsigned int nx, 35 const unsigned int ny) 36 { 37 return i + 1 + (nx + 2) * (j + 1) + (nx + 2) * (ny + 2) * (k + 1); 38 } 39 40 /* Translate a i,j index in grid with dimensions nx, ny into a linear 41 * index */ 42 unsigned int 43 idx2(const unsigned int i, const unsigned int j, const unsigned int nx) 44 { 45 return i + nx * j; 46 } 47 48 /* Translate a i,j index in grid with dimensions nx, ny and a padding of 49 * single ghost nodes into a linear index */ 50 unsigned int 51 idx2g(const unsigned int i, const unsigned int j, const unsigned int nx) 52 { 53 return i + 1 + (nx + 2) * (j + 1); 54 } 55 56 /* Translate a i index in grid with a padding of single into a linear 57 * index */ 58 unsigned int 59 idx1g(const unsigned int i) 60 { 61 return i + 1; 62 } 63 64 /* Return an array of `n` linearly spaced values in the range [lower; 65 * upper] */ 66 double * 67 linspace(const double lower, const double upper, const int n) 68 { 69 int i; 70 double *x; 71 double dx; 72 73 check_magnitude(__func__, 1, n); 74 if (!(x = calloc(n, sizeof(double)))) 75 err(1, "%s: x calloc", __func__); 76 dx = (upper - lower) / (double) (n - 1); 77 for (i = 0; i < n; ++i) 78 x[i] = lower + dx * i; 79 80 return x; 81 } 82 83 /* Return an array of `n-1` values with the intervals between `x` values */ 84 double * 85 spacing(const double *x, const int n) 86 { 87 int i; 88 double *dx; 89 90 check_magnitude(__func__, 2, n); 91 if (!(dx = calloc((n - 1), sizeof(double)))) 92 err(1, "%s: dx calloc", __func__); 93 for (i = 0; i < n - 1; ++i) 94 dx[i] = x[i + 1] - x[i]; 95 96 return dx; 97 } 98 99 /* Return an array of `n` values with the value 0.0 */ 100 double * 101 zeros(const int n) 102 { 103 double *x; 104 105 check_magnitude(__func__, 1, n); 106 if (!(x = calloc(n, sizeof(double)))) 107 err(1, "%s: x calloc", __func__); 108 109 return x; 110 } 111 112 /* Return an array of `n` values with the value 1.0 */ 113 double * 114 ones(const int n) 115 { 116 int i; 117 double *x; 118 119 check_magnitude(__func__, 1, n); 120 if (!(x = malloc(n * sizeof(double)))) 121 err(1, "%s: x malloc", __func__); 122 for (i = 0; i < n; ++i) 123 x[i] = 1.0; 124 125 return x; 126 } 127 128 /* Return an array of `n` values with a specified value */ 129 double * 130 initval(const double value, const int n) 131 { 132 int i; 133 double *x; 134 135 check_magnitude(__func__, 1, n); 136 if (!(x = malloc(n * sizeof(double)))) 137 err(1, "%s: x malloc", __func__); 138 for (i = 0; i < n; ++i) 139 x[i] = value; 140 141 return x; 142 } 143 144 /* Return an array of `n` uninitialized values */ 145 double * 146 empty(const int n) 147 { 148 double *out; 149 150 check_magnitude(__func__, 1, n); 151 152 if (!(out = malloc(n * sizeof(double)))) 153 err(1, "%s: malloc", __func__); 154 155 return out; 156 } 157 158 /* Return largest value in array `a` with size `n` */ 159 double 160 max(const double *a, const int n) 161 { 162 int i; 163 double maxval; 164 165 check_magnitude(__func__, 1, n); 166 maxval = -INFINITY; 167 for (i = 0; i < n; ++i) 168 if (a[i] > maxval) 169 maxval = a[i]; 170 171 return maxval; 172 } 173 174 /* Return maximum value in array, with early exit if threshold exceeded */ 175 double 176 max_with_threshold(const double *a, const int n, const double threshold) 177 { 178 int i; 179 double maxval; 180 181 check_magnitude(__func__, 1, n); 182 maxval = -INFINITY; 183 for (i = 0; i < n; ++i) { 184 if (a[i] > maxval) 185 maxval = a[i]; 186 if (maxval > threshold) 187 return maxval; 188 } 189 190 return maxval; 191 } 192 193 /* Return smallest value in array `a` with size `n` */ 194 double 195 min(const double *a, const int n) 196 { 197 int i; 198 double minval; 199 200 check_magnitude(__func__, 1, n); 201 minval = +INFINITY; 202 for (i = 0; i < n; ++i) 203 if (a[i] < minval) 204 minval = a[i]; 205 206 return minval; 207 } 208 209 void 210 print_array(const double *a, const int n) 211 { 212 int i; 213 214 check_magnitude(__func__, 1, n); 215 for (i = 0; i < n; ++i) 216 printf("%.17g\n", a[i]); 217 } 218 219 void 220 print_arrays(const double *a, const double *b, const int n) 221 { 222 int i; 223 224 check_magnitude(__func__, 1, n); 225 for (i = 0; i < n; ++i) 226 printf("%.17g\t%.17g\n", a[i], b[i]); 227 } 228 229 void 230 print_arrays_2nd_normalized(const double *a, const double *b, const int n) 231 { 232 int i; 233 double max_b; 234 235 check_magnitude(__func__, 1, n); 236 max_b = max(b, n); 237 for (i = 0; i < n; ++i) 238 printf("%.17g\t%.17g\n", a[i], b[i] / max_b); 239 } 240 241 void 242 print_three_arrays(const double *a, 243 const double *b, 244 const double *c, 245 const int n) 246 { 247 int i; 248 249 check_magnitude(__func__, 1, n); 250 for (i = 0; i < n; ++i) 251 printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); 252 } 253 254 void 255 fprint_arrays(FILE *fp, const double *a, const double *b, const int n) 256 { 257 int i; 258 259 check_magnitude(__func__, 1, n); 260 for (i = 0; i < n; ++i) 261 fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]); 262 } 263 264 void 265 fprint_three_arrays(FILE *fp, 266 const double *a, 267 const double *b, 268 const double *c, 269 const int n) 270 { 271 int i; 272 273 check_magnitude(__func__, 1, n); 274 for (i = 0; i < n; ++i) 275 fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); 276 } 277 278 void 279 copy_values(const double *in, double *out, const int n) 280 { 281 int i; 282 283 check_magnitude(__func__, 1, n); 284 for (i = 0; i < n; ++i) 285 out[i] = in[i]; 286 } 287 288 double * 289 copy(const double *in, const int n) 290 { 291 double *out; 292 293 check_magnitude(__func__, 1, n); 294 out = empty(n); 295 copy_values(in, out, n); 296 return out; 297 } 298 299 double * 300 normalize(const double *in, const int n) 301 { 302 int i; 303 double max_val; 304 double *out; 305 306 check_magnitude(__func__, 1, n); 307 out = empty(n); 308 copy_values(in, out, n); 309 max_val = max(out, n); 310 311 if (max_val == 0.0) 312 max_val = 1.0; 313 314 for (i = 0; i < n; ++i) 315 out[i] /= max_val; 316 317 return out; 318 } 319 320 double 321 euclidean_norm(const double *a, const int n) 322 { 323 int i; 324 double out = 0.0; 325 326 for (i = 0; i < n; i++) 327 out += a[i] * a[i]; 328 329 return sqrt(out); 330 } 331 332 double 333 euclidean_distance(const double *a, const double *b, const int n) 334 { 335 int i; 336 double out = 0.0; 337 338 for (i = 0; i < n; i++) 339 out += (b[i] - a[i]) * (b[i] - a[i]); 340 341 return sqrt(out); 342 } 343 344 double 345 dot(const double *a, const double *b, const int n) 346 { 347 int i; 348 double out = 0.0; 349 350 for (i = 0; i < n; i++) 351 out += a[i] * b[i]; 352 353 return out; 354 } 355 356 double * 357 cross(const double a[3], const double b[3]) 358 { 359 double *out = empty(3); 360 361 out[0] = a[1]*b[2] - a[2]*b[1]; 362 out[1] = a[2]*b[0] - a[0]*b[2]; 363 out[2] = a[0]*b[1] - a[1]*b[0]; 364 365 return out; 366 }