BlockIP
SparseChol.h
Go to the documentation of this file.
1 #ifndef SPCHOL_H
2 #define SPCHOL_H
3 
4 #include <iostream>
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include "ExceptionBlockIP.h"
8 
9 using namespace std;
10 
11 // global enums so they can be reused by classes that include this one
12 // without any castings
13 
24 
26 class SparseChol {
27 
28  // data representation
29 
30 
31  private:
32 
34  typedef struct {
35  int src;
36  int dst;
37  int arc;
38  } SRC_DST_ARC;
39 
43 
45  static constexpr double MIN_PIVOT= 0.0;
47  static constexpr double POSITIVE_PIVOT= 1.0E+128;
48 
49  int m;
50  int njka;
51  int maxfillin;
52  int maxlnz;
53  int k;
54 
57 
58 
59  // for diagonal matrices
60  double *d1;
61  double *d2;
62 
76  int *pfa;
77  int *ipfa;
79  double *permrhs;
80 
81 
105  int *ia;
106  int *ja;
107  int *ka;
108  int *la;
109  int nla;
110  double *va;
111 
112 
123  int nnu;
124  int nar;
125  int *inik;
126  int *dst2;
127  int *arck_l;
128  int *inil;
129  int *src2;
130  int *arcl_k;
131 
132 
133 
153  int *xadj;
154  int *ajcncy;
155  int *workspak;
156  int *xlnz;
157  int *xlindx;
158  int *lindx;
159  int maxsub;
160  int *colcnt;
161  int *snode;
162  int maxsuper;
163  int *xsuper;
164  int iwsiz;
165  int *split;
166  int tmpsiz;
167  double *lnz;
168  //
170 
174  double **plnz;
176 
191  int *ilnz;
192  int *ifillin;
193  int maxiarc;
194  int *iarc;
195  int *ini_arc;
196 
197 
198  // ********** if chol_solver is cholmod:
199  // not yet implemented
200 
201 
202  // member functions
203 
204 public: //interface routines
205 
206  SparseChol(); // Constructor
207  ~SparseChol(); // Destructor
208  void reset(CHOL_SOLVER chol_solver= (CHOL_SOLVER)NULL, TYPE_MATRIX type_matrix= (TYPE_MATRIX)NULL,
209  TYPE_ORIENTATION type_orientation= ORIENTED); // frees memory and set to 0 (like calling destructor+constructor)
210  void initialize(CHOL_SOLVER chol_solver= (CHOL_SOLVER)NULL, TYPE_MATRIX type_matrix= (TYPE_MATRIX)NULL,
211  TYPE_ORIENTATION type_orientation= ORIENTED);
212 
213  // symbolic_fact_MMt() is overloaded function
214  // this is for GENERAL matrices
215  void symbolic_fact_MMt(int m, int n, int nz, int *icola, int *inirowa, double *a,
216  int k=1, int *prov_pfa= NULL, int *prov_ipfa= NULL,
217  NUMBERING prov_pfa_numbering= NOT_COMPUTED);
218  // this is for NETWORK matrices (oriented and nonoriented)
219  void symbolic_fact_MMt(int nnu, int nar, int *src, int *dst, int k=1,
220  int *prov_pfa= NULL, int *prov_ipfa= NULL,
221  NUMBERING prov_pfa_numbering_= NOT_COMPUTED);
222  // this is for IDENTITY or IDTY_IDTY matrices
223  void symbolic_fact_MMt(int m, int k=1);
224  // this is for DIAGONAL matrices
225  void symbolic_fact_MMt(int m, double *d1_in, int k=1);
226  // this is for DIAG_DIAG matrices
227  void symbolic_fact_MMt(int m, double *d1_in, double *d2_in, int k=1);
228 
229  void numeric_fact_MMt(double *Theta, int i_k=0);
230  void numeric_solve_MMt(double *rhs, int i_k= 0, WHO_PERMUTES whoperm= (WHO_PERMUTES)NULL);
231 
232  // routines for positive semidefinite matrices M (not M*M')
233  // symbolic_fact_M() is overloaded function
234  // this is for GEN_SYM_UPTR matrices
235  void symbolic_fact_M(int m, int nz, int *icola, int *inirowa,
236  int *prov_pfa= NULL, int *prov_ipfa= NULL,
237  NUMBERING prov_pfa_numbering= NOT_COMPUTED);
238  // this is for DIAGONAL matrices
239  void symbolic_fact_M(int m);
240 
241  void numeric_fact_M(double *a);
242  void numeric_solve_M(double *rhs, WHO_PERMUTES whoperm= (WHO_PERMUTES)NULL);
243 
244  // access to data in class
245  int get_pfa(int i); // return pfa[i] with no check
246  int get_ipfa(int i); // return ipfa[i] with no check
247  int get_maxlnz(); // return maxlnz with no check
248  int get_maxfillin(); // return maxfillin with no check
249  int get_njka(); // return njka with no check
250  int get_num_zero_pivots(); // returns num_zero_pivots with no check
251  int get_num_semidef_matrix(); // returns num_semidef_matrix with no check
252 
253 private:
254 
255  void free_mem();
256 
257  void symbolic_fact_MMt_sprsblkllt_general(int n, int nz, int *icola, int *inirowa, double *a );
258  void symbolic_fact_MMt_sprsblkllt_network(int *src, int *dst);
259 
260  void get_indices_a_general(int *icola, int *inirowa, double *a);
261  static bool comp_field1(const SRC_DST_ARC & i, const SRC_DST_ARC & j);
262  static bool comp_field2(const SRC_DST_ARC & i, const SRC_DST_ARC & j);
263  void get_ipk_ipl_network(int *src, int *dst);
264  void get_pfa_ipfa_network();
265  void get_pfa_ipfa_general(int *inp_ia, int *inp_ja);
266  void symbolic_AThetaAt_A(); //common to general and network matrices A*A' and symmetric matrices A
267  void get_ilnz_network();
268  void get_ilnz_ifillin_general(int *inp_ia, int *inp_ja);
269 
270  void numeric_fact_MMt_sprsblkllt_general(double *Theta, int i_k);
271  void numeric_fact_MMt_sprsblkllt_network(double *Theta, int i_k);
272  void numeric_fact_MMt_identity(double *Theta, int i_k);
273  void numeric_fact_MMt_idty_idty(double *Theta, int i_k);
274  void numeric_fact_MMt_diagonal(double *Theta, int i_k);
275  void numeric_fact_MMt_diag_diag(double *Theta, int i_k);
276 
277  void numeric_solve_MMt_sprsblkllt(double *rhs, int i_k, WHO_PERMUTES whoperm); //common to general and network matrices
278  void numeric_solve_MMt_diag(double *rhs, int i_k);
279 
280  void symbolic_fact_M_sprsblkllt_general(int nz, int *icola, int *inirowa);
281 
282  void numeric_fact_M_sprsblkllt_general(double *a);
283  void numeric_fact_M_diagonal(double *d1);
284 
285  void numeric_solve_M_sprsblkllt_general(double *rhs, WHO_PERMUTES whoperm);
286  void numeric_solve_M_diagonal(double *rhs);
287 };
288 
289 inline void SparseChol::symbolic_fact_MMt(int m, int n, int nz, int *icola, int *inirowa,
290  double *a, int k, int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)
292 
309 {
310  this->m= m;
311  this->k= k;
312 /*
313  // not yet implemented: the computation of pfa and ipfa computes other
314  // needed vectors, so it is not simple
315 
316  if ((pfa != NULL) && (ipfa != NULL)) {
317  for(int i=0; i<m; i++) {
318  this->pfa[i]= pfa[i];
319  }
320  for(int i=0; i<m; i++) {
321  this->ipfa[i]= ipfa[i];
322  }
323  }
324 */
325  switch(chol_solver) {
326  case SPRSBLKLLT:
327  symbolic_fact_MMt_sprsblkllt_general(n,nz,icola,inirowa,a);
328  break;
329  case CHOLMOD:
330  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
331  break;
332  }
333 }
334 
335 inline void SparseChol::symbolic_fact_MMt(int nnu, int nar, int *src, int *dst, int k, int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)
337 
353 {
354  this->nnu= nnu;
355  this->nar= nar;
356  m= nnu-1; //last node removed
357  this->k= k;
358 /*
359  // not yet implemented: the computation of pfa and ipfa computes other
360  // needed vectors, so it is not simple
361 
362  if ((pfa_ != NULL) && (ipfa_ != NULL)) {
363  for(int i=0; i<m; i++) {
364  pfa[i]= pfa_[i];
365  }
366  for(int i=0; i<m; i++) {
367  ipfa[i]= ipfa_[i];
368  }
369  }
370 */
371  switch(chol_solver) {
372  case SPRSBLKLLT:
373  symbolic_fact_MMt_sprsblkllt_network(src,dst);
374  break;
375  case CHOLMOD:
376  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
377  break;
378  }
379 }
380 
381 inline void SparseChol::numeric_fact_MMt(double *Theta, int i_k)
383 
388 {
389  switch (type_matrix) {
390  // the first type of matrices involve a diagonal M*Mt matrix
391  case IDENTITY:
392  numeric_fact_MMt_identity(Theta,i_k);
393  break;
394  case IDTY_IDTY:
395  numeric_fact_MMt_idty_idty(Theta,i_k);
396  break;
397  case DIAGONAL:
398  numeric_fact_MMt_diagonal(Theta,i_k);
399  break;
400  case DIAG_DIAG:
401  numeric_fact_MMt_diag_diag(Theta,i_k);
402  break;
403 
404  // the second type of matrices need a Cholesky solver
405  case GENERAL:
406  switch(chol_solver) {
407  case SPRSBLKLLT:
408  numeric_fact_MMt_sprsblkllt_general(Theta,i_k);
409  break;
410  case CHOLMOD:
411  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
412  break;
413  }
414  break;
415  case NETWORK:
416  switch(chol_solver) {
417  case SPRSBLKLLT:
418  numeric_fact_MMt_sprsblkllt_network(Theta,i_k);
419  break;
420  case CHOLMOD:
421  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
422  break;
423  }
424  break;
425  }
426 }
427 
428 inline void SparseChol::numeric_solve_MMt(double *rhs, int i_k, WHO_PERMUTES whoperm)
430 
442 {
443  switch (type_matrix) {
444  // the first type of matrices involve a diagonal M*Mt matrix
445  case IDENTITY:
446  case IDTY_IDTY:
447  case DIAGONAL:
448  case DIAG_DIAG:
449  numeric_solve_MMt_diag(rhs,i_k);
450  break;
451 
452  // the second type of matrices need a Cholesky solver
453  case GENERAL:
454  case NETWORK:
455  switch(chol_solver) {
456  case SPRSBLKLLT:
457  numeric_solve_MMt_sprsblkllt(rhs,i_k,whoperm);
458  break;
459  case CHOLMOD:
460  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
461  break;
462  }
463  break;
464  }
465 }
466 
467 inline void SparseChol::symbolic_fact_M(int m, int nz, int *icola, int *inirowa,
468  int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)
470 
484 {
485  this->m= m;
486 /*
487  // not yet implemented: the computation of pfa and ipfa computes other
488  // needed vectors, so it is not simple
489 
490  if ((pfa != NULL) && (ipfa != NULL)) {
491  for(int i=0; i<m; i++) {
492  this->pfa[i]= pfa[i];
493  }
494  for(int i=0; i<m; i++) {
495  this->ipfa[i]= ipfa[i];
496  }
497  }
498 */
499  switch(chol_solver) {
500  case SPRSBLKLLT:
501  symbolic_fact_M_sprsblkllt_general(nz,icola,inirowa);
502  break;
503  case CHOLMOD:
504  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
505  break;
506  }
507 }
508 
509 inline void SparseChol::numeric_fact_M(double *a)
511 
516 {
517  switch (type_matrix) {
518  case DIAGONAL:
519  numeric_fact_M_diagonal(a);
520  break;
521  case GEN_SYM_UPTR:
522  switch(chol_solver) {
523  case SPRSBLKLLT:
524  numeric_fact_M_sprsblkllt_general(a);
525  break;
526  case CHOLMOD:
527  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
528  break;
529  }
530  }
531 }
532 
533 inline void SparseChol::numeric_solve_M(double *rhs, WHO_PERMUTES whoperm)
535 
545 {
546  switch (type_matrix) {
547  case DIAGONAL:
548  numeric_solve_M_diagonal(rhs);
549  break;
550  case GEN_SYM_UPTR:
551  switch(chol_solver) {
552  case SPRSBLKLLT:
553  numeric_solve_M_sprsblkllt_general(rhs,whoperm);
554  break;
555  case CHOLMOD:
556  throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
557  break;
558  }
559  }
560 }
561 
562 inline int SparseChol::get_pfa(int i)
563 
564 
568 {
569  return(pfa[i]);
570 }
571 
572 inline int SparseChol::get_ipfa(int i)
574 
578 {
579  return(ipfa[i]);
580 }
581 
584 
588 {
589  return (maxlnz);
590 }
591 
594 
598 {
599  return (maxfillin);
600 }
601 
604 
608 {
609  return (njka);
610 }
611 
614 
618 {
619  return (num_zero_pivots);
620 }
621 
624 
628 {
629  return (num_semidef_matrix);
630 }
631 
632 #endif //SPCHOL_H