Added methods for local 2 global matrix

src/fd.h

@@ -526,6 +526,8 @@ float average_matrix(float ** matrix);
 float global_maximum(float ** gradiant_1);
 void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]);
 float matrix_product(float ** matrix1, float **matrix2);
+void get_local_from_global_matrix(float ** global_matrix,float ** local_matrix);
+float ** get_global_from_local_matrix(float ** local_matrix);
 
 /* L-BFGS */
 void lbfgs(float **grad1, float **grad2, float **grad3,float Vs_avg,float rho_avg,float Vp_avg, float *bfgsscale, float **bfgsmod, float **bfgsgrad,int bfgsnum,int bfgspar, int iteration, int * LBFGS_iter_start);

src/matrix_operations.c

@@ -21,7 +21,7 @@
 
 #include "fd.h"
 
-void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]){
+void write_matrix_disk(float ** local_matrix,char path_name[STRING_SIZE]){
     char joint[225];
     FILE *FPjoint;
     extern int POS[3],MYID;
@@ -32,7 +32,7 @@ void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]){
     
     for (i=1;i<=NX;i=i+IDX){
         for (j=1;j<=NY;j=j+IDY){
-            fwrite(&gradient[j][i],sizeof(float),1,FPjoint);
+            fwrite(&local_matrix[j][i],sizeof(float),1,FPjoint);
         }
     }
     
@@ -122,6 +122,72 @@ float matrix_product(float ** matrix1, float **matrix2) {
     return global_sum;
 }
 
+float ** get_global_from_local_matrix(float ** local_matrix) {
     
+    extern int NXG, NYG;
+    extern int NX,NY;
+    extern int POS[3];
+    
+    
+    float ** global_matrix=NULL,** global_matrix_temp=NULL;
+    int i=0,j=0;
+    int ii=0, jj=0;
+    
+    /* Allocate global matrix temp */
+    global_matrix_temp=matrix(1,NYG,1,NXG);
+    if(global_matrix_temp==NULL) {
+        declare_error("Allocation of global_matrix_temp in get_global_from_local_matrix failed!");
+    }
+    
+    /* Allocate global matrix */
+    global_matrix=matrix(1,NYG,1,NXG);
+    if(global_matrix==NULL) {
+        declare_error("Allocation of global_matrix in get_global_from_local_matrix failed!");
+    }
+    
+    /* Store local matrix in global matrix */
+    for (i=1;i<=NXG;i++){
+        for (j=1;j<=NYG;j++){
+            
+            if ( (POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY)) ) {
+                ii=i-POS[1]*NX;
+                jj=j-POS[2]*NY;
+                
+                global_matrix_temp[j][i]=local_matrix[jj][ii];
+            }
+        }
+    }
+    
+    MPI_Allreduce(&global_matrix_temp[1][1],&global_matrix[1][1],NXG*NYG,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
+    /* or (NXG+1)*(NYG+1) ? */
     
+    free_matrix(global_matrix_temp,1,NYG,1,NXG);
+    
+    return global_matrix;
+}
+
+void get_local_from_global_matrix(float ** global_matrix,float ** local_matrix) {
+    
+    extern int NXG, NYG;
+    extern int NX,NY;
+    extern int POS[3];
+    
+    int i=0,j=0;
+    int ii=0, jj=0;
+    
+    /* Store local matrix in global matrix */
+    for (i=1;i<=NXG;i++){
+        for (j=1;j<=NYG;j++){
+            
+            if ( (POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY)) ) {
+                ii=i-POS[1]*NX;
+                jj=j-POS[2]*NY;
+                
+                local_matrix[jj][ii]=global_matrix[j][i];
+            }
+        }
+    }
+    
+    free_matrix(global_matrix,1,NYG,1,NXG);
+}