Changeset 17328

Timestamp:

02/20/14 20:14:39 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added GslTriple solve to have a faster solution of Kriging

Location:

issm/trunk-jpl/src/c

Files:

• classes/kriging/Observations.cpp (modified) (1 diff)
• toolkits/gsl/DenseGslSolve.cpp (modified) (2 diffs)
• toolkits/gsl/gslincludes.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/kriging/Observations.cpp

@@ -401 +401 @@
         /*Solve the three linear systems*/
 #if _HAVE_GSL_
-        DenseGslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
-        DenseGslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
-        DenseGslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
+        DenseGslTripleSolve(&GinvG0,&Ginv1,&GinvZ,Gamma,Gamma,ones,obs,n_obs);
+        //DenseGslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
+        //DenseGslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
+        //DenseGslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
 #else
         _error_("GSL is required");

issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp

@@ -32 +32 @@
 }
 /*}}}*/
-void DenseGslSolve(/*output*/ IssmPDouble** px,/*stiffness matrix:*/ IssmPDouble* Kff, int Kff_M, int Kff_N, /*right hand side load vector: */ IssmPDouble* pf, int pf_M, Parameters* parameters){ /*{{{*/
+void DenseGslSolve(IssmPDouble** px,IssmPDouble* Kff,int Kff_M,int Kff_N,IssmPDouble* pf,int pf_M,Parameters* parameters){ /*{{{*/
 
         /*Intermediary: */
@@ -44 +44 @@
         /*allocate output pointers: */
         *px=x;
+}
+/*}}}*/
+void DenseGslTripleSolve(IssmPDouble** pX1,IssmPDouble** pX2,IssmPDouble** pX3,IssmPDouble* A,IssmPDouble* B1,IssmPDouble* B2,IssmPDouble* B3,int n){/*{{{*/
+
+        /*Allocate output: */
+        IssmPDouble *X1 = xNew<IssmPDouble>(n);
+        IssmPDouble *X2 = xNew<IssmPDouble>(n);
+        IssmPDouble *X3 = xNew<IssmPDouble>(n);
+
+#ifdef _HAVE_GSL_
+        /*GSL Matrices and vectors: */
+        int s;
+
+        double* Acopy = xNew<double>(n*n);
+        xMemCpy(Acopy,A,n*n);
+
+        /*Initialize gsl matrices and vectors: */
+        gsl_matrix_view a = gsl_matrix_view_array(Acopy,n,n);
+        gsl_vector_view b1= gsl_vector_view_array(B1,n);
+        gsl_vector_view b2= gsl_vector_view_array(B2,n);
+        gsl_vector_view b3= gsl_vector_view_array(B3,n);
+        gsl_vector_view x1= gsl_vector_view_array(X1,n);
+        gsl_vector_view x2= gsl_vector_view_array(X2,n);
+        gsl_vector_view x3= gsl_vector_view_array(X3,n);
+
+        /*Run LU and solve: */
+        gsl_permutation* p = gsl_permutation_alloc (n);
+        gsl_linalg_LU_decomp(&a.matrix,p,&s);
+
+        gsl_linalg_LU_solve(&a.matrix,p,&b1.vector,&x1.vector);
+        gsl_linalg_LU_solve(&a.matrix,p,&b2.vector,&x2.vector);
+        gsl_linalg_LU_solve(&a.matrix,p,&b3.vector,&x3.vector);
+
+        /*Clean up and assign output pointer*/
+        xDelete(Acopy);
+        gsl_permutation_free(p);
+#endif
+
+        /*allocate output pointers: */
+        *pX1=X1;
+        *pX2=X2;
+        *pX3=X3;
 }
 /*}}}*/

issm/trunk-jpl/src/c/toolkits/gsl/gslincludes.h

@@ -21 +21 @@
 
 void DenseGslSolve(IssmPDouble** pX,IssmPDouble* A,IssmPDouble* B, int n);
-void DenseGslSolve(IssmDouble** px, IssmDouble* Kff,int Kff_M, int Kff_N, IssmDouble* pf, int pf_M, Parameters* parameters);
+void DenseGslSolve(IssmDouble** pX,IssmDouble* Kff,int Kff_M,int Kff_N,IssmDouble* pf,int pf_M,Parameters* parameters);
+void DenseGslTripleSolve(IssmPDouble** pX1,IssmPDouble** pX2,IssmPDouble** pX3,IssmPDouble* A,IssmPDouble* B1,IssmPDouble* B2,IssmPDouble* B3,int n);
 
 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);