Changeset 12417

Timestamp:

06/15/12 11:40:32 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added simpler gsl slover

Location:

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

Files:

• Solverx.h (modified) (1 diff)
• SolverxGsl.cpp (modified) (3 diffs)

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

@@ -19 +19 @@
 #ifdef _HAVE_PETSC_
 void    SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);
-void    DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
+void  DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
 #endif
 
 #ifdef _HAVE_GSL_
-void    SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf);
+void SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf);
+void SolverxGsl(double** pX,double* A,double* B,int n);
 #endif
 

issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp

@@ -13 +13 @@
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
 #endif
- 
 #include <gsl/gsl_linalg.h>
 
-void    SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf){
+void SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf){/*{{{*/
 
-        /*intermediary: */
-        SeqMat* KffCopy=NULL;
-
-        /*Output: */
-        SeqVec*        uf               = NULL;
-
+        /*Intermediaries*/
+        SeqMat* KffCopy = NULL;
+        SeqVec*  uf = NULL;
         
         /*GSL Matrices and vectors: */
@@ -33 +29 @@
         /*We are going to do an in place LU decomp, so we need to save the matrix with its original structure: */
         KffCopy=Kff->Duplicate();
-
 
         /*Intermediary: */
@@ -63 +58 @@
         gsl_permutation_free (p);
         gsl_vector_free (x);
-
         delete KffCopy;
-
 
         /*Assign output pointers:*/
         *puf=uf;
-}
+}/*}}}*/
+void SolverxGsl(double** pX,double* A,double* B,int n){/*{{{*/
+
+        /*GSL Matrices and vectors: */
+        int              s;
+        gsl_matrix_view  a;
+        gsl_vector_view  b;
+        gsl_vector      *x = NULL;
+        gsl_permutation *p = NULL;
+
+        /*A will be modified by LU decomposition. Use copy*/
+        double* Acopy = (double*)xmalloc(n*n*sizeof(double));
+        memcpy(Acopy,A,n*n*sizeof(double));
+
+        /*Initialize gsl matrices and vectors: */
+        a = gsl_matrix_view_array (Acopy,n,n);
+        b = gsl_vector_view_array (B,n);
+        x = gsl_vector_alloc (n);
+
+        /*Run LU and solve: */
+        p = gsl_permutation_alloc (n);
+        gsl_linalg_LU_decomp (&a.matrix, p, &s);
+        gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
+
+        //printf ("x = \n");
+        //gsl_vector_fprintf (stdout, x, "%g");
+
+        /*Copy result*/
+        double* X = (double*)xmalloc(n*sizeof(double));
+        memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
+
+        /*Clean up and assign output pointer*/
+        xfree((void**)&Acopy);
+        gsl_permutation_free(p);
+        gsl_vector_free(x);
+        *pX=X;
+}/*}}}*/