Parameters handling for solvers

Petsc understands common parameters. Support for different option files
for Petsc and Trilinos. Please see files in examples/MatSolve.

examples/MatSolve/database.txt

+PETSc: petsc_options.txt
+Trilinos_Ifpack: trilinos_ifpack_options.xml
+Trilinos_ML: 
+Trilinos_Aztec: 
+Trilinos_Belos: 
+

examples/MatSolve/main.cpp

@@ -15,7 +15,7 @@ int main(int argc, char ** argv)
 	int rank,procs;
 	if( argc < 3 )
 	{
-		std::cout << "Usage: " << argv[0] << " method_number<0:INNER_ILU2,1:INNER_MLILUC,2:ANI,3:PETSC> matrix.mtx [right_hand_side.rhs] [solver_options.txt]" << std::endl;
+		std::cout << "Usage: " << argv[0] << " method_number<0:INNER_ILU2,1:INNER_MLILUC,2:ANI,3:PETSc,4:Trilinos_Aztec,5:Trilinos_Belos,6:Trilinos_Ifpack,7:Trilinos_ML> matrix.mtx [right_hand_side.rhs] [solver_options.txt]" << std::endl;
 		return -1;
 	}
 	Solver::Type type;
@@ -24,7 +24,11 @@ int main(int argc, char ** argv)
 		case 0: type = Solver::INNER_ILU2; break;
 		case 1: type = Solver::INNER_MLILUC; break;
 		case 2: type = Solver::ANI; break;
-		case 3: type = Solver::PETSC; break;
+		case 3: type = Solver::PETSc; break;
+		case 4: type = Solver::Trilinos_Aztec; break;
+		case 5: type = Solver::Trilinos_Belos; break;
+		case 6: type = Solver::Trilinos_Ifpack; break;
+		case 7: type = Solver::Trilinos_ML; break;
 	}
 	Solver::Initialize(&argc,&argv,argc > 4 ? argv[4] : NULL); // Initialize the linear solver in accordance with args
 	{

examples/MatSolve/trilinos_ifpack_options.xml

+<ParameterList name="top">
+  <ParameterList name="test">
+    <ParameterList name="AztecOO">
+      <Parameter name="Max Iterations" type="int" value="2500" />
+      <Parameter name="Tolerance" type="double" value="1.e-6" />
+      <ParameterList name="AztecOO Settings">
+        <Parameter name="Aztec Solver" type="string" value="GMRES" />
+        <Parameter name="AZ_output" type="int" value="0" />
+      </ParameterList>
+    </ParameterList>
+    <ParameterList name="Ifpack">
+      <Parameter name="Prec Type" type="string" value="ILUT" />
+      <Parameter name="Overlap" type="int" value="2" />
+      <ParameterList name="Ifpack Settings">
+        <Parameter name="fact: ilut level-of-fill" type="double" value="2" />
+        <Parameter name="fact: drop tolerance" type="double" value="1.0e-6" />
+      </ParameterList>
+    </ParameterList>
+  </ParameterList>
+</ParameterList>
+  

examples/Solver/database.txt

+PETSc: petsc_options.txt
+Trilinos_Ifpack: trilinos_ifpack_options.xml
+Trilinos_ML: 
+Trilinos_Aztec: 
+Trilinos_Belos: 
+

examples/Solver/main.cpp

@@ -147,7 +147,7 @@ Storage::real tensor_K_mat_rotOxy(Storage::integer mat)
 int main(int argc,char ** argv)
 {
 	
-	Solver::Initialize(&argc,&argv,"");
+	Solver::Initialize(&argc,&argv,"database.txt");
 #if defined(USE_PARTITIONER)
 	Partitioner::Initialize(&argc,&argv);
 #endif

examples/Solver/petsc_options.txt

+#-pc_type ilu
+#-pc_factor_leves 5
+-ksp_max_it 2500
+-ksp_type bcgsl
+-ksp_rtol 1e-13
+-ksp_atol 1e-7
+-ksp_divtol 1e+100
+-ksp_monitor

examples/Solver/trilinos_ifpack_options.xml

+<ParameterList name="top">
+  <ParameterList name="test">
+    <ParameterList name="AztecOO">
+      <Parameter name="Max Iterations" type="int" value="2500" />
+      <Parameter name="Tolerance" type="double" value="1.e-6" />
+      <ParameterList name="AztecOO Settings">
+        <Parameter name="Aztec Solver" type="string" value="GMRES" />
+        <Parameter name="AZ_output" type="int" value="0" />
+      </ParameterList>
+    </ParameterList>
+    <ParameterList name="Ifpack">
+      <Parameter name="Prec Type" type="string" value="ILUT" />
+      <Parameter name="Overlap" type="int" value="2" />
+      <ParameterList name="Ifpack Settings">
+        <Parameter name="fact: ilut level-of-fill" type="double" value="2" />
+        <Parameter name="fact: drop tolerance" type="double" value="1.0e-6" />
+      </ParameterList>
+    </ParameterList>
+  </ParameterList>
+</ParameterList>
+  

inmost_solver.h

@@ -46,7 +46,7 @@ namespace INMOST
 			Trilinos_Belos, ///< external Solver Belos from Trilinos package, currently without preconditioner
 			Trilinos_ML,    ///< external Solver AztecOO with ML preconditioner
 			Trilinos_Ifpack,///< external Solver AztecOO with Ifpack preconditioner
-			PETSC,          ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/.
+			PETSc,          ///< external Solver PETSc, @see http://www.mcs.anl.gov/petsc/.
 			ANI             ///< external Solver from ANI3D based on ILU2 (sequential Fortran version).
 		};
 
@@ -535,11 +535,18 @@ namespace INMOST
 		/// @see Solver::Solve
 		INMOST_DATA_REAL_TYPE Residual();
 		/// Set the matrix and construct the preconditioner.
+		/// @param A Matrix A in linear problem Ax = b
 		/// @param OldPreconditioner If this parameter is set to true,
 		/// then the previous preconditioner will be used,
 		/// otherwise the new preconditioner will be constructed. 
+		/// Preconditioner will be constructed on call to this function
+		/// for INNER_*, PETSc and ANI packages
+		/// for Trilinos preconditioner will be constructed each time Solver::Solve is called
+		/// Any changes to preconditioner parameters should happen before that point
+		/// If you increase gmres_substep after this point, inner methods most likely will fail
 		void SetMatrix(Matrix & A, bool OldPreconditioner = false);
 		/// Solver the linear system: A*x = b.
+		/// Prior to this call you should call SetMatrix
 		/// @param RHS The right-hand side Vector b.
 		/// @param SOL The initial guess to the solution on input and the solution Vector x on return.
 		/// It is assumed that the coefficient matrix A have been set
@@ -551,6 +558,8 @@ namespace INMOST
 		std::string GetReason();
 
 		/// Main constructor of the solver.
+		/// Solver name provided here is used to extract options from database file
+		/// for PETSc and Trilinos packages.
 		/// @param pack The package Type to be used for solution.
 		/// @param _name The user specified name of the current solver.
 		/// @param _comm Communicator for parallel data exchanges, MPI_COMM_WORLD by default.
@@ -562,12 +571,23 @@ namespace INMOST
 		~Solver();
 		/// Initialize the stage of parallel solution.
 		/// If MPI is not initialized yet, then it will be initialized.
+		/// database file is used to pass parameters to PETSc and Trilinos packages.
+		/// if database file for is provided any changes through SetParameterEnum,
+		/// SetParameterReal would not be effective for PETSc and Trilinos packages.
+		/// Currently this database file provides directions for package-specific
+		/// files. In future it is supposed to set up parameters for internal solvers.
 		/// @param argc The number of arguments transmitted to the function main.
 		/// @param argv The pointer to arguments transmitted to the function main.
 		/// @param database Usually the name of the file with the Solver parameters.
 		/// The shortest call to this function with the default solver parameters is the following: Initialize(NULL,NULL,"");
 		/// @see Solver::Finalize
 		/// @see Solver::isInitialized
+		/// Example of contents of the database file:
+		/// PETSc: petsc_options.txt
+		/// Trilinos_Ifpack: trilinos_ifpack_options.xml
+		/// Trilinos_ML: trilinos_ml_options.xml
+		/// Trilinos_Aztec: trilinos_aztec_options.xml
+		/// Trilinos_Belos: trilinos_belos_options.xml
 		static void Initialize(int * argc, char *** argv, const char * database = "");
 		/// Finalize the stage of parallel solution.
 		/// If MPI was initialized in Solver::Initialize, then it will be finalized.

solver_bcgsl.hpp

@@ -48,7 +48,11 @@ namespace INMOST
 				if (prec != NULL) return prec->EnumParameter(name.substr(1, name.size() - 1));
 			}
 			if (name == "maxits") return maxits;
-			else if (name == "levels") return l;
+			else if (name == "levels") 
+			{
+				if( init ) throw - 1; //solver was already initialized, value should not be changed
+				return l;
+			}
 			else if (prec != NULL) return prec->EnumParameter(name);
 			throw - 1;
 		}

solver_petsc.cpp

@@ -326,6 +326,57 @@ double SolverResidualNormPetsc(void * data)
 	return norm;
 }
 
+
+void SolverSetTolerancesPetsc(void * data, double rtol, double atol, double divtol, int maxits)
+{
+	KSP * ksp = static_cast<KSP *>(data);
+	PetscErrorCode ierr = KSPSetTolerances(*ksp,rtol,atol,divtol,maxits);
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+}
+
+
+void SolverSetOverlapPetsc(void * data, int levels)
+{
+	KSP * ksp = static_cast<KSP *>(data);
+	PetscErrorCode ierr;
+	PC pc; PCType pc_type;
+	ierr = KSPGetPC(*ksp,&pc);
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+	ierr = PCGetType(pc,&pc_type);
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+	if( !strcmp(pc_type,"asm") )
+	{
+		ierr = PCASMSetOverlap(pc,levels);
+	}
+	else if( !strcmp(pc_type,"gasm") )
+	{
+		ierr = PCGASMSetOverlap(pc,levels);
+	}
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+}
+
+void SolverSetDropTolerancePetsc(void * data, double dtol)
+{
+	KSP * ksp = static_cast<KSP *>(data);
+	PetscErrorCode ierr;
+	PC pc;
+	ierr = KSPGetPC(*ksp,&pc);
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+	ierr = PCFactorSetDropTolerance(pc,dtol,0.01,10000); //2 other parameters are set to extreme
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+}
+
+void SolverSetFillLevelPetsc(void * data, double lfill)
+{
+	KSP * ksp = static_cast<KSP *>(data);
+	PetscErrorCode ierr;
+	PC pc;
+	ierr = KSPGetPC(*ksp,&pc);
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+	ierr = PCFactorSetLevels(pc,static_cast<PetscInt>(lfill));
+	if( ierr != PETSC_SUCCESS ) throw INMOST::ErrorInSolver;
+}
+
 const char * SolverConvergedReasonPetsc(void * data)
 {
 	static char reason_str[4096];

solver_petsc.h

@@ -36,4 +36,9 @@ int SolverIterationNumberPetsc(void * data);
 double SolverResidualNormPetsc(void * data);
 const char * SolverConvergedReasonPetsc(void * data);
 
+void SolverSetTolerancesPetsc(void * data, double rtol, double atol, double divtol, int maxits);
+void SolverSetOverlapPetsc(void * data, int levels);
+void SolverSetDropTolerancePetsc(void * data, double dtol);
+void SolverSetFillLevelPetsc(void * data, double lfill);
+
 #endif