SuperLU update in OpenNL + debug in renderTopoSurface plugin

CMakeLists.txt

@@ -33,7 +33,7 @@ find_package(Boost COMPONENTS regex thread system REQUIRED)
 find_package(ZLIB REQUIRED)
 find_package(LibXml2 REQUIRED)
 find_package(GLEW REQUIRED)
-find_package(SuiteSparse REQUIRED)
+#find_package(SuiteSparse REQUIRED)
 
 IF (DEFINED ASSERTON)
 	add_definitions(-DCGOGN_ASSERT_BOOL=${ASSERTON})
@@ -73,7 +73,7 @@ SET (CGoGN_EXT_INCLUDES
 	${ZLIB_INCLUDE_DIRS}
 	${LIBXML2_INCLUDE_DIR}
 	${Boost_INCLUDE_DIRS}
-
+#	${SUITESPARSE_INCLUDE_DIRS}
 )
 
 # define libs for external libs
@@ -86,7 +86,7 @@ SET (CGoGN_EXT_LIBS
 	${Boost_SYSTEM_LIBRARY}
 	${Boost_REGEX_LIBRARY}
 	${Boost_THREAD_LIBRARY}
-	${SUITESPARSE_LIBRARIES}
+#	${SUITESPARSE_LIBRARIES}
 )
 
 #optionnal libs

SCHNApps/CMakeLists.txt

@@ -18,6 +18,8 @@ find_package(Qt4 REQUIRED)
 find_package(QGLViewer REQUIRED)
 find_package(PythonLibs 2.7 REQUIRED)
 find_package(SuiteSparse REQUIRED)
+find_package(SuperLU REQUIRED)
+
 
 SET( QT_USE_QTOPENGL TRUE )
 SET( QT_USE_QTXML TRUE )
@@ -68,6 +70,7 @@ SET (EXT_INCLUDES
 	${QGLVIEWER_INCLUDE_DIR}
 	${PYTHON_INCLUDE_DIRS}
 	${SUITESPARSE_INCLUDE_DIRS}
+	${SUPERLU_INCLUDE_DIRS}
 )
 
 # define libs for external libs
@@ -85,6 +88,7 @@ SET (EXT_LIBS
 	${QGLVIEWER_LIBRARIES}
 	${PYTHON_LIBRARIES}
 	${SUITESPARSE_LIBRARIES}
+	${SUPERLU_LIBRARIES}
 )
 
 

SCHNApps/Plugins/differentialProperties/src/computeNormalDialog.cpp

@@ -28,7 +28,6 @@ ComputeNormalDialog::ComputeNormalDialog(Window* w) :
 	{
 		QListWidgetItem* item = new QListWidgetItem(map->getName(), mapList);
 		item->setCheckState(Qt::Unchecked);
-		connect(map, SIGNAL(attributeModified(unsigned int, QString)), this, SLOT(attributeModified(unsigned int, QString)));
 	}
 }
 

SCHNApps/Plugins/renderTopoSurface/src/renderTopoSurface.cpp

@@ -95,6 +95,8 @@ void RenderTopoSurfacePlugin::viewLinked(View* view, Plugin* plugin)
 		const QList<MapHandlerGen*>& maps = view->getLinkedMaps();
 		foreach(MapHandlerGen* mh, maps)
 		{
+			connect(mh, SIGNAL(attributeModified(unsigned int, QString)), this, SLOT(attributeModified(unsigned int, QString)));
+			connect(mh, SIGNAL(connectivityModified()), this, SLOT(connectivityModified()));
 			PerMapParameterSet* p = new PerMapParameterSet(mh);
 			registerShader(p->m_renderTopo->shader1());
 			registerShader(p->m_renderTopo->shader2());
@@ -119,6 +121,8 @@ void RenderTopoSurfacePlugin::viewUnlinked(View* view, Plugin* plugin)
 		QHash<QString, PerMapParameterSet*>::const_iterator i = params->perMap.constBegin();
 		while (i != params->perMap.constEnd())
 		{
+			disconnect(i.value()->mh, SIGNAL(attributeModified(unsigned int, QString)), this, SLOT(attributeModified(unsigned int, QString)));
+			disconnect(i.value()->mh, SIGNAL(connectivityModified()), this, SLOT(connectivityModified()));
 			unregisterShader(i.value()->m_renderTopo->shader1());
 			unregisterShader(i.value()->m_renderTopo->shader2());
 			delete i.value();
@@ -146,10 +150,10 @@ void RenderTopoSurfacePlugin::mapLinked(MapHandlerGen* m)
 	connect(m, SIGNAL(connectivityModified()), this, SLOT(connectivityModified()));
 
 	ParameterSet* params = h_viewParams[view];
-	PerMapParameterSet* p = new PerMapParameterSet(m);
-	registerShader(p->m_renderTopo->shader1());
-	registerShader(p->m_renderTopo->shader2());
-	params->perMap.insert(m->getName(), p);
+	PerMapParameterSet* perMap = new PerMapParameterSet(m);
+	registerShader(perMap->m_renderTopo->shader1());
+	registerShader(perMap->m_renderTopo->shader2());
+	params->perMap.insert(m->getName(), perMap);
 	if(params->selectedMap == NULL || params->perMap.count() == 1)
 		changeSelectedMap(view, m);
 	else
@@ -165,7 +169,10 @@ void RenderTopoSurfacePlugin::mapUnlinked(MapHandlerGen* m)
 	disconnect(m, SIGNAL(connectivityModified()), this, SLOT(connectivityModified()));
 
 	ParameterSet* params = h_viewParams[view];
-	delete params->perMap[m->getName()];
+	PerMapParameterSet* perMap = params->perMap[m->getName()];
+	unregisterShader(perMap->m_renderTopo->shader1());
+	unregisterShader(perMap->m_renderTopo->shader2());
+	delete perMap;
 	params->perMap.remove(m->getName());
 
 	if(params->selectedMap == m)

SCHNApps/bin/init.py

 importPlugin = schnapps.loadPlugin("ImportSurface");
 renderPlugin = schnapps.loadPlugin("Render");
-renderVectorPlugin = schnapps.loadPlugin("RenderVector");
-renderTopoSurfacePlugin = schnapps.loadPlugin("RenderTopoSurface");
-differentialPropertiesPlugin = schnapps.loadPlugin("DifferentialProperties");
-subdivisionPlugin = schnapps.loadPlugin("SubdivideSurface");
+#renderVectorPlugin = schnapps.loadPlugin("RenderVector");
+#renderTopoSurfacePlugin = schnapps.loadPlugin("RenderTopoSurface");
+#differentialPropertiesPlugin = schnapps.loadPlugin("DifferentialProperties");
+#subdivisionPlugin = schnapps.loadPlugin("SubdivideSurface");
 surfaceDeformationPlugin = schnapps.loadPlugin("SurfaceDeformation");
 
 obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/lowRes/iphi_good_9k.off");
@@ -11,11 +11,11 @@ obj = importPlugin.importFromFile("/home/kraemer/Media/Data/surface/lowRes/iphi_
 v = schnapps.getView("view_0");
 
 schnapps.linkViewAndPlugin(v.getName(), renderPlugin.getName());
-schnapps.linkViewAndPlugin(v.getName(), renderVectorPlugin.getName());
-schnapps.linkViewAndPlugin(v.getName(), renderTopoSurfacePlugin.getName());
+#schnapps.linkViewAndPlugin(v.getName(), renderVectorPlugin.getName());
+#schnapps.linkViewAndPlugin(v.getName(), renderTopoSurfacePlugin.getName());
 schnapps.linkViewAndPlugin(v.getName(), surfaceDeformationPlugin.getName());
 
 schnapps.linkViewAndMap(v.getName(), obj.getName());
 
-differentialPropertiesPlugin.computeNormal(obj.getName());
-differentialPropertiesPlugin.computeCurvature(obj.getName());
+#differentialPropertiesPlugin.computeNormal(obj.getName());
+#differentialPropertiesPlugin.computeCurvature(obj.getName());

SCHNApps/src/colorComboBox.cpp

@@ -23,18 +23,15 @@ void ColorComboBox::populateList()
 	//QStringList colorNames = QColor::colorNames();
 	QStringList colorNames;
 	colorNames <<
-		"darkGreen" <<
-		"green" <<
-		"gray" <<
 		"red" <<
-		"white" <<
+		"green" <<
 		"blue" <<
 		"cyan "<<
-		"darkMagenta" <<
+		"magenta" <<
 		"yellow" <<
-		"darkRed" <<
-		"black" <<
-		"magenta";
+		"gray" <<
+		"white" <<
+		"black";
 
 	for (int i = 0; i < colorNames.size(); ++i) {
 		QColor color(colorNames[i]);

ThirdParty/OpenNL/CMakeOptions.txt

@@ -55,7 +55,8 @@ SET(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/OpenNL/CMakeModules)
 # (with a single underscore) rather than -DAdd__ (with two
 # underscores). If you got two underscores, remove the 
 # second one, make clean and re-make SuperLU.
-#SET(USE_SUPERLU TRUE)
+SET(USE_SUPERLU TRUE)
+find_package(SuperLU REQUIRED)
 
 # Edit these lines to match your installation
 #SET(SUPERLU_HOME /home/yourself/SuperLU_4.0 )

ThirdParty/OpenNL/src/CMakeLists.txt

@@ -20,8 +20,8 @@ ENDIF(USE_ATLAS)
 
 IF(USE_SUPERLU)
     ADD_DEFINITIONS(-DNL_USE_SUPERLU)
-    SET(SUPERLU_INCS ${SUPERLU_HOME}/SRC)
-    INCLUDE_DIRECTORIES(${SUPERLU_INCS})
+#   SET(SUPERLU_INCS ${SUPERLU_HOME}/SRC)
+	INCLUDE_DIRECTORIES(${SUPERLU_INCLUDE_DIRS})
 ENDIF(USE_SUPERLU)
 
 IF(USE_CHOLMOD)

ThirdParty/OpenNL/src/NL/nl.h

@@ -121,9 +121,6 @@ typedef void* NLContext ;
 #define NL_ELAPSED_TIME     0x10a
 #define NL_PRECONDITIONER   0x10b
 
-//#define NL_DIRECT           0x10c
-//#define NL_FACTORIZED       0x10d
-
 /* Solvers */
 
 #define NL_CG                    0x200
@@ -133,7 +130,6 @@ typedef void* NLContext ;
 #define NL_PERM_SUPERLU_EXT      0x211
 #define NL_SYMMETRIC_SUPERLU_EXT 0x212
 #define NL_SOLVER_USER           0x213
-
 #define NL_CHOLMOD_EXT           0x214
 
 #define NL_CNC_FLOAT_CRS         0x220

ThirdParty/OpenNL/src/NL/nl_api.c

@@ -517,8 +517,11 @@ void nlEndMatrix() {
         ) ;
     }
 
-	if(nlCurrentContext->solver == NL_CHOLMOD_EXT && // or any other direct solver
-	   nlCurrentContext->direct_solver_context == NULL) {
+	if((nlCurrentContext->solver == NL_CHOLMOD_EXT || // or any other direct solver
+		nlCurrentContext->solver == NL_SUPERLU_EXT ||
+		nlCurrentContext->solver == NL_PERM_SUPERLU_EXT ||
+		nlCurrentContext->solver == NL_SYMMETRIC_SUPERLU_EXT) &&
+			nlCurrentContext->direct_solver_context == NULL) {
 		nlCurrentContext->factorize_func() ;
 	}
 }
@@ -567,9 +570,9 @@ void nlEndRow() {
     NLRowColumn*    al = &nlCurrentContext->al ;
     NLRowColumn*    xl = &nlCurrentContext->xl ;
     NLSparseMatrix* M  = &nlCurrentContext->M  ;
-    NLdouble* b        = nlCurrentContext->b ;
-    NLuint nf          = af->size ;
-    NLuint nl          = al->size ;
+	NLdouble*       b  = nlCurrentContext->b ;
+	NLuint          nf = af->size ;
+	NLuint          nl = al->size ;
     NLuint current_row = nlCurrentContext->current_row ;
     NLuint i ;
     NLuint j ;
@@ -674,15 +677,15 @@ void nlEnd(NLenum prim) {
 /* nlSolve() driver routine */
 
 NLboolean nlSolve() {
-    NLboolean result ;
-    NLdouble start_time = nlCurrentTime() ; 
-    nlCheckState(NL_STATE_SYSTEM_CONSTRUCTED) ;
-    nlCurrentContext->elapsed_time = 0 ;
+	NLboolean result ;
+	NLdouble start_time = nlCurrentTime() ;
+	nlCheckState(NL_STATE_SYSTEM_CONSTRUCTED) ;
+	nlCurrentContext->elapsed_time = 0 ;
 	result = nlCurrentContext->solver_func() ;
-    nlVectorToVariables() ;
+	nlVectorToVariables() ;
 	nlCurrentContext->elapsed_time = nlCurrentTime() - start_time ;
-    nlTransition(NL_STATE_SYSTEM_CONSTRUCTED, NL_STATE_SOLVED) ;
-    return result ;
+	nlTransition(NL_STATE_SYSTEM_CONSTRUCTED, NL_STATE_SOLVED) ;
+	return result ;
 }
 
 /************************************************************************/
@@ -717,7 +720,10 @@ void nlReset(NLboolean keep_matrix) {
 				nlUnlockVariable(i);
 			}
 			nlCurrentContext->matrix_already_set = NL_FALSE ;
-			if(nlCurrentContext->solver == NL_CHOLMOD_EXT) { // or any other direct solver
+			if( nlCurrentContext->solver == NL_CHOLMOD_EXT || // or any other direct solver
+				nlCurrentContext->solver == NL_SUPERLU_EXT ||
+				nlCurrentContext->solver == NL_PERM_SUPERLU_EXT ||
+				nlCurrentContext->solver == NL_SYMMETRIC_SUPERLU_EXT) {
 				nlCurrentContext->clear_factor_func() ;
 			}
 		}

ThirdParty/OpenNL/src/NL/nl_cholmod.c

@@ -57,8 +57,8 @@
 
 typedef struct {
 	cholmod_common c ;
-	cholmod_factor* cL ;
-	cholmod_dense* cb ;
+	cholmod_factor* cL ; // factor
+	cholmod_dense* cb ;  // right-hand side
 } cholmod_context ;
 
 
@@ -66,9 +66,9 @@ typedef struct {
 NLboolean nlFactorize_CHOLMOD() {
 
 	/* OpenNL Context */
-	NLSparseMatrix* M = &(nlCurrentContext->M) ;
-	NLuint n   = nlCurrentContext->n ;
-	NLuint nnz = nlSparseMatrixNNZ(M) ; /* Number of Non-Zero coeffs */
+	NLSparseMatrix* M   = &(nlCurrentContext->M) ;
+	NLuint          n   = nlCurrentContext->n ;
+	NLuint          nnz = nlSparseMatrixNNZ(M) ; /* Number of Non-Zero coeffs */
 
 	cholmod_context* context = (cholmod_context*)(nlCurrentContext->direct_solver_context) ;
 	if(context == NULL) {
@@ -87,13 +87,13 @@ NLboolean nlFactorize_CHOLMOD() {
 	nl_assert(M->storage & NL_MATRIX_STORE_COLUMNS) ;
 	nl_assert(M->m == M->n) ;
 
+	cholmod_start(&(context->c)) ;
+
 	/*
 	 * Step 1: convert matrix M into CHOLMOD compressed column representation
 	 * ----------------------------------------------------------------------
 	 */
 
-	cholmod_start(&(context->c)) ;
-
 	cA = cholmod_allocate_sparse(n, n, nnz, NL_FALSE, NL_TRUE, -1, CHOLMOD_REAL, &(context->c)) ;
 	int* colptr = (int*)(cA->p) ;
 	int* rowind = (int*)(cA->i) ;
@@ -137,7 +137,7 @@ NLboolean nlSolve_CHOLMOD() {
 	/* OpenNL Context */
 	NLdouble* b = nlCurrentContext->b ;
 	NLdouble* x = nlCurrentContext->x ;
-	NLuint n    = nlCurrentContext->n ;
+	NLuint    n = nlCurrentContext->n ;
 
 	cholmod_context* context = (cholmod_context*)(nlCurrentContext->direct_solver_context) ;
 	nl_assert(context != NULL) ;
@@ -154,9 +154,8 @@ NLboolean nlSolve_CHOLMOD() {
 	 */
 
 	double* cbx = (double*)(context->cb->x) ;
-	for(i = 0; i < n; i++) {
+	for(i = 0; i < n; i++)
 		cbx[i] = b[i] ;
-	}
 
 	/*
 	 * Step 2: solve
@@ -207,4 +206,9 @@ NLboolean nlSolve_CHOLMOD() {
 	return NL_FALSE ;
 }
 
+NLboolean nlClear_CHOLMOD() {
+	nl_assert_not_reached ;
+	return NL_FALSE ;
+}
+
 #endif

ThirdParty/OpenNL/src/NL/nl_context.c

@@ -108,7 +108,10 @@ void nlDeleteContext(NLContext context_in) {
     if(context->alloc_b) {
         NL_DELETE_ARRAY(context->b) ;
     }
-	if(context->solver == NL_CHOLMOD_EXT) { // or any other direct solver
+	if( context->solver == NL_CHOLMOD_EXT || // or any other direct solver
+		context->solver == NL_SUPERLU_EXT ||
+		context->solver == NL_PERM_SUPERLU_EXT ||
+		context->solver == NL_SYMMETRIC_SUPERLU_EXT) {
 		context->clear_factor_func() ;
 	}
 
@@ -268,10 +271,12 @@ NLboolean nlDefaultFactorize() {
 	case NL_CNC_FLOAT_ELL:
 	case NL_CNC_DOUBLE_ELL:
 	case NL_CNC_FLOAT_HYB:
-	case NL_CNC_DOUBLE_HYB:
+	case NL_CNC_DOUBLE_HYB: break ;
 	case NL_SUPERLU_EXT:
 	case NL_PERM_SUPERLU_EXT:
-	case NL_SYMMETRIC_SUPERLU_EXT: break ;
+	case NL_SYMMETRIC_SUPERLU_EXT: {
+			result = nlFactorize_SUPERLU() ;
+	} break ;
 	case NL_CHOLMOD_EXT: {
 		result = nlFactorize_CHOLMOD() ;
 	} break ;
@@ -293,10 +298,12 @@ void nlDefaultClearFactor() {
 	case NL_CNC_FLOAT_ELL:
 	case NL_CNC_DOUBLE_ELL:
 	case NL_CNC_FLOAT_HYB:
-	case NL_CNC_DOUBLE_HYB:
+	case NL_CNC_DOUBLE_HYB: break ;
 	case NL_SUPERLU_EXT:
 	case NL_PERM_SUPERLU_EXT:
-	case NL_SYMMETRIC_SUPERLU_EXT: break ;
+	case NL_SYMMETRIC_SUPERLU_EXT: {
+		nlClear_SUPERLU() ;
+	} break ;
 	case NL_CHOLMOD_EXT: {
 		nlClear_CHOLMOD() ;
 	} break ;

ThirdParty/OpenNL/src/NL/nl_superlu.c

@@ -54,175 +54,232 @@
 #include <slu_cdefs.h>
 #include <supermatrix.h>
 
-/* Note: SuperLU is difficult to call, but it is worth it.    */
-/* Here is a driver inspired by A. Sheffer's "cow flattener". */
+
+typedef struct {
+	superlu_options_t options ;
+	SuperLUStat_t stat ;
+	SuperMatrix L ;
+	SuperMatrix U ;
+	NLint* perm_c ;
+	NLint* perm_r;
+} superlu_context ;
+
+
+
+NLboolean nlFactorize_SUPERLU() {
+
+	/* OpenNL Context */
+	NLSparseMatrix* M   = &(nlCurrentContext->M) ;
+	NLuint          n   = nlCurrentContext->n ;
+	NLuint          nnz = nlSparseMatrixNNZ(M) ; /* Number of Non-Zero coeffs */
+
+	superlu_context* context = (superlu_context*)(nlCurrentContext->direct_solver_context) ;
+	if(context == NULL) {
+		nlCurrentContext->direct_solver_context = malloc(sizeof(superlu_context)) ;
+		context = (superlu_context*)(nlCurrentContext->direct_solver_context) ;
+	}
+
+	/* SUPERLU variables */
+	NLint info ;
+	SuperMatrix A, AC ;
+
+	/* Temporary variables */
+	NLRowColumn* Ci = NULL ;
+	NLuint       i,j,count ;
+
+	/* Sanity checks */
+	nl_assert(!(M->storage & NL_MATRIX_STORE_SYMMETRIC)) ;
+	nl_assert(M->storage & NL_MATRIX_STORE_ROWS) ;
+	nl_assert(M->m == M->n) ;
+
+	set_default_options(&(context->options)) ;
+
+	switch(nlCurrentContext->solver) {
+	case NL_SUPERLU_EXT: {
+		context->options.ColPerm = NATURAL ;
+	} break ;
+	case NL_PERM_SUPERLU_EXT: {
+		context->options.ColPerm = COLAMD ;
+	} break ;
+	case NL_SYMMETRIC_SUPERLU_EXT: {
+		context->options.ColPerm = MMD_AT_PLUS_A ;
+		context->options.SymmetricMode = YES ;
+	} break ;
+	default: {
+		nl_assert_not_reached ;
+	} break ;
+	}
+
+	StatInit(&(context->stat)) ;
+
+	/*
+	 * Step 1: convert matrix M into SUPERLU compressed column representation
+	 * ----------------------------------------------------------------------
+	 */
+
+	NLint*    xa   = NL_NEW_ARRAY(NLint, n+1) ;
+	NLdouble* a    = NL_NEW_ARRAY(NLdouble, nnz) ;
+	NLint*    asub = NL_NEW_ARRAY(NLint, nnz) ;
+
+	count = 0 ;
+	for(i = 0; i < n; i++) {
+		Ci = &(M->row[i]) ;
+		xa[i] = count ;
+		for(j = 0; j < Ci->size; j++) {
+			a[count]    = Ci->coeff[j].value ;
+			asub[count] = Ci->coeff[j].index ;
+			count++ ;
+		}
+	}
+	xa[n] = nnz ;
+
+	dCreate_CompCol_Matrix(
+		&A, n, n, nnz, a, asub, xa,
+		SLU_NR, /* Row wise     */
+		SLU_D,  /* doubles         */
+		SLU_GE  /* general storage */
+	);
+
+	/*
+	 * Step 2: factorize matrix
+	 * ------------------------
+	 */
+
+	context->perm_c = NL_NEW_ARRAY(NLint, n) ;
+	context->perm_r = NL_NEW_ARRAY(NLint, n) ;
+	NLint* etree    = NL_NEW_ARRAY(NLint, n) ;
+
+	get_perm_c(context->options.ColPerm, &A, context->perm_c) ;
+	sp_preorder(&(context->options), &A, context->perm_c, etree, &AC) ;
+
+	int panel_size = sp_ienv(1) ;
+	int relax = sp_ienv(2) ;
+
+	dgstrf(&(context->options),
+		   &AC,
+		   relax,
+		   panel_size,
+		   etree,
+		   NULL,
+		   0,
+		   context->perm_c,
+		   context->perm_r,
+		   &(context->L),
+		   &(context->U),
+		   &(context->stat),
+		   &info) ;
+
+	/*
+	 * Step 3: cleanup
+	 * ---------------
+	 */
+
+	NL_DELETE_ARRAY(xa) ;
+	NL_DELETE_ARRAY(a) ;
+	NL_DELETE_ARRAY(asub) ;
+	NL_DELETE_ARRAY(etree) ;
+	Destroy_SuperMatrix_Store(&A);
+	Destroy_CompCol_Permuted(&AC);
+	StatFree(&(context->stat));
+
+	return NL_TRUE ;
+}
+
 NLboolean nlSolve_SUPERLU() {
 
-    /* OpenNL Context */
-    NLSparseMatrix* M  = &(nlCurrentContext->M) ;
-    NLdouble* b          = nlCurrentContext->b ;
-    NLdouble* x          = nlCurrentContext->x ;
-
-    /* Compressed Row Storage matrix representation */
-    NLuint    n      = nlCurrentContext->n ;
-    NLuint    nnz    = nlSparseMatrixNNZ(M) ; /* Number of Non-Zero coeffs */
-    NLint*    xa     = NL_NEW_ARRAY(NLint, n+1) ;
-    NLdouble* rhs    = NL_NEW_ARRAY(NLdouble, n) ;
-    NLdouble* a      = NL_NEW_ARRAY(NLdouble, nnz) ;
-    NLint*    asub   = NL_NEW_ARRAY(NLint, nnz) ;
-
-    /* Permutation vector */
-    NLint*    perm_r  = NL_NEW_ARRAY(NLint, n) ;
-    NLint*    perm    = NL_NEW_ARRAY(NLint, n) ;
-
-    /* SuperLU variables */
-    SuperMatrix A, B ; /* System       */
-    SuperMatrix L, U ; /* Inverse of A */
-    NLint info ;       /* status code  */
-    DNformat *vals = NULL ; /* access to result */
-    double *rvals  = NULL ; /* access to result */
-
-    /* SuperLU options and stats */
-    superlu_options_t options ;
-    SuperLUStat_t     stat ;
-
-    /* Temporary variables */
-    NLRowColumn* Ri = NULL ;
-    NLuint         i,jj,count ;
-    
-    /* Sanity checks */
-    nl_assert(!(M->storage & NL_MATRIX_STORE_SYMMETRIC)) ;
-    nl_assert(M->storage & NL_MATRIX_STORE_ROWS) ;
-    nl_assert(M->m == M->n) ;
-    
-    /*
-     * Step 1: convert matrix M into SuperLU compressed column 
-     *   representation.
-     * -------------------------------------------------------
-     */
-
-    count = 0 ;
-    for(i=0; i<n; i++) {
-        Ri = &(M->row[i]) ;
-        xa[i] = count ;
-        for(jj=0; jj<Ri->size; jj++) {
-            a[count]    = Ri->coeff[jj].value ;
-            asub[count] = Ri->coeff[jj].index ;
-            count++ ;
-        }
-    }
-    xa[n] = nnz ;
-
-    /* Save memory for SuperLU */
-    nlSparseMatrixClear(M) ;
-
-
-    /*
-     * Rem: SuperLU does not support symmetric storage.
-     * In fact, for symmetric matrix, what we need 
-     * is a SuperLLt algorithm (SuperNodal sparse Cholesky),
-     * but it does not exist, anybody wants to implement it ?
-     * However, this is not a big problem (SuperLU is just
-     * a superset of what we really need.
-     */
-    dCreate_CompCol_Matrix(
-        &A, n, n, nnz, a, asub, xa, 
-        SLU_NR,              /* Row_wise, no supernode */
-        SLU_D,               /* doubles                */ 
-        SLU_GE               /* general storage        */
-    );
-
-    /* Step 2: create vector */
-    dCreate_Dense_Matrix(
-        &B, n, 1, b, n, 
-        SLU_DN, /* Fortran-type column-wise storage */
-        SLU_D,  /* doubles                          */
-        SLU_GE  /* general                          */
-    );
-            
-
-    /* Step 3: set SuperLU options 
-     * ------------------------------
-     */
-
-    set_default_options(&options) ;
-
-    switch(nlCurrentContext->solver) {
-    case NL_SUPERLU_EXT: {
-        options.ColPerm = NATURAL ;
-    } break ;
-    case NL_PERM_SUPERLU_EXT: {
-        options.ColPerm = COLAMD ;