Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 12416)
+++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 12417)
@@ -19,9 +19,10 @@
 #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
 
Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp	(revision 12416)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp	(revision 12417)
@@ -13,15 +13,11 @@
 #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,5 +29,4 @@
 	/*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,9 +58,43 @@
 	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;
+}/*}}}*/