SolverILU2.cpp 4.71 KB
Newer Older
Dmitry Bagaev's avatar
Dmitry Bagaev committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
#include "SolverILU2.h"

namespace INMOST {

	    SolverILU2::SolverILU2() {
	    	additive_schwartz_overlap = DEFAULT_ADDITIVE_SCHWARTZ_OVERLAP;
			maximum_iterations = DEFAULT_MAXIMUM_ITERATIONS;
			absolute_tolerance = DEFAULT_ABSOLUTE_TOLERANCE;
			relative_tolerance = DEFAULT_RELATIVE_TOLERANCE;
			divergence_tolerance = DEFAULT_DIVERGENCE_TOLERANCE;

			preconditioner_drop_tolerance = DEFAULT_PRECONDITIONER_DROP_TOLERANCE;
			preconditioner_reuse_tolerance = DEFAULT_PRECONDITIONER_REUSE_TOLERANCE;
			preconditioner_rescale_iterations = DEFAULT_PRECONDITIONER_RESCALE_ITERS;
			preconditioner_fill_level = DEFAULT_PRECONDITIONER_FILL_LEVEL;

	    	Method *prec = new ILU2_preconditioner(info);
	    	solver = new BCGS_solver(prec, info);
	    	matrix = NULL;
	    }

        SolverILU2::SolverILU2(const SolverInterface* other) {
        	//You should not really want to copy solver's information
        	throw INMOST::SolverUnsupportedOperation;
        }

		void SolverILU2::Assign(const SolverInterface* other) {
			//You should not really want to copy solver's information
        	throw INMOST::SolverUnsupportedOperation;
		}

        void SolverILU2::Initialize(int *argc, char ***argv, const char *parameters_file, std::string prefix) {

        }

		void SolverILU2::SetMatrix(Sparse::Matrix & A, bool ModifiedPattern, bool OldPreconditioner) {
			Sparse::Matrix *m = new Sparse::Matrix(A);
			info.PrepareMatrix(*m, additive_schwartz_overlap);
			solver->ReplaceMAT(*m);
			if (matrix != NULL) {
				delete matrix;
			}
			matrix = m;

			solver->RealParameter(":tau") = preconditioner_drop_tolerance;
			solver->RealParameter(":tau2") = preconditioner_reuse_tolerance;
			solver->EnumParameter(":scale_iters") = preconditioner_rescale_iterations;
			solver->EnumParameter(":fill") = static_cast<INMOST_DATA_ENUM_TYPE>(preconditioner_fill_level);

			if (!solver->isInitialized()) {
				solver->Initialize();
			}
		}

        bool SolverILU2::Solve(Sparse::Vector &RHS, Sparse::Vector &SOL) {
			solver->EnumParameter("maxits") = maximum_iterations;
			solver->RealParameter("rtol") = relative_tolerance;
			solver->RealParameter("atol") = absolute_tolerance;
			solver->RealParameter("divtol") = divergence_tolerance;
			
			return solver->Solve(RHS,SOL);
        }

        bool SolverILU2::Clear() {
        	info.Clear();
        	if (matrix != NULL) {
        		delete matrix;
        		matrix = NULL;
        	}
        	if (solver != NULL) {
        		delete solver;
        		solver = NULL;
        	}
        	return true;
        }

    	bool SolverILU2::isMatrixSet() {
    		return matrix != NULL;
    	}

		INMOST_DATA_REAL_TYPE SolverILU2::GetPropertyReal(std::string property) const {
			return solver->RealParameter(property);
		}

		INMOST_DATA_ENUM_TYPE SolverILU2::GetPropertyEnum(std::string property) const {
			return solver->EnumParameter(property);
		}

		void SolverILU2::SetPropertyReal(std::string property, INMOST_DATA_REAL_TYPE value) {
			//Maybe we should use explicit names?, like maxits ..etc
        	//like solver->RealParameter(property) = value;
			if (property == "absolute_tolerance") {
				absolute_tolerance = value;	
			} else if (property == "relative_tolerance") {
				relative_tolerance = value;
			} else if (property == "divergence_tolerance") {
				divergence_tolerance = value;
			} else if (property == "drop_tolerance") {
				preconditioner_drop_tolerance = value;
			} else if (property == "reuse_tolerance") {
				preconditioner_reuse_tolerance = value;
			} else if (property == "fill_level") {
				preconditioner_fill_level = value;
			} else {
				std::cout << "Parameter " << property << " of real type is unknown" << std::endl;
			}
		}

        void SolverILU2::SetPropertyEnum(std::string property, INMOST_DATA_ENUM_TYPE value) {
        	//Maybe we should use explicit names?, like maxits ..etc
        	//like solver->EnumParameter(property) = value;
        	if (property == "maximum_iterations") {
				maximum_iterations = value;
        	} else if (property == "rescale_iterations") {
        		preconditioner_rescale_iterations = value;	
        	} else if (property == "schwartz_overlap") {
        		additive_schwartz_overlap = value;	
        	} else std::cout << "Parameter " << property << " of integral type is unknown" << std::endl;
        }

		const INMOST_DATA_ENUM_TYPE SolverILU2::Iterations() const {
			return solver->GetIterations();
		}

		const INMOST_DATA_REAL_TYPE SolverILU2::Residual() const {
			return solver->GetResidual();
		}

		const std::string SolverILU2::ReturnReason() const {
			return solver->GetReason();
		}

		const std::string SolverILU2::SolverName() const {
			return "inner_ilu2";
		}

		void SolverILU2::Finalize() {

		}

        SolverILU2::~SolverILU2() {
        	Clear();
        }

}