util.c 12.9 KB
Newer Older
tilman.metz's avatar
tilman.metz committed
1
/*-----------------------------------------------------------------------------------------
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
tilman.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz 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.metz's avatar
tilman.metz committed
90 91 92
}

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

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

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


unsigned char *cvector(int ni, int nj){
126 127 128 129
    
    unsigned char *a;
    
    a=(unsigned char *)malloc((size_t) ((nj-ni+1+SHIFT_IND)*sizeof(unsigned char)));
130
    if (!a) declare_error("util.c: allocation failure in function cvector()");
131
    return a-ni+SHIFT_IND;
tilman.metz's avatar
tilman.metz committed
132 133 134 135
}


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

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

float **fmatrix(int mrl, int mrh, int mcl, int mch){
158 159 160 161 162
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    float **ma;
    
    ma=(float **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(float*)));
163
    if (!ma) declare_error("util.c: allocation failure 1 in function fmatrix() ");
164 165 166 167
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(float *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(float)));
168
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function fmatrix() ");
169 170 171 172 173 174 175 176 177
    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.metz's avatar
tilman.metz committed
178 179 180
}

float **matrix(int mrl, int mrh, int mcl, int mch){
181 182 183 184 185
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    float **ma;
    
    ma=(float **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(float*)));
186
    if (!ma) declare_error("util.c: allocation failure 1 in function matrix() ");
187 188 189 190
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(float *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(float)));
191
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function matrix() ");
192 193 194 195 196 197 198 199 200
    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.metz's avatar
tilman.metz committed
201 202 203 204
}


double **dmatrix(int mrl, int mrh, int mcl, int mch){
205 206 207 208 209
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    double **ma;
    
    ma=(double **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(double*)));
210
    if (!ma) declare_error("util.c: allocation failure 1 in function matrix() ");
211 212 213 214
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(double *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(double)));
215
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function dmatrix() ");
216 217 218 219 220 221 222 223 224
    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.metz's avatar
tilman.metz committed
225 226 227
}

int **imatrix(int mrl, int mrh, int mcl, int mch){
228 229 230 231 232
    
    int k,l, mrow=mrh-mrl+1,mcol=mch-mcl+1;
    int **ma;
    
    ma=(int **) malloc((size_t) ((mrow+SHIFT_IND)*sizeof(int*)));
233
    if (!ma) declare_error("util.c: allocation failure 1 in function imatrix() ");
234 235 236 237
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(int *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(int)));
238
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function imatrix() ");
239 240 241 242 243 244 245 246 247
    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.metz's avatar
tilman.metz committed
248 249 250
}

unsigned short int **usmatrix(int mrl, int mrh, int mcl, int mch){
251 252 253 254 255
    
    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*)));
256
    if (!ma) declare_error("util.c: allocation failure 1 in function usmatrix() ");
257 258 259 260
    ma += SHIFT_IND;
    ma -= mrl;
    
    ma[mrl]=(unsigned short int *) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(unsigned short int)));
261
    if (!ma[mrl]) declare_error("util.c: allocation failure 2 in function usmatrix() ");
262 263 264 265 266 267 268 269 270
    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.metz's avatar
tilman.metz committed
271 272 273
}

float ***f3tensor(int mrl, int mrh, int mcl, int mch,int mdl, int mdh){
274 275 276 277 278
    
    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**)));
279
    if (!te) declare_error("util.c: allocation failure 1 in function f3tensor() ");
280 281 282 283
    te += SHIFT_IND;
    te -= mrl;
    
    te[mrl]=(float **) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(float*)));
284
    if (!te[mrl]) declare_error("util.c: allocation failure 2 in function f3tensor() ");
285 286 287 288
    te[mrl] += SHIFT_IND;
    te[mrl] -= mcl;
    
    te[mrl][mcl]=(float *) malloc((size_t)((mrow*mcol*mdep+SHIFT_IND)*sizeof(float)));
289
    if (!te[mrl][mcl]) declare_error("util.c: allocation failure 3 in function f3tensor() ");
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304
    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.metz's avatar
tilman.metz committed
305 306 307 308
}


int ***i3tensor(int mrl, int mrh, int mcl, int mch,int mdl, int mdh){
309 310 311 312 313
    
    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**)));
314
    if (!te) declare_error("util.c: allocation failure 1 in function i3tensor() ");
315 316 317 318
    te += SHIFT_IND;
    te -= mrl;
    
    te[mrl]=(int **) malloc((size_t)((mrow*mcol+SHIFT_IND)*sizeof(int*)));
319
    if (!te[mrl]) declare_error("util.c: allocation failure 2 in function i3tensor() ");
320 321 322 323
    te[mrl] += SHIFT_IND;
    te[mrl] -= mcl;
    
    te[mrl][mcl]=(int *) malloc((size_t)((mrow*mcol*mdep+SHIFT_IND)*sizeof(int)));
324
    if (!te[mrl][mcl]) declare_error("util.c: allocation failure 3 in function i3tensor() ");
325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
    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.metz's avatar
tilman.metz committed
341 342 343 344
}


void free_vector(float *a, int ni, int nj){
345
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
tilman.metz's avatar
tilman.metz committed
346 347 348
}

void free_ivector(int *a, int ni, int nj){
349
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
tilman.metz's avatar
tilman.metz committed
350 351 352
}

void free_cvector(char *a, int ni, int nj){
353
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
tilman.metz's avatar
tilman.metz committed
354 355 356
}

void free_dvector(double *a, int ni, int nj){
357 358 359
    free((FREE_ARGUMENT) (a+ni-SHIFT_IND));
}

tilman.metz's avatar
tilman.metz committed
360
void free_matrix(float **ma, int mrl, int mrh, int mcl, int mch){
361 362
    free((FREE_ARGUMENT) (ma[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (ma+mrl-SHIFT_IND));
tilman.metz's avatar
tilman.metz committed
363 364 365
}

void free_imatrix(int **ma, int mrl, int mrh, int mcl, int mch){
366 367
    free((FREE_ARGUMENT) (ma[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (ma+mrl-SHIFT_IND));
tilman.metz's avatar
tilman.metz committed
368 369 370
}

void free_usmatrix(unsigned short int **ma, int mrl, int mrh, int mcl, int mch){
371 372
    free((FREE_ARGUMENT) (ma[mrl]+mcl-SHIFT_IND));
    free((FREE_ARGUMENT) (ma+mrl-SHIFT_IND));
tilman.metz's avatar
tilman.metz committed
373 374 375
}

void free_f3tensor(float ***te, int mrl, int mrh, int mcl, int mch, int mdl, int mdh){
376 377 378
    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.metz's avatar
tilman.metz committed
379 380 381
}

void free_i3tensor(int ***te, int mrl, int mrh, int mcl, int mch, int mdl, int mdh){
382 383 384
    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.metz's avatar
tilman.metz committed
385 386 387
}

void zero(float *A, int u_max){
388 389
    int u=0;
    while(++u<=u_max) *(A++)=0.0;
tilman.metz's avatar
tilman.metz committed
390 391 392
}

void normalize_data(float **data, int ntr, int ns){
393 394 395 396 397 398 399 400 401
    
    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.metz's avatar
tilman.metz committed
402 403
      		max=0.0;
      		min=0.0;
404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424
        
        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.metz's avatar
tilman.metz committed
425 426 427
}

/*  quickSort
428 429
 * This public-domain C implementation by Darel Rex Finley .
 * http://alienryderflex.com/quicksort */
tilman.metz's avatar
tilman.metz committed
430 431

void quicksort(float *arr, int dummy, int elements) {
432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459
    
    /* 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.metz's avatar
tilman.metz committed
460
}