Changeset 12174

Timestamp:

05/02/12 12:56:30 (13 years ago)

Author:

Mathieu Morlighem

Message:

Final version of Kriging (still not working)

Location:

issm/trunk-jpl/src/c/modules/Krigingx

Files:

• Krigingx.cpp (modified) (1 diff)
• Krigingx.h (modified) (1 diff)

issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp

@@ -15 +15 @@
         double *predictions = NULL;
 
-        /*Allocate output*/
-        predictions=(double*)xmalloc(n_interp*sizeof(double));
+        /*Intermediaries*/
+        int     i,j;
+        double  numerator,denominator,ratio;
+        double *Gamma     = NULL;
+        double *GinvG0    = NULL;
+        double *Ginv1     = NULL;
+        double *GinvZ     = NULL;
+        double *gamma0    = NULL;
+        double *ones      = NULL;
 
-        /*Assign output pointer*/
+        /*Memory allocation*/
+        predictions =(double*)xmalloc(n_interp*sizeof(double));
+        Gamma       =(double*)xmalloc(n_obs*n_obs*sizeof(double));
+        gamma0      =(double*)xmalloc(n_obs*sizeof(double));
+        ones        =(double*)xmalloc(n_obs*sizeof(double));
+
+        /*First: Create semivariogram matrix for observations*/
+        for(i=0;i<n_obs;i++){
+                for(j=0;j<=i;j++){
+                        SemiVariogram(&Gamma[i*n_obs+j],x[i],y[i],x[j],y[j]);
+                        Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
+                }
+        }
+        for(i=0;i<n_obs;i++) ones[i]=1;
+
+        /*Loop over all interpolations*/
+        for(int idx=0;idx<n_interp;idx++){
+
+                /*Get semivariogram vector associated to this location*/
+                for(i=0;i<n_obs;i++) SemiVariogram(&gamma0[i],x[i],y[i],x_interp[idx],y_interp[idx]);
+
+                /*Solve the three linear systems*/
+                MumpsSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
+                MumpsSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
+                MumpsSolve(&GinvZ, Gamma,observations,n_obs); // Gamma^-1 Z
+
+                /*Prepare predictor*/
+                numerator=0.; denominator=0.;
+                for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
+                for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
+                ratio=numerator/denominator;
+
+                predictions[idx]=0;
+                for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
+
+                /*clean-up*/
+                xfree((void**)&GinvG0);
+                xfree((void**)&Ginv1);
+                xfree((void**)&GinvZ);
+        }
+
+        /*clean-up and Assign output pointer*/
+        xfree((void**)&Gamma);
+        xfree((void**)&gamma0);
+        xfree((void**)&ones);
         *ppredictions=predictions;
 }
+
+void SemiVariogram(double* gamma,double x1,double y1,double x2,double y2){
+
+        /*Calculate distance*/
+        double r=sqrt(pow(x1-x2,2.)+pow(y1-y2,2.));
+
+        /*Switch between variogram models*/
+        switch(1){
+                case 1:{ /*Exponential*/
+                        double c0=0.2;
+                        double c1=0.8;
+                        double a =1;
+                        *gamma = c0 + c1*(1-exp(-r/a));
+                        return;}
+                default:
+                        _error_("Not implemented yet");
+        }
+}
+
+void MumpsSolve(double** Xdouble,double* Adouble,double* Bdouble,int n){
+#ifdef _HAVE_PETSC_
+
+        /*Intermediaries*/
+        Mat A = NULL;
+        Vec B = NULL;
+        Vec X = NULL;
+        KSP ksp = NULL;
+        PC  pc  = NULL;
+
+        /*Create parameters and select MUMPS solver*/
+        PetscOptionsSetFromOptions();
+        PetscOptionsClear();
+        PetscOptionsInsertMultipleString("-mat_type mpiaij \
+                                -ksp_type preonly \
+                                -pc_type lu \
+                                -pc_factor_mat_solver_package mumps \
+                                -mat_mumps_icntl_14 120 \
+                                -pc_factor_shift_positive_definite true");
+
+        /*Create indexing and ones*/
+        int* idx=(int*)xmalloc(n*sizeof(int));
+        for(int i=0;i<n;i++) idx[i]=i;
+
+        /*Convert matrix to Dense PETSc matrix*/
+        MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
+        MatSetValues(A,n,idx,n,idx,Adouble,INSERT_VALUES);
+        MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); 
+        MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
+
+        /*Convert vector to PETSc vector*/
+        VecCreateSeq(PETSC_COMM_SELF,n,&B);
+        VecSetValues(B,n,idx,Bdouble,INSERT_VALUES);
+        VecAssemblyBegin(B);
+        VecAssemblyEnd(B);
+
+        /*Allocate output*/
+        VecCreateSeq(PETSC_COMM_SELF,n,&X);
+
+        /*Go Solve*/
+        KSPCreate(MPI_COMM_WORLD,&ksp);
+        KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
+        KSPSetFromOptions(ksp);
+        KSPGetPC(ksp,&pc);
+        PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
+        KSPSolve(ksp,B,X);
+        KSPFree(&ksp);
+
+        /*Serialize vector*/
+        VecToMPISerial(Xdouble,X);
+#else
+        error_("PETSc support required");
+#endif
+}

issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h

@@ -9 +9 @@
 #include "../../toolkits/toolkits.h"
 
-
 int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp);
+void SemiVariogram(double* gamma,double x1,double y1,double x2,double y2);
+void MumpsSolve(double** X,double* A,double* B,int n); //Solve X for the linear system AX = B
 
 #endif /* _KRIGINGX_H */
-