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

Timestamp:

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

Author:

Mathieu Morlighem

Message:

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

File:

• issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp (modified) (2 diffs)

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;
 }
 /*}}}*/