util.c 12.9 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
/*-----------------------------------------------------------------------------------------
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
Tilman Steinweg's avatar
Tilman Steinweg committed
3
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
4
 * This file is part of IFOS.
5
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
6
 * IFOS is free software: you can redistribute it and/or modify
Tilman Steinweg's avatar
Tilman Steinweg committed
7 8
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
9
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
10
 * IFOS is distributed in the hope that it will be useful,
Tilman Steinweg's avatar
Tilman Steinweg committed
11 12 13
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
14
 *
Tilman Steinweg's avatar
Tilman Steinweg committed
15
 * You should have received a copy of the GNU General Public License
Florian Wittkamp's avatar
Florian Wittkamp committed
16
 * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
17
 -----------------------------------------------------------------------------------------*/
Tilman Steinweg's avatar
Tilman Steinweg committed
18 19 20 21 22 23 24 25


#define SHIFT_IND 1
#define FREE_ARGUMENT char *

#include "fd.h"


26
void declare_error(char err_text[]){
27 28 29 30 31 32 33 34 35
    extern int MYID;
    
    fprintf(stdout,"Message from PE %d\n",MYID);
    fprintf(stdout,"R U N - T I M E  E R R O R: \n");
    fprintf(stdout,"%s\n",err_text);
    fprintf(stdout,"...now exiting to system.\n");
    
    MPI_Abort(MPI_COMM_WORLD, 1);
    exit(1);
Tilman Steinweg's avatar
Tilman Steinweg committed
36 37 38
}

void warning(char warn_text[]){
39 40 41
    /* standard warnings handler */
    fprintf(stdout,"W A R N I N G   M E S S A G E: \n");
    fprintf(stdout,"%s\n",warn_text);
Tilman Steinweg's avatar
Tilman Steinweg committed
42 43 44 45
}


double maximum(float **a, int nx, int ny){
46 47 48 49 50 51 52 53
    double maxi=0.0;
    int i, j;
    
    
    for (j=1;j<=ny;j++)
        for (i=1;i<=nx;i++)
            if (fabs((double) a[i][j])>maxi) maxi=fabs((double)a[i][j]);
    return maxi;
Tilman Steinweg's avatar
Tilman Steinweg committed
54 55 56 57
}

float minimum_m(float **mat, int nx, int ny)
{
58 59 60 61 62 63 64 65 66 67 68 69 70 71
    float minm;
    int i,j;
    
    minm = mat[1][1];
    for (i=1;i<=nx;i++)
        for (j=1;j<=ny;j++)
        {
            if (i*j==1) continue;
            if (mat[j][i] < minm)
            {
                minm = mat[j][i];
            }
        }
    return minm;
Tilman Steinweg's avatar
Tilman Steinweg committed
72 73 74 75
}

float maximum_m(float **mat, int nx, int ny)
{
76 77 78 79 80 81 82 83 84 85 86 87 88 89
    float maxm;
    int i,j;
    
    maxm = mat[1][1];
    for (i=1;i<=nx;i++)
        for (j=1;j<=ny;j++)
        {
            if (i*j==1) continue;
            if (mat[j][i] > maxm)
            {
                maxm = mat[j][i];
            }
        }
    return maxm;
Tilman Steinweg's avatar
Tilman Steinweg committed
90 91 92 93 94 95
}




float *vector(int ni, int nj){
96 97 98 99
    float *a;
    int k;
    
    a=(float *)malloc((size_t) ((nj-ni+1+SHIFT_IND)*sizeof(float)));
100
    if (!a) declare_error("util.c: allocation failure in function vector()");
101 102
    for (k=0;k<(nj-ni+1+SHIFT_IND);k++) a[k]=0.0;
    return a-ni+SHIFT_IND;
Tilman Steinweg's avatar
Tilman Steinweg committed
103 104 105 106
}


int *ivector(int ni, int nj){
107 108 109 110 111
    
    int *a;
    int k;
    
    a=(int *)malloc((size_t) ((nj-ni+1+SHIFT_IND)*sizeof(int)));
112
    if (!a) declare_error("util.c: allocation failure in function ivector()");
113 114
    for (k=0;k<(nj-ni+1+SHIFT_IND);k++) a[k]=0;
    return a-ni+SHIFT_IND;
Tilman Steinweg's avatar
Tilman Steinweg committed
115 116 117
}

unsigned short int *usvector(int ni, int nj){
118 119 120 121 122
    
    unsigned short int *a;
    int k;
    
    a=(unsigned short int *)malloc((size_t) ((nj-ni+1+SHIFT_IND)*sizeof(unsigned short int)));
123
    if (!a) declare_error("util.c: allocation failure in function usvector()");
124 125
    for (k=0;k<(nj-ni+1+SHIFT_IND);k++) a[k]=0;
    return a-ni+SHIFT_IND;
Tilman Steinweg's avatar
Tilman Steinweg committed
126 127 128 129
}


unsigned char *cvector(int ni, int nj){
130 131 132 133
    
    unsigned char *a;
    
    a=(unsigned char *)malloc((size_t) ((nj-ni+1+SHIFT_IND)*sizeof(unsigned char)));
134
    if (!a) declare_error("util.c: allocation failure in function cvector()");
135
    return a-ni+SHIFT_IND;
Tilman Steinweg's avatar
Tilman Steinweg committed
136 137 138 139
}


unsigned long *lvector(int ni, int nj){
140 141 142 143 144
    
    unsigned long *a;
    int k;
    
    a=(unsigned long *)malloc((size_t) ((nj-ni+1+SHIFT_IND)*sizeof(unsigned long)));
145
    if (!a) declare_error("util.c: allocation failure in function lvector()");
146 147
    for (k=0;k<(nj-ni+1+SHIFT_IND);k++) a[k]=0;
    return a-ni+SHIFT_IND;
Tilman Steinweg's avatar
Tilman Steinweg committed
148 149 150
}

double *dvector(int ni, int nj){
151 152 153 154 155
    
    double *a;
    int k;
    
    a=(double *)malloc((size_t) ((nj-ni+1+SHIFT_IND)*sizeof(double)));
156
    if (!a) declare_error("util.c: allocation failure in function dvector()");
157 158
    for (k=0;k<(nj-ni+1+SHIFT_IND);k++) a[k]=0.0;
    return a-ni+SHIFT_IND;
Tilman Steinweg's avatar
Tilman Steinweg committed
159 160 161
}

float **fmatrix(int mrl, int mrh, int mcl, int mch){
162 163 164 165 166
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    float **ma;
    
    ma=(float **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(float*)));
167
    if (!ma) declare_error("util.c: allocation failure 1 in function fmatrix() ");
168 169 170 171
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(float *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(float)));
172
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function fmatrix() ");
173 174 175 176 177 178 179 180 181
    ma[mrl] += SHIFT_IND;
    ma[mrl] -= mcl;
    
    for (k=mrl+1;k<=mrh;k++) ma[k]=ma[k-1]+mcol;
    
    for (k=mrl;k<=mrh;k++)
        for (l=mcl;l<=mch;l++) ma[k][l]=0.0;
    
    return ma;
Tilman Steinweg's avatar
Tilman Steinweg committed
182 183 184
}

float **matrix(int mrl, int mrh, int mcl, int mch){
185 186 187 188 189
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    float **ma;
    
    ma=(float **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(float*)));
190
    if (!ma) declare_error("util.c: allocation failure 1 in function matrix() ");
191 192 193 194
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(float *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(float)));
195
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function matrix() ");
196 197 198 199 200 201 202 203 204
    ma[mrl] += SHIFT_IND;
    ma[mrl] -= mcl;
    
    for (k=mrl+1;k<=mrh;k++) ma[k]=ma[k-1]+mcol;
    
    for (k=mrl;k<=mrh;k++)
        for (l=mcl;l<=mch;l++) ma[k][l]=0.0;
    
    return ma;
Tilman Steinweg's avatar
Tilman Steinweg committed
205 206 207 208
}


double **dmatrix(int mrl, int mrh, int mcl, int mch){
209 210 211 212 213
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    double **ma;
    
    ma=(double **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(double*)));
214
    if (!ma) declare_error("util.c: allocation failure 1 in function matrix() ");
215 216 217 218
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(double *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(double)));
219
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function dmatrix() ");
220 221 222 223 224 225 226 227 228
    ma[mrl] += SHIFT_IND;
    ma[mrl] -= mcl;
    
    for (k=mrl+1;k<=mrh;k++) ma[k]=ma[k-1]+mcol;
    
    for (k=mrl;k<=mrh;k++)
        for (l=mcl;l<=mch;l++) ma[k][l]=0.0;
    
    return ma;
Tilman Steinweg's avatar
Tilman Steinweg committed
229 230 231
}

int **imatrix(int mrl, int mrh, int mcl, int mch){
232 233 234 235 236
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    int **ma;
    
    ma=(int **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(int*)));
237
    if (!ma) declare_error("util.c: allocation failure 1 in function imatrix() ");
238 239 240 241
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(int *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(int)));
242
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function imatrix() ");
243 244 245 246 247 248 249 250 251
    ma[mrl] += SHIFT_IND;
    ma[mrl] -= mcl;
    
    for (k=mrl+1;k<=mrh;k++) ma[k]=ma[k-1]+mcol;
    
    for (k=mrl;k<=mrh;k++)
        for (l=mcl;l<=mch;l++) ma[k][l]=0;
    
    return ma;
Tilman Steinweg's avatar
Tilman Steinweg committed
252 253 254
}

unsigned short int **usmatrix(int mrl, int mrh, int mcl, int mch){
255 256 257 258 259
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    unsigned short int **ma;
    
    ma=(unsigned short int **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(unsigned short int*)));
260
    if (!ma) declare_error("util.c: allocation failure 1 in function usmatrix() ");
261 262 263 264
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(unsigned short int *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(unsigned short int)));
265
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function usmatrix() ");
266 267 268 269 270 271 272 273 274
    ma[mrl] += SHIFT_IND;
    ma[mrl] -= mcl;
    
    for (k=mrl+1;k<=mrh;k++) ma[k]=ma[k-1]+mcol;
    
    for (k=mrl;k<=mrh;k++)
        for (l=mcl;l<=mch;l++) ma[k][l]=0;
    
    return ma;
Tilman Steinweg's avatar
Tilman Steinweg committed
275 276 277
}

float ***f3tensor(int mrl, int mrh, int mcl, int mch,int mdl, int mdh){
278 279 280 281 282
    
    int w,e,r, mrow=mrh-mrl+1,mcol=mch-mcl+1,mdep=mdh-mdl+1;
    float ***te;
    
    te=(float ***) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(float**)));
283
    if (!te) declare_error("util.c: allocation failure 1 in function f3tensor() ");
284 285 286 287
    te += SHIFT_IND;
    te -= mrl;
    
    te[mrl]=(float **) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(float*)));
288
    if (!te[mrl]) declare_error("util.c: allocation failure 2 in function f3tensor() ");
289 290 291 292
    te[mrl] += SHIFT_IND;
    te[mrl] -= mcl;
    
    te[mrl][mcl]=(float *) malloc((size_t)((mrow*mcol*mdep+SHIFT_IND)*sizeof(float)));
293
    if (!te[mrl][mcl]) declare_error("util.c: allocation failure 3 in function f3tensor() ");
294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
    te[mrl][mcl] += SHIFT_IND;
    te[mrl][mcl] -= mdl;
    
    for (e=mcl+1;e<=mch;e++) te[mrl][e]=te[mrl][e-1]+mdep;
    for (w=mrl+1;w<=mrh;w++){
        te[w]=te[w-1]+mcol;
        te[w][mcl]=te[w-1][mcl]+mcol*mdep;
        for (e=mcl+1;e<=mch;e++) te[w][e]=te[w][e-1]+mdep;
    }
    
    for (w=mrl;w<=mrh;w++)
        for (e=mcl;e<=mch;e++)
            for (r=mdl;r<=mdh;r++) te[w][e][r]=0.0;
    
    return te;
Tilman Steinweg's avatar
Tilman Steinweg committed
309 310 311 312
}


int ***i3tensor(int mrl, int mrh, int mcl, int mch,int mdl, int mdh){
313 314 315 316 317
    
    int w,e,r, mrow=mrh-mrl+1,mcol=mch-mcl+1,mdep=mdh-mdl+1;
    int ***te;
    
    te=(int ***) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(int**)));
318
    if (!te) declare_error("util.c: allocation failure 1 in function i3tensor() ");
319 320 321 322
    te += SHIFT_IND;
    te -= mrl;
    
    te[mrl]=(int **) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(int*)));
323
    if (!te[mrl]) declare_error("util.c: allocation failure 2 in function i3tensor() ");
324 325 326 327
    te[mrl] += SHIFT_IND;
    te[mrl] -= mcl;
    
    te[mrl][mcl]=(int *) malloc((size_t)((mrow*mcol*mdep+SHIFT_IND)*sizeof(int)));
328
    if (!te[mrl][mcl]) declare_error("util.c: allocation failure 3 in function i3tensor() ");
329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344
    te[mrl][mcl] += SHIFT_IND;
    te[mrl][mcl] -= mdl;
    
    for (e=mcl+1;e<=mch;e++) te[mrl][e]=te[mrl][e-1]+mdep;
    for (w=mrl+1;w<=mrh;w++){
        te[w]=te[w-1]+mcol;
        te[w][mcl]=te[w-1][mcl]+mcol*mdep;
        for (e=mcl+1;e<=mch;e++) te[w][e]=te[w][e-1]+mdep;
    }
    
    for (w=mrl;w<=mrh;w++)
        for (e=mcl;e<=mch;e++)
            for (r=mdl;r<=mdh;r++) te[w][e][r]=0.0;
    
    
    return te;
Tilman Steinweg's avatar
Tilman Steinweg committed
345 346 347 348
}


void free_vector(float *a, int ni, int nj){
349
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
350 351 352
}

void free_ivector(int *a, int ni, int nj){
353
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
354 355 356
}

void free_cvector(char *a, int ni, int nj){
357
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
358 359 360
}

void free_dvector(double *a, int ni, int nj){
361 362 363
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
}

Tilman Steinweg's avatar
Tilman Steinweg committed
364
void free_matrix(float **ma, int mrl, int mrh, int mcl, int mch){
365 366
    free((FREE_ARGUMENT) (ma[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (ma+mrl-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
367 368 369
}

void free_imatrix(int **ma, int mrl, int mrh, int mcl, int mch){
370 371
    free((FREE_ARGUMENT) (ma[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (ma+mrl-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
372 373 374
}

void free_usmatrix(unsigned short int **ma, int mrl, int mrh, int mcl, int mch){
375 376
    free((FREE_ARGUMENT) (ma[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (ma+mrl-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
377 378 379
}

void free_f3tensor(float ***te, int mrl, int mrh, int mcl, int mch, int mdl, int mdh){
380 381 382
    free((FREE_ARGUMENT) (te[mrl][mcl]+mdl-SHIFT_IND));
    free((FREE_ARGUMENT) (te[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (te+mrl-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
383 384 385
}

void free_i3tensor(int ***te, int mrl, int mrh, int mcl, int mch, int mdl, int mdh){
386 387 388
    free((FREE_ARGUMENT) (te[mrl][mcl]+mdl-SHIFT_IND));
    free((FREE_ARGUMENT) (te[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (te+mrl-SHIFT_IND));
Tilman Steinweg's avatar
Tilman Steinweg committed
389 390 391
}

void zero(float *A, int u_max){
392 393
    int u=0;
    while(++u<=u_max) *(A++)=0.0;
Tilman Steinweg's avatar
Tilman Steinweg committed
394 395 396
}

void normalize_data(float **data, int ntr, int ns){
397 398 399 400 401 402 403 404 405
    
    float *max_min_trace=NULL;
    float max, min, tmp=0.0;
    int i,j;
    
    
    max_min_trace = vector(1,ntr);
    for(i=1;i<=ntr;i++){
        
Tilman Steinweg's avatar
Tilman Steinweg committed
406 407
      		max=0.0;
      		min=0.0;
408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
        
        for(j=2;j<=ns;j++){
            /* Looking for max and min */
            if(data[i][j]>max){max = data[i][j];}
            if(data[i][j]<min){min = data[i][j];}
            
            if (max>(fabs(min))){tmp=max;}
            else{tmp=fabs(min);}
        }
        
        max_min_trace[i]=tmp;
        
        /* set maximum_values=0.0 to 1.0 */
        if(max_min_trace[i]==0.0){max_min_trace[i]=1.0;}
        
        
        for(j=2;j<=ns;j++){
            data[i][j] = data[i][j]/max_min_trace[i];
        }
        
    }
Tilman Steinweg's avatar
Tilman Steinweg committed
429 430 431
}

/*  quickSort
432 433
 * This public-domain C implementation by Darel Rex Finley .
 * http://alienryderflex.com/quicksort */
Tilman Steinweg's avatar
Tilman Steinweg committed
434 435

void quicksort(float *arr, int dummy, int elements) {
436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463
    
    /* number of levels: 64 is enough for 1.8e19 elements to be sorted */
#define  MAX_LEVELS  64
    
    int	beg[MAX_LEVELS], end[MAX_LEVELS], i=0, L, R, swap ;
    float	piv;
    
    beg[0]=0;
    end[0]=elements;
    
    while (i>=0) {
        L=beg[i]; R=end[i]-1;
        if (L<R) {
            piv=arr[L];
            while (L<R) {
                while (arr[R]>=piv && L<R) R--;
                if (L<R) arr[L++]=arr[R];
                while (arr[L]<=piv && L<R) L++;
                if (L<R) arr[R--]=arr[L];
            }
            arr[L]=piv; beg[i+1]=L+1; end[i+1]=end[i]; end[i++]=L;
            if (end[i]-beg[i]>end[i-1]-beg[i-1]) {
                swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap;
                swap=end[i]; end[i]=end[i-1]; end[i-1]=swap;
            }
        }
        else	i--;
    }
Tilman Steinweg's avatar
Tilman Steinweg committed
464
}