matrix_operations.c 3.28 KB
Newer Older
1
/*-----------------------------------------------------------------------------------------
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
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
7 8 9
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
10
 * IFOS is distributed in the hope that it will be useful,
11 12 13 14 15
 * 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.
 *
 * 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 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
 -----------------------------------------------------------------------------------------*/
/*
 * This function provides several matrix operations
 */

#include "fd.h"

void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]){
    char joint[225];
    FILE *FPjoint;
    extern int POS[3],MYID;
    extern int IDX, IDY,NX,NY;
    int i,j;
    sprintf(joint,"%s.bin.%i.%i",path_name,POS[1],POS[2]);
    FPjoint=fopen(joint,"wb");
    
    for (i=1;i<=NX;i=i+IDX){
        for (j=1;j<=NY;j=j+IDY){
            fwrite(&gradient[j][i],sizeof(float),1,FPjoint);
        }
    }
    
    fclose(FPjoint);
    
    MPI_Barrier(MPI_COMM_WORLD);
    
    /* merge gradient file */
    sprintf(joint,"%s.bin",path_name);
    if (MYID==0) mergemod(joint,3);
    
    MPI_Barrier(MPI_COMM_WORLD);
    sprintf(joint,"%s.bin.%i.%i",path_name,POS[1],POS[2]);
    remove(joint);
};

/* Return global maximum of Matrix */
float global_maximum(float ** gradiant_1) {
    int i, j;
    float max=0.0,max1=0.0;
    extern int NX, NY;
57 58
    for (j=1;j<=NY;j++){
        for (i=1;i<=NX;i++){
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
            if((i*j == 1) || (max == 0.0)) {
                max = fabs(gradiant_1[j][i]);
            } else {
                if(fabs(gradiant_1[j][i]) > max) {
                    max = fabs(gradiant_1[j][i]);
                }
            }
        }
    }
    MPI_Allreduce(&max,  &max1,  1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
    return max1;
};

/* Return average  of Matrix */
float average_matrix(float ** matrix){
    
    extern int POS[3],MYID;
    extern int IDX, IDY,NX,NY;
    extern int NXG, NYG;
    
    int i,j;
    float local_sum=0;
    float global_sum;
    float buf1=0, buf2=0;
    float average;
    
85 86
    for (j=1;j<=NY;j=j+IDY){
        for (i=1;i<=NX;i=i+IDX){
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
            local_sum+=matrix[j][i];
        }
    }
    
    buf1=local_sum;
    buf2=0;
    MPI_Allreduce(&buf1,&buf2, 1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
    global_sum=buf2;
    
    average=global_sum/(NXG*NYG);
    
    MPI_Bcast(&average,1,MPI_FLOAT,0,MPI_COMM_WORLD);
    
    return average;
};

float matrix_product(float ** matrix1, float **matrix2) {
    
    extern int IDX, IDY,NX,NY;
    
    float local_sum=0.0;
    float global_sum=0.0;
    float buf1=0, buf2=0;
    int i,j;
    
112 113
    for (j=1;j<=NY;j=j+IDY){
        for (i=1;i<=NX;i=i+IDX){
114 115 116 117 118 119 120 121 122 123 124 125 126 127
            local_sum+=(matrix1[j][i]*matrix2[j][i]);
        }
    }
    
    buf1=local_sum; buf2=0;
    MPI_Allreduce(&buf1,&buf2, 1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
    global_sum=buf2;
    
    return global_sum;
}