9 #include "lilc_matrix.h" 54 template<
class el_type>
55 bool save_vector(
const std::vector<el_type>& vec, std::string filename) {
56 std::ofstream out(filename.c_str(), std::ios::out | std::ios::binary);
60 out.flags(std::ios_base::scientific);
62 std::string header =
"%%MatrixMarket matrix coordinate real general";;
64 out << header << std::endl;
65 out << vec.size() <<
" " << 1 <<
" " << vec.size() <<
"\n";
67 for(
int i = 0; i < (int) vec.size(); i++) {
68 out << i+1 <<
" " << 1 <<
" " << vec[i] <<
"\n";
79 template<
class el_type>
80 bool read_vector(std::vector<el_type>& vec, std::string filename,
int msg_lvl = message_level::STATISTICS) {
81 std::ifstream input(filename.c_str(), std::ios::in);
83 if(!input)
return false;
85 const int maxBuffersize = 2048;
86 char buffer[maxBuffersize];
88 bool readsizes =
false;
91 int i = 0, n_rows, n_cols;
92 while(input.getline(buffer, maxBuffersize)) {
95 if(buffer[0]==
'%')
continue;
97 std::stringstream line(buffer);
100 line >> n_rows >> n_cols;
101 if(n_rows > 0 && n_cols > 0) {
103 vec.resize(std::max(n_rows, n_cols));
111 if (i != std::max(n_rows, n_cols)) {
112 std::cerr <<
"Expected " << std::max(n_rows, n_cols) <<
" elems but read " << i <<
"." << std::endl;
115 if (msg_lvl) std::cout <<
"Load succeeded. " <<
"Vector file " << filename <<
" was loaded." << std::endl;
124 template<
class el_type,
class mat_type = lilc_matrix<el_type> >
127 typedef typename mat_type::pivot_type pivot_type;
152 msg_lvl = message_level::STATISTICS;
153 piv_type = pivot_type::ROOK;
154 reorder_type = reordering_type::AMD;
155 equil_type = equilibration_type::BUNCH;
156 solve_type = solver_type::SQMR;
158 perform_inplace =
false;
164 void load(std::string filename) {
165 bool result = A.load(filename);
167 if (msg_lvl) printf(
"A is %d by %d with %d non-zeros.\n", A.n_rows(), A.n_cols(), A.nnz() );
172 void load(
const std::vector<int>& ptr,
const std::vector<int>& row,
const std::vector<el_type>& val) {
173 bool result = A.load(ptr, row, val);
175 if (msg_lvl) printf(
"A is %d by %d with %d non-zeros.\n", A.n_rows(), A.n_cols(), A.nnz() );
180 void load(
const int* ptr,
const int* row,
const el_type* val,
int dim) {
181 bool result = A.load(ptr, row, val, dim);
183 if (msg_lvl) printf(
"A is %d by %d with %d non-zeros.\n", A.n_rows(), A.n_cols(), A.nnz() );
193 if (msg_lvl) printf(
"Right hand side has %d entries.\n", rhs.size() );
199 if (strcmp(ordering,
"rcm") == 0) {
200 reorder_type = reordering_type::RCM;
201 }
else if (strcmp(ordering,
"amd") == 0) {
202 reorder_type = reordering_type::AMD;
203 }
else if (strcmp(ordering,
"none") == 0) {
204 reorder_type = reordering_type::NONE;
211 if (strcmp(equil,
"bunch") == 0) {
212 equil_type = equilibration_type::BUNCH;
213 }
else if (strcmp(equil,
"none") == 0) {
214 equil_type = equilibration_type::NONE;
221 if (strcmp(solver,
"minres") == 0) {
222 solve_type = solver_type::MINRES;
223 }
else if (strcmp(solver,
"sqmr") == 0) {
224 solve_type = solver_type::SQMR;
225 }
else if (strcmp(solver,
"full") == 0) {
226 solve_type = solver_type::FULL;
227 }
else if (strcmp(solver,
"none") == 0) {
228 solve_type = solver_type::NONE;
235 if (strcmp(msg,
"none") == 0) {
236 msg_lvl = message_level::NONE;
237 }
else if (strcmp(msg,
"statistics") == 0) {
238 msg_lvl = message_level::STATISTICS;
239 }
else if (strcmp(msg,
"debug") == 0) {
240 msg_lvl = message_level::DEBUG;
247 perform_inplace = inplace;
253 if (strcmp(pivot,
"rook") == 0) {
254 piv_type = pivot_type::ROOK;
255 }
else if (strcmp(pivot,
"bunch") == 0) {
256 piv_type = pivot_type::BKP;
269 void solve(
double fill_factor,
double tol,
double pp_tol,
int max_iter = -1,
double minres_tol = 1e-6,
double shift = 0.0) {
271 if (solve_type == solver_type::FULL) {
273 fill_factor = A.n_cols();
276 perm.reserve(A.n_cols());
277 cout << std::fixed << std::setprecision(3);
279 double dif, total = 0;
282 if (equil_type == equilibration_type::BUNCH) {
285 dif = clock() - start; total += dif;
286 if (msg_lvl) printf(
" Equilibration:\t\t%.3f seconds.\n", dif/CLOCKS_PER_SEC);
289 if (reorder_type != reordering_type::NONE) {
291 std::string perm_name;
292 switch (reorder_type) {
293 case reordering_type::AMD:
297 case reordering_type::RCM:
303 dif = clock() - start; total += dif;
304 if (msg_lvl) printf(
" %s:\t\t\t\t%.3f seconds.\n", perm_name.c_str(), dif/CLOCKS_PER_SEC);
308 dif = clock() - start; total += dif;
309 if (msg_lvl) printf(
" Permutation:\t\t\t%.3f seconds.\n", dif/CLOCKS_PER_SEC);
312 for (
int i = 0; i < A.n_cols(); i++) {
318 if (perform_inplace) {
319 A.ildl_inplace(D, perm, fill_factor, tol, pp_tol, piv_type);
321 A.ildl(L, D, perm, fill_factor, tol, pp_tol, piv_type);
323 dif = clock() - start; total += dif;
325 std::string pivot_name;
326 if (piv_type == pivot_type::BKP) {
328 }
else if (piv_type == pivot_type::ROOK) {
332 if (msg_lvl) printf(
" Factorization (%s pivoting):\t%.3f seconds.\n", pivot_name.c_str(), dif/CLOCKS_PER_SEC);
333 if (msg_lvl) printf(
"Total time:\t\t\t%.3f seconds.\n", total/CLOCKS_PER_SEC);
334 if (perform_inplace) {
335 if (msg_lvl) printf(
"L is %d by %d with %d non-zeros.\n", A.n_rows(), A.n_cols(), A.nnz() );
337 if (msg_lvl) printf(
"L is %d by %d with %d non-zeros.\n", L.n_rows(), L.n_cols(), L.nnz() );
339 if (msg_lvl) printf(
"\n");
346 if (perform_inplace) {
347 if (msg_lvl) printf(
"Inplace factorization cannot be used with the solver. Please try again without -inplace.\n");
356 for (
int i = 0; i < A.n_cols(); i++) {
357 rhs[i] = A.S[i]*rhs[i];
362 vector<el_type> tmp(A.n_cols());
363 for (
int i = 0; i < A.n_cols(); i++) {
364 tmp[i] = rhs[perm[i]];
368 if (solve_type == solver_type::FULL) {
369 if (msg_lvl) printf(
"Solving matrix with direct solver...\n");
370 sol_vec.resize(A.n_cols(), 0);
374 L.backsolve(rhs, sol_vec);
375 D.
solve(sol_vec, tmp);
376 L.forwardsolve(tmp, sol_vec);
380 if (solve_type == solver_type::MINRES) {
384 L.backsolve(rhs, tmp);
387 if (msg_lvl) printf(
"Solving matrix with MINRES...\n");
389 minres(max_iter, minres_tol, shift);
399 L.forwardsolve(tmp, sol_vec);
400 }
else if (solve_type == solver_type::SQMR) {
401 if (msg_lvl) printf(
"Solving matrix with SQMR...\n");
402 sqmr(max_iter, minres_tol);
407 for (
int i = 0; i < A.n_cols(); i++) {
408 tmp[perm[i]] = sol_vec[i];
413 for (
int i = 0; i < A.n_cols(); i++) {
414 sol_vec[i] = A.S[i]*sol_vec[i];
416 dif = clock() - start;
417 if (msg_lvl) printf(
"Solve time:\t%.3f seconds.\n", dif/CLOCKS_PER_SEC);
418 if (msg_lvl) printf(
"\n");
423 if (msg_lvl) printf(
"Solution saved to output_matrices/outsol.mtx.\n");
424 save_vector(sol_vec,
"output_matrices/outsol.mtx");
437 void minres(
int max_iter = 1000,
double stop_tol = 1e-6,
double shift = 0.0);
444 void sqmr(
int max_iter = 1000,
double stop_tol = 1e-6);
451 if (msg_lvl) cout <<
"Saving matrices..." << endl;
452 if (!perform_inplace) {
453 A.save(
"output_matrices/outB.mtx",
true);
454 L.save(
"output_matrices/outL.mtx",
false);
456 A.save(
"output_matrices/outL.mtx",
false);
459 A.S.save(
"output_matrices/outS.mtx");
460 save_vector(perm,
"output_matrices/outP.mtx");
462 D.
save(
"output_matrices/outD.mtx");
463 if (msg_lvl) cout <<
"Save complete." << endl;
469 #ifdef SYM_ILDL_DEBUG 470 if (perform_inplace) {
476 cout << perm << endl;
481 #include "solver_minres.h" 482 #include "solver_sqmr.h" mat_type A
The matrix to be factored.
Definition: solver.h:129
void set_reorder_scheme(const char *ordering)
Sets the reordering scheme for the solver.
Definition: solver.h:198
bool perform_inplace
Set to true if we are factoring the matrix A inplace.
Definition: solver.h:142
bool has_rhs
Set to true if we have a right hand side that we expect to solve.
Definition: solver.h:141
void display()
Prints the L and D factors to stdout.
Definition: solver.h:468
void solve(const elt_vector_type &b, elt_vector_type &x)
Solves the system Dx = b.
Definition: block_diag_matrix.h:199
void save()
Save results of factorization (automatically saved into the output_matrices folder).
Definition: solver.h:450
block_diag_matrix< el_type > D
The diagonal factor of A.
Definition: solver.h:133
void sqrt_solve(const elt_vector_type &b, elt_vector_type &x, bool transposed=false)
Solves the preconditioned problem |D| = Q|V|Q', where QVQ' is the eigendecomposition of D...
Definition: block_diag_matrix.h:143
void load(const std::vector< int > &ptr, const std::vector< int > &row, const std::vector< el_type > &val)
Loads the matrix A into solver. A must be of CSC format.
Definition: solver.h:172
void set_pivot(const char *pivot)
Decides the kind of partial pivoting we should use.
Definition: solver.h:252
void set_message_level(const char *msg)
Controls how much information gets printed to stdout.
Definition: solver.h:234
bool save(std::string filename) const
Definition: block_diag_matrix_save.h:6
void load(const int *ptr, const int *row, const el_type *val, int dim)
Loads the matrix A into solver. A must be of CSC format.
Definition: solver.h:180
void set_equil(const char *equil)
Decides whether we should use equilibriation on the matrix or not.
Definition: solver.h:210
int piv_type
Set to 0 for rook, 1 for bunch.
Definition: solver.h:135
void load(std::string filename)
Loads the matrix A into solver. A must be of matrix market format.
Definition: solver.h:164
int reorder_type
Set to to 0 for AMD, 1 for RCM, 2 for no reordering.
Definition: solver.h:134
A quick implementation of a diagonal matrix with 1x1 and 2x2 blocks.
Definition: block_diag_matrix.h:39
int equil_type
The equilibration method used. Set to 1 for max-norm equilibriation.
Definition: solver.h:137
solver()
Solver constructor, initializes default reordering scheme.
Definition: solver.h:151
bool save_sol
Set to true if we want to save the solution to a file.
Definition: solver.h:144
void set_inplace(bool inplace)
Decides whether we perform the factorization inplace or not.
Definition: solver.h:246
Set of tools that facilitates conversion between different matrix formats. Also contains solver metho...
Definition: solver.h:125
mat_type L
The lower triangular factor of A.
Definition: solver.h:130
void set_solver(const char *solver)
Decides whether we perform a full solve or not.
Definition: solver.h:220
vector< el_type > rhs
The right hand side we'll solve for.
Definition: solver.h:146
void set_rhs(vector< el_type > b)
Loads a right hand side b into the solver.
Definition: solver.h:190
vector< el_type > sol_vec
The solution vector.
Definition: solver.h:147
vector< int > perm
A permutation vector containing all permutations on A.
Definition: solver.h:132
void solve(double fill_factor, double tol, double pp_tol, int max_iter=-1, double minres_tol=1e-6, double shift=0.0)
Factors the matrix A into P' * S * A * S * P = LDL' in addition to printing some timing data to scree...
Definition: solver.h:269
int msg_lvl
Controls the amount of output to stdout.
Definition: solver.h:139