sym-ildl  1.2
Incomplete LDL' factorizations of indefinite symmetric and skew-symmetric matrices.
solver.h
1 #ifndef _SOLVER_H
2 #define _SOLVER_H
3 
4 #include <iostream>
5 #include <cstring>
6 #include <ctime>
7 #include <iomanip>
8 
9 #include "lilc_matrix.h"
10 
11 namespace symildl {
12 
13 // Using struct'd enums to achieve a C++11 style enum class without C++11
15  enum {
16  NONE,
17  AMD,
18  RCM,
19  MC64
20  };
21 };
22 
24  enum {
25  NONE,
26  BUNCH,
27  RUIZ,
28  MC64
29  };
30 };
31 
32 struct solver_type {
33  enum {
34  NONE,
35  MINRES,
36  SQMR,
37  FULL
38  };
39 };
40 
41 struct message_level {
42  enum {
43  NONE,
44  STATISTICS,
45  DEBUG
46  };
47 };
48 
49 
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);
57  if(!out)
58  return false;
59 
60  out.flags(std::ios_base::scientific);
61  out.precision(12);
62  std::string header = "%%MatrixMarket matrix coordinate real general";;
63 
64  out << header << std::endl;
65  out << vec.size() << " " << 1 << " " << vec.size() << "\n";
66 
67  for(int i = 0; i < (int) vec.size(); i++) {
68  out << i+1 << " " << 1 << " " << vec[i] << "\n";
69  }
70 
71  out.close();
72  return true;
73 }
74 
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);
82 
83  if(!input) return false;
84 
85  const int maxBuffersize = 2048;
86  char buffer[maxBuffersize];
87 
88  bool readsizes = false;
89  el_type value;
90 
91  int i = 0, n_rows, n_cols;
92  while(input.getline(buffer, maxBuffersize)) {
93  // skip comments
94  //NOTE An appropriate test should be done on the header to get the symmetry
95  if(buffer[0]=='%') continue;
96 
97  std::stringstream line(buffer);
98 
99  if(!readsizes) {
100  line >> n_rows >> n_cols;
101  if(n_rows > 0 && n_cols > 0) {
102  readsizes = true;
103  vec.resize(std::max(n_rows, n_cols));
104  }
105  } else {
106  line >> value;
107  vec[i++] = value;
108  }
109  }
110 
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;
113  }
114 
115  if (msg_lvl) std::cout << "Load succeeded. " << "Vector file " << filename << " was loaded." << std::endl;
116  input.close();
117  return true;
118 }
119 
124 template<class el_type, class mat_type = lilc_matrix<el_type> >
125 class solver {
126  public:
127  typedef typename mat_type::pivot_type pivot_type;
128 
129  mat_type A;
130  mat_type L;
131 
132  vector<int> perm;
135  int piv_type;
136 
138 
139  int msg_lvl;
140  int solve_type; //<The type of solver used to solve the right hand side.
141  bool has_rhs;
143  // TODO: refactor this away
144  bool save_sol;
145 
146  vector<el_type> rhs;
147  vector<el_type> sol_vec;
148 
151  solver() {
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;
157  has_rhs = false;
158  perform_inplace = false;
159  }
160 
164  void load(std::string filename) {
165  bool result = A.load(filename);
166  assert(result);
167  if (msg_lvl) printf("A is %d by %d with %d non-zeros.\n", A.n_rows(), A.n_cols(), A.nnz() );
168  }
169 
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);
174  assert(result);
175  if (msg_lvl) printf("A is %d by %d with %d non-zeros.\n", A.n_rows(), A.n_cols(), A.nnz() );
176  }
177 
180  void load(const int* ptr, const int* row, const el_type* val, int dim) {
181  bool result = A.load(ptr, row, val, dim);
182  assert(result);
183  if (msg_lvl) printf("A is %d by %d with %d non-zeros.\n", A.n_rows(), A.n_cols(), A.nnz() );
184  }
185 
186 
190  void set_rhs(vector<el_type> b) {
191  rhs = b;
192  has_rhs = true;
193  if (msg_lvl) printf("Right hand side has %d entries.\n", rhs.size() );
194  }
195 
198  void set_reorder_scheme(const char* ordering) {
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;
205  }
206  }
207 
210  void set_equil(const char* equil) {
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;
215  }
216  }
217 
220  void set_solver(const char* solver) {
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;
229  }
230  }
231 
234  void set_message_level(const char* msg) {
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;
241  }
242  }
243 
246  void set_inplace(bool inplace) {
247  perform_inplace = inplace;
248  }
249 
252  void set_pivot(const char* pivot) {
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;
257  }
258  }
259 
269  void solve(double fill_factor, double tol, double pp_tol, int max_iter = -1, double minres_tol = 1e-6, double shift = 0.0) {
270  // A full factorization is equivalent to a fill factor of n and tol of 0
271  if (solve_type == solver_type::FULL) {
272  tol = 0.0;
273  fill_factor = A.n_cols();
274  }
275 
276  perm.reserve(A.n_cols());
277  cout << std::fixed << std::setprecision(3);
278 
279  double dif, total = 0;
280  clock_t start;
281 
282  if (equil_type == equilibration_type::BUNCH) {
283  start = clock();
284  A.sym_equil();
285  dif = clock() - start; total += dif;
286  if (msg_lvl) printf(" Equilibration:\t\t%.3f seconds.\n", dif/CLOCKS_PER_SEC);
287  }
288 
289  if (reorder_type != reordering_type::NONE) {
290  start = clock();
291  std::string perm_name;
292  switch (reorder_type) {
293  case reordering_type::AMD:
294  A.sym_amd(perm);
295  perm_name = "AMD";
296  break;
297  case reordering_type::RCM:
298  A.sym_rcm(perm);
299  perm_name = "RCM";
300  break;
301  }
302 
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);
305 
306  start = clock();
307  A.sym_perm(perm);
308  dif = clock() - start; total += dif;
309  if (msg_lvl) printf(" Permutation:\t\t\t%.3f seconds.\n", dif/CLOCKS_PER_SEC);
310  } else {
311  // no permutation specified, store identity permutation instead.
312  for (int i = 0; i < A.n_cols(); i++) {
313  perm.push_back(i);
314  }
315  }
316 
317  start = clock();
318  if (perform_inplace) {
319  A.ildl_inplace(D, perm, fill_factor, tol, pp_tol, piv_type);
320  } else {
321  A.ildl(L, D, perm, fill_factor, tol, pp_tol, piv_type);
322  }
323  dif = clock() - start; total += dif;
324 
325  std::string pivot_name;
326  if (piv_type == pivot_type::BKP) {
327  pivot_name = "BK";
328  } else if (piv_type == pivot_type::ROOK) {
329  pivot_name = "Rook";
330  }
331 
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() );
336  } else {
337  if (msg_lvl) printf("L is %d by %d with %d non-zeros.\n", L.n_rows(), L.n_cols(), L.nnz() );
338  }
339  if (msg_lvl) printf("\n");
340  fflush(stdout);
341 
342  // if there is a right hand side, it means the user wants a solve.
343  // TODO: refactor this solve to be in its own method, and separate
344  // factoring/minres solve phase
345  if (has_rhs) {
346  if (perform_inplace) {
347  if (msg_lvl) printf("Inplace factorization cannot be used with the solver. Please try again without -inplace.\n");
348  } else {
349  // start timer in case we're doing a full solve
350  start = clock();
351 
352  // we've permuted and equilibrated the matrix, so we gotta apply
353  // the same permutation and equilibration to the right hand side.
354  // i.e. rhs = P'S*rhs
355  // 0. apply S
356  for (int i = 0; i < A.n_cols(); i++) {
357  rhs[i] = A.S[i]*rhs[i];
358  }
359 
360  // 1. apply P' (takes rhs[perm[i]] to rhs[i], i.e. inverse of perm,
361  // where perm takes i to perm[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]];
365  }
366  rhs = tmp;
367 
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);
371  // MINRES uses the preconditioned solver that
372  // splits the block D into |D|^(1/2).
373  // For the full solver we'll just solve D directly.
374  L.backsolve(rhs, sol_vec);
375  D.solve(sol_vec, tmp);
376  L.forwardsolve(tmp, sol_vec);
377  } else {
378  start = clock();
379 
380  if (solve_type == solver_type::MINRES) {
381  // finally, since we're preconditioning with M = L|D|^(1/2), we have
382  // to multiply M^(-1) to the rhs and solve the system
383  // M^(-1) * B * M'^(-1) y = M^(-1)P'*S*b
384  L.backsolve(rhs, tmp);
385  D.sqrt_solve(tmp, rhs, false);
386 
387  if (msg_lvl) printf("Solving matrix with MINRES...\n");
388  // solve the equilibrated, preconditioned, and permuted linear system
389  minres(max_iter, minres_tol, shift);
390 
391  // now we've solved M^(-1)*B*M'^(-1)y = M^(-1)P'*S*b
392  // where B = P'SASPy.
393 
394  // but the actual solution is y = M' * P'S^(-1)*x
395  // so x = S*P*M'^(-1)*y
396 
397  // 0. apply M'^(-1)
398  D.sqrt_solve(sol_vec, tmp, true);
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);
403  }
404  }
405 
406  // 1. apply P
407  for (int i = 0; i < A.n_cols(); i++) {
408  tmp[perm[i]] = sol_vec[i];
409  }
410  sol_vec = tmp;
411 
412  // 2. apply S
413  for (int i = 0; i < A.n_cols(); i++) {
414  sol_vec[i] = A.S[i]*sol_vec[i];
415  }
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");
419 
420  if (save_sol) {
421  // save results
422  // TODO: refactor this to be in its own method
423  if (msg_lvl) printf("Solution saved to output_matrices/outsol.mtx.\n");
424  save_vector(sol_vec, "output_matrices/outsol.mtx");
425  }
426 
427  }
428  }
429  }
430 
437  void minres(int max_iter = 1000, double stop_tol = 1e-6, double shift = 0.0);
438 
444  void sqmr(int max_iter = 1000, double stop_tol = 1e-6);
445 
450  void save() { // TODO: refactor this as a "save factors" method
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);
455  } else {
456  A.save("output_matrices/outL.mtx", false);
457  }
458 
459  A.S.save("output_matrices/outS.mtx");
460  save_vector(perm, "output_matrices/outP.mtx");
461 
462  D.save("output_matrices/outD.mtx");
463  if (msg_lvl) cout << "Save complete." << endl;
464  }
465 
468  void display() {
469 #ifdef SYM_ILDL_DEBUG
470  if (perform_inplace) {
471  cout << A << endl;
472  } else {
473  cout << L << endl;
474  }
475  cout << D << endl;
476  cout << perm << endl;
477 #endif
478  }
479 };
480 
481 #include "solver_minres.h"
482 #include "solver_sqmr.h"
483 
484 }
485 
486 #endif
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&#39;, where QVQ&#39; 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
Definition: solver.h:32
void set_inplace(bool inplace)
Decides whether we perform the factorization inplace or not.
Definition: solver.h:246
Definition: solver.h:11
Definition: solver.h:23
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
Definition: solver.h:41
vector< el_type > rhs
The right hand side we&#39;ll solve for.
Definition: solver.h:146
Definition: solver.h:14
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&#39; * S * A * S * P = LDL&#39; 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