Preprocessing of singular matrices for SuperLU

Examples/ADMFD/main.cpp

@@ -556,8 +556,8 @@ int main(int argc,char ** argv)
                 
                 if( R.Norm() < 1.0e-4 ) break;
                 
-                Solver S(Solver::INNER_MPTILUC);
-				//Solver S(Solver::SUPERLU);
+                //Solver S(Solver::INNER_MPTILUC);
+				Solver S(Solver::SUPERLU);
 				S.SetParameterReal("relative_tolerance", 1.0e-14);
                 S.SetParameterReal("absolute_tolerance", 1.0e-12);
                 S.SetParameterReal("drop_tolerance", 1.0e-2);

Source/Solver/solver.cpp

@@ -1752,24 +1752,50 @@ namespace INMOST
 			double * a;
 			int mbeg = A.GetFirstIndex();
 			int mend = A.GetLastIndex();
-			for(int k = 0; k < mend-mbeg; ++k) nnz += A[k].Size();
-			ia = (int *)malloc(sizeof(int)*(A.Size()+1));
+			int size = 0;
+			int * remap = new int[mend-mbeg];
+			for(int k = 0; k < mend-mbeg; ++k) 
+			{
+				Sparse::Row & r = A[k+mbeg];
+				if( r.Size() )
+				{
+					double nrm = 0;
+					for(int l = 0; l < r.Size(); ++l)
+						nrm += r.GetValue(l)*r.GetValue(l);
+					if( nrm )
+					{
+						std::sort(r.Begin(),r.End());
+						nnz += r.Size();
+						remap[k] = size;
+						size++;
+					}
+					else remap[k] = -1;
+				}
+				else remap[k] = -1;
+			}
+			ia = (int *)malloc(sizeof(int)*(size+1));
 			ja = (int *)malloc(sizeof(int)*nnz);
 			a = (double *)malloc(sizeof(double)*nnz);
-			int q = 0;
+			int q = 0, f = 0;
 			ia[0] = 0;
-			for(int k = 0; k < mend-mbeg; ++k)
+			for(int k = 0; k < mend-mbeg; ++k) if( remap[k] != -1 )
 			{
 				Sparse::Row & r = A[k+mbeg];
 				for(int l = 0; l < r.Size(); ++l)
 				{
-					ja[q] = r.GetIndex(l)-mbeg;
-					a[q] = r.GetValue(l);
-					++q;
+					if( remap[r.GetIndex(l)-mbeg] != -1 )
+					{
+						ja[q] = remap[r.GetIndex(l)-mbeg];
+						a[q] = r.GetValue(l);
+						++q;
+					}
+					else //if( fabs(a[q]) > 1.0e-9 )
+						std::cout << "Matrix has connections to singular rows" << std::endl;
 				}
-				ia[k+1] = q;
+				ia[f+1] = q;
+				f++;
 			}
-			MatrixFillSuperLU(solver_data,mend-mbeg,nnz,ia,ja,a);
+			MatrixFillSuperLU(solver_data,size,nnz,ia,ja,a,remap);
 			//arrays are freed inside SuperLU
 			//delete [] ia;
 			//delete [] ja;
@@ -2554,7 +2580,14 @@ namespace INMOST
 #if defined(USE_SOLVER_SUPERLU)
 		if(_pack == SUPERLU )
 		{
-			bool ret = SolverSolveSuperLU(solver_data,&*RHS.Begin(),&*SOL.Begin());
+			int size = MatrixSizeSuperLU(solver_data);
+			double * inout = new double[size];
+			int * remap = MatrixRemapArraySuperLU(solver_data);
+			int mbeg = RHS.GetFirstIndex(), mend = RHS.GetLastIndex();
+			for(int k = 0; k < mend - mbeg; ++k) if( remap[k] != -1 ) inout[remap[k]] = RHS[k+mbeg];
+			bool ret = SolverSolveSuperLU(solver_data,inout);
+			for(int k = 0; k < mend - mbeg; ++k) if( remap[k] != -1 ) SOL[k+mbeg] = inout[remap[k]];
+			delete [] inout;
 			last_it = 1;
 			last_resid = 0;
 			return_reason = SolverConvergedReasonSuperLU(solver_data);

Source/Solver/solver_superlu.cpp

@@ -10,6 +10,7 @@ struct SuperLU_data
 	SuperMatrix A, L, U;
 	int * perm_r;
 	int * perm_c;
+	int * remap;
 	int size, info;
 	superlu_options_t options;
 	SuperLUStat_t stat;
@@ -21,15 +22,26 @@ void SolverInitializeSuperLU(int * argc,char *** argv, const char * file_options
 }
 
 
-void MatrixFillSuperLU(void * data, int size, int nnz, int * col_ranges, int * col_positions, double * col_values)
+void MatrixFillSuperLU(void * data, int size, int nnz, int * col_ranges, int * col_positions, double * col_values, int * remap)
 {
 	SuperLU_data * m = static_cast<SuperLU_data *>(data);
 	dCreate_CompCol_Matrix(&m->A,size,size,nnz,col_values,col_positions,col_ranges,SLU_NR,SLU_D,SLU_GE);
+	//dPrint_CompCol_Matrix("A", &m->A);
 	m->size = size;
+	m->remap = remap;
 	m->perm_c = new int [size];
 	m->perm_r = new int [size];
 }
 
+int * MatrixRemapArraySuperLU(void * data)
+{
+	return static_cast<SuperLU_data *>(data)->remap;
+}
+
+int MatrixSizeSuperLU(void * data)
+{
+	return static_cast<SuperLU_data *>(data)->size;
+}
 
 void SolverDestroyDataSuperLU(void ** data)
 {
@@ -40,6 +52,7 @@ void SolverDestroyDataSuperLU(void ** data)
 	StatFree(&m->stat);
 	if( m->perm_c != NULL ) delete [] m->perm_c;
 	if( m->perm_r != NULL ) delete [] m->perm_r;
+	if( m->remap != NULL ) delete [] m->remap; //allocated outside
 	delete m;
 	*data = NULL;
 }
@@ -59,13 +72,12 @@ void SolverInitDataSuperLU(void ** data)
 }
 
 
-bool SolverSolveSuperLU(void * data, void * rhs_data, void * sol_data)
+bool SolverSolveSuperLU(void * data, void * rhs_in_sol_out_data)
 {
 	SuperLU_data * m = static_cast<SuperLU_data *>(data);
 	assert(m->size != 0);
 	SuperMatrix B;
-	memcpy(sol_data,rhs_data,sizeof(double)*m->size);
-	dCreate_Dense_Matrix(&B,m->size,1,(double *)sol_data,m->size,SLU_DN,SLU_D,SLU_GE);
+	dCreate_Dense_Matrix(&B,m->size,1,(double *)rhs_in_sol_out_data,m->size,SLU_DN,SLU_D,SLU_GE);
 	dgssv(&m->options,&m->A,m->perm_c,m->perm_r,&m->L,&m->U,&B,&m->stat,&m->info);
 	Destroy_SuperMatrix_Store(&B);
 	return m->info == 0;

Source/Solver/solver_superlu.h

@@ -4,10 +4,12 @@
 //IF CERTAIN FUNCTIONALITY IS NOT AVAILABLE, PLEASE THROW INMOST::NotImplemented EXCEPTION
 
 void SolverInitializeSuperLU(int * argc,char *** argv, const char * file_options);
-void MatrixFillSuperLU(void * data, int size, int nnz, int * col_ranges, int * col_positions, double * col_values);
+void MatrixFillSuperLU(void * data, int size, int nnz, int * col_ranges, int * col_positions, double * col_values,int * remap);
+int * MatrixRemapArraySuperLU(void * data);
+int MatrixSizeSuperLU(void * data);
 void SolverDestroyDataSuperLU(void ** data);
 void SolverInitDataSuperLU(void ** data);
-bool SolverSolveSuperLU(void * data, void * rhs_data, void * sol_data);
+bool SolverSolveSuperLU(void * data, void * rhs_in_sol_out_data);
 const char * SolverConvergedReasonSuperLU(void * data);
 
 #endif