Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp	(revision 11725)
+++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp	(revision 11726)
@@ -27,24 +27,28 @@
 
 	#ifdef _HAVE_PETSC_
-	Vec uf0_vector = NULL;
-	Vec df_vector  = NULL;
-	if(uf0) uf0_vector = uf0->vector;
-	if(df)  df_vector  = df->vector;
+		Vec uf0_vector = NULL;
+		Vec df_vector  = NULL;
+		if(uf0) uf0_vector = uf0->vector;
+		if(df)  df_vector  = df->vector;
 
-	/*In serial mode, the Petsc Options database has not been initialized properly: */
-	#ifdef _SERIAL_
-	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	PetscOptionsFromAnalysis(parameters,analysis_type);
-	#endif
+		/*In serial mode, the Petsc Options database has not been initialized properly: */
+		#ifdef _SERIAL_
+		parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+		PetscOptionsFromAnalysis(parameters,analysis_type);
+		#endif
 
-	SolverxPetsc(&uf->vector,Kff->matrix,pf->vector,uf0_vector,df_vector,parameters);
-	if(uf->vector == NULL){
-		uf->M = 0;
-	}
-	else{
-		VecGetSize(uf->vector,&uf->M);
-	}
+		SolverxPetsc(&uf->vector,Kff->matrix,pf->vector,uf0_vector,df_vector,parameters);
+		if(uf->vector == NULL){
+			uf->M = 0;
+		}
+		else{
+			VecGetSize(uf->vector,&uf->M);
+		}
 	#else
-	_error_("not supported yet!");
+		#ifdef _HAVE_GSL_
+		SolverxGsl(&uf->vector,Kff->matrix,pf->vector);
+		#else
+			_error_("not supported yet!");
+		#endif
 	#endif
 
Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 11725)
+++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 11726)
@@ -22,4 +22,8 @@
 #endif
 
+#ifdef _HAVE_GSL_
+void	SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf);
+#endif
+
 #endif  /* _SOLVERX_H */
 
Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp	(revision 11726)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp	(revision 11726)
@@ -0,0 +1,67 @@
+/*!\file SolverxGsl
+ * \brief Gsl implementation of solver
+ */
+
+#include "./Solverx.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../io/io.h"
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#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){
+
+	/*intermediary: */
+	SeqMat* KffCopy=NULL;
+
+	/*Output: */
+	SeqVec*        uf               = NULL;
+
+	/*We are going to do an in place LU decomp, so we need to save the matrix with its original structure: */
+	KffCopy=Kff->Duplicate();
+
+	/*GSL Matrices and vectors: */
+	gsl_matrix_view m;
+	gsl_vector_view b;
+	gsl_vector* x=NULL;
+	gsl_permutation* p=NULL;
+	
+	/*Intermediary: */
+	int M,N,N2,s;
+
+	Kff->GetSize(&M,&N);
+	pf->GetSize(&N2);
+
+	if(N!=N2)_error_("Right hand side vector of size %i, when matrix is of size %i-%i !",N2,M,N);
+	if(M!=N)_error_("Stiffness matrix should be square!");
+
+	/*Initialize gsl matrices and vectors: */
+	m = gsl_matrix_view_array (KffCopy->matrix, M, N);
+	b = gsl_vector_view_array (pf->vector, N);
+	x = gsl_vector_alloc (N);
+
+	/*Run LU and solve: */
+	p = gsl_permutation_alloc (N);
+	gsl_linalg_LU_decomp (&m.matrix, p, &s);
+	gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
+
+	//printf ("x = \n");
+	//gsl_vector_fprintf (stdout, x, "%g");
+
+	/*Get uf initialized with the results: */
+	uf=new SeqVec(gsl_vector_ptr(x,0),M);
+
+	/*Free resources:*/
+	gsl_permutation_free (p);
+	gsl_vector_free (x);
+	delete KffCopy;
+
+	/*Assign output pointers:*/
+	*puf=uf;
+}