SparseChol.h

#ifndef SPCHOL_H
#define SPCHOL_H

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include "ExceptionBlockIP.h"

using namespace std;

// global enums so they can be reused by classes that include this one
// without any castings

enum CHOL_SOLVER {CHOLMOD, SPRSBLKLLT};
enum WHO_PERMUTES {CHOLESKY, USER_OF_CHOLESKY};
enum TYPE_MATRIX {GENERAL, NETWORK, IDENTITY, DIAGONAL, IDTY_IDTY, DIAG_DIAG, GEN_SYM_UPTR};
enum TYPE_ORIENTATION {ORIENTED, NONORIENTED};
enum NUMBERING {C_0, FORTRAN_1, NOT_COMPUTED};
 
class SparseChol {
    
  // data representation

   
 private:

  typedef struct { 
    int src; 
    int dst; 
    int arc; 
  } SRC_DST_ARC;

  CHOL_SOLVER chol_solver; 
  TYPE_MATRIX type_matrix; 
  TYPE_ORIENTATION type_orientation; 

  static constexpr double MIN_PIVOT= 0.0; 
  static constexpr double POSITIVE_PIVOT= 1.0E+128; 
  
  int m; 
  int njka; 
  int maxfillin; 
  int maxlnz; 
  int k; 
  
  int num_zero_pivots; 
  int num_semidef_matrix; 

  
  // for diagonal matrices
  double *d1; 
  double *d2; 
  
  int *pfa; 
  int *ipfa; 
  NUMBERING pfa_numbering; 
  double *permrhs; 


  int *ia; 
  int *ja; 
  int *ka; 
  int *la; 
  int nla; 
  double *va; 


  int nnu; 
  int nar; 
  int *inik; 
  int *dst2; 
  int *arck_l; 
  int *inil; 
  int *src2; 
  int *arcl_k; 

 
 
  int *xadj; 
  int *ajcncy; 
  int *workspak; 
  int *xlnz; 
  int *xlindx; 
  int *lindx; 
  int maxsub; 
  int *colcnt; 
  int *snode; 
  int maxsuper; 
  int *xsuper; 
  int iwsiz; 
  int *split;
  int tmpsiz; 
  double *lnz; 
  //

  double **plnz;

  int *ilnz; 
  int *ifillin; 
  int maxiarc; 
  int *iarc; 
  int *ini_arc; 

  
  // ********** if chol_solver is cholmod:
  // not yet implemented


  // member functions
    
public: //interface routines

    SparseChol(); // Constructor
    ~SparseChol(); // Destructor
    void reset(CHOL_SOLVER chol_solver= (CHOL_SOLVER)NULL, TYPE_MATRIX type_matrix= (TYPE_MATRIX)NULL, 
           TYPE_ORIENTATION type_orientation= ORIENTED); // frees memory and set to 0 (like calling destructor+constructor)
    void initialize(CHOL_SOLVER chol_solver= (CHOL_SOLVER)NULL, TYPE_MATRIX type_matrix= (TYPE_MATRIX)NULL, 
            TYPE_ORIENTATION type_orientation= ORIENTED);

    // symbolic_fact_MMt() is overloaded function
    // this is for GENERAL matrices
    void symbolic_fact_MMt(int m, int n, int nz, int *icola, int *inirowa, double *a, 
              int k=1, int *prov_pfa= NULL, int *prov_ipfa= NULL, 
              NUMBERING prov_pfa_numbering= NOT_COMPUTED); 
    // this is for NETWORK matrices (oriented and nonoriented)
    void symbolic_fact_MMt(int nnu, int nar, int *src, int *dst, int k=1, 
              int *prov_pfa= NULL, int *prov_ipfa= NULL, 
              NUMBERING prov_pfa_numbering_= NOT_COMPUTED);
    // this is for IDENTITY or IDTY_IDTY matrices
    void symbolic_fact_MMt(int m, int k=1);
    // this is for DIAGONAL matrices
    void symbolic_fact_MMt(int m, double *d1_in, int k=1);
    // this is for DIAG_DIAG matrices
    void symbolic_fact_MMt(int m, double *d1_in, double *d2_in, int k=1);

    void numeric_fact_MMt(double *Theta, int i_k=0);
    void numeric_solve_MMt(double *rhs, int i_k= 0, WHO_PERMUTES whoperm= (WHO_PERMUTES)NULL);

    // routines for positive semidefinite matrices M (not M*M')
    // symbolic_fact_M() is overloaded function
    // this is for GEN_SYM_UPTR matrices
    void symbolic_fact_M(int m, int nz, int *icola, int *inirowa, 
             int *prov_pfa= NULL, int *prov_ipfa= NULL, 
             NUMBERING prov_pfa_numbering= NOT_COMPUTED); 
    // this is for DIAGONAL matrices
    void symbolic_fact_M(int m);

    void numeric_fact_M(double *a);
    void numeric_solve_M(double *rhs, WHO_PERMUTES whoperm= (WHO_PERMUTES)NULL);

    // access to data in class
    int get_pfa(int i); // return pfa[i] with no check 
    int get_ipfa(int i); // return ipfa[i] with no check
    int get_maxlnz(); // return maxlnz with no check
    int get_maxfillin(); // return maxfillin with no check
    int get_njka(); // return njka with no check
    int get_num_zero_pivots(); // returns num_zero_pivots with no check
    int get_num_semidef_matrix(); // returns num_semidef_matrix with no check

private:

    void free_mem();

    void symbolic_fact_MMt_sprsblkllt_general(int n, int nz, int *icola, int *inirowa, double *a );
    void symbolic_fact_MMt_sprsblkllt_network(int *src, int *dst);

    void get_indices_a_general(int *icola, int *inirowa, double *a); 
    static bool comp_field1(const SRC_DST_ARC & i, const SRC_DST_ARC & j);
    static bool comp_field2(const SRC_DST_ARC & i, const SRC_DST_ARC & j);
    void get_ipk_ipl_network(int *src, int *dst); 
    void get_pfa_ipfa_network();
    void get_pfa_ipfa_general(int *inp_ia, int *inp_ja);
    void symbolic_AThetaAt_A(); //common to general and network matrices A*A' and symmetric matrices A
    void get_ilnz_network();
    void get_ilnz_ifillin_general(int *inp_ia, int *inp_ja);

    void numeric_fact_MMt_sprsblkllt_general(double *Theta, int i_k);
    void numeric_fact_MMt_sprsblkllt_network(double *Theta, int i_k);
    void numeric_fact_MMt_identity(double *Theta, int i_k);
    void numeric_fact_MMt_idty_idty(double *Theta, int i_k);
    void numeric_fact_MMt_diagonal(double *Theta, int i_k);
    void numeric_fact_MMt_diag_diag(double *Theta, int i_k);

    void numeric_solve_MMt_sprsblkllt(double *rhs, int i_k, WHO_PERMUTES whoperm); //common to general and network matrices
    void numeric_solve_MMt_diag(double *rhs, int i_k);

    void symbolic_fact_M_sprsblkllt_general(int nz, int *icola, int *inirowa);

    void numeric_fact_M_sprsblkllt_general(double *a);
    void numeric_fact_M_diagonal(double *d1);

    void numeric_solve_M_sprsblkllt_general(double *rhs, WHO_PERMUTES whoperm); 
    void numeric_solve_M_diagonal(double *rhs);
};

inline void SparseChol::symbolic_fact_MMt(int m, int n, int nz, int *icola, int *inirowa, 
                  double *a, int k, int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)

{
    this->m= m;
    this->k= k;
/*
  // not yet implemented: the computation of pfa and ipfa computes other 
  // needed vectors, so it is not simple

  if ((pfa != NULL) && (ipfa != NULL)) {
  for(int i=0; i<m; i++) {
  this->pfa[i]= pfa[i];
  }
  for(int i=0; i<m; i++) {
  this->ipfa[i]= ipfa[i];
  }
  }
*/
  switch(chol_solver) {
    case SPRSBLKLLT:
      symbolic_fact_MMt_sprsblkllt_general(n,nz,icola,inirowa,a);
      break;
    case CHOLMOD:
      throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
      break;
  }
}

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)

{
  this->nnu= nnu;
  this->nar= nar;
  m= nnu-1; //last node removed
  this->k= k;
/*
  // not yet implemented: the computation of pfa and ipfa computes other 
  // needed vectors, so it is not simple

  if ((pfa_ != NULL) && (ipfa_ != NULL)) {
    for(int i=0; i<m; i++) {
      pfa[i]= pfa_[i];
    }
    for(int i=0; i<m; i++) {
      ipfa[i]= ipfa_[i];
    }
  }
*/
  switch(chol_solver) {
    case SPRSBLKLLT:
      symbolic_fact_MMt_sprsblkllt_network(src,dst);
      break;
    case CHOLMOD:
      throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
      break;
  }
}

inline void SparseChol::numeric_fact_MMt(double *Theta, int i_k)

{
  switch (type_matrix) {
    // the first type of matrices involve a diagonal M*Mt matrix
    case IDENTITY:
      numeric_fact_MMt_identity(Theta,i_k);
      break;
    case IDTY_IDTY:
      numeric_fact_MMt_idty_idty(Theta,i_k);
      break;
    case DIAGONAL:
      numeric_fact_MMt_diagonal(Theta,i_k);
      break;
    case DIAG_DIAG:
      numeric_fact_MMt_diag_diag(Theta,i_k);
      break;

    // the second type of matrices need a Cholesky solver
    case GENERAL:
      switch(chol_solver) {
        case SPRSBLKLLT:
          numeric_fact_MMt_sprsblkllt_general(Theta,i_k);
          break;
        case CHOLMOD:
          throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
          break;
      }
      break;
    case NETWORK:
      switch(chol_solver) {
        case SPRSBLKLLT:
          numeric_fact_MMt_sprsblkllt_network(Theta,i_k);
          break;
        case CHOLMOD:
          throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
          break;
        }
      break;
  }
}

inline void SparseChol::numeric_solve_MMt(double *rhs, int i_k, WHO_PERMUTES whoperm)

{
  switch (type_matrix) {
    // the first type of matrices involve a diagonal M*Mt matrix
    case IDENTITY:
    case IDTY_IDTY:
    case DIAGONAL:
    case DIAG_DIAG:
      numeric_solve_MMt_diag(rhs,i_k);
      break;

    // the second type of matrices need a Cholesky solver
    case GENERAL:
    case NETWORK:
      switch(chol_solver) {
      case SPRSBLKLLT:
        numeric_solve_MMt_sprsblkllt(rhs,i_k,whoperm);
        break;
      case CHOLMOD:
        throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
        break;
      }
      break;
  }
}

inline void SparseChol::symbolic_fact_M(int m, int nz, int *icola, int *inirowa, 
                    int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)

{
    this->m= m;
/*
  // not yet implemented: the computation of pfa and ipfa computes other 
  // needed vectors, so it is not simple

  if ((pfa != NULL) && (ipfa != NULL)) {
  for(int i=0; i<m; i++) {
  this->pfa[i]= pfa[i];
  }
  for(int i=0; i<m; i++) {
  this->ipfa[i]= ipfa[i];
  }
  }
*/
  switch(chol_solver) {
    case SPRSBLKLLT:
      symbolic_fact_M_sprsblkllt_general(nz,icola,inirowa);
      break;
    case CHOLMOD:
      throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
      break;
  }
}

inline void SparseChol::numeric_fact_M(double *a)

{
  switch (type_matrix) {
  case DIAGONAL:
    numeric_fact_M_diagonal(a);
    break;
  case GEN_SYM_UPTR:
    switch(chol_solver) {
    case SPRSBLKLLT:
      numeric_fact_M_sprsblkllt_general(a);
      break;
    case CHOLMOD:
      throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
      break;
    }
  }
}

inline void SparseChol::numeric_solve_M(double *rhs, WHO_PERMUTES whoperm)

{ 
  switch (type_matrix) {
  case DIAGONAL:
    numeric_solve_M_diagonal(rhs);
    break;
  case GEN_SYM_UPTR:
    switch(chol_solver) {
    case SPRSBLKLLT:
      numeric_solve_M_sprsblkllt_general(rhs,whoperm);
      break;
    case CHOLMOD:
      throw ExceptionBlockIP(CHOLMOD_NOT_IMPLEMENTED, __FILE__, __LINE__);
      break;
    }
  }
}

inline int SparseChol::get_pfa(int i) 


{
    return(pfa[i]);
}

inline int SparseChol::get_ipfa(int i)

{
    return(ipfa[i]);
}

inline int SparseChol::get_maxlnz()

{
    return (maxlnz);
}

inline int SparseChol::get_maxfillin()

{
    return (maxfillin);
}

inline int SparseChol::get_njka()

{
    return (njka);
}

inline int SparseChol::get_num_zero_pivots()

{
    return (num_zero_pivots);
}

inline int SparseChol::get_num_semidef_matrix()

{
    return (num_semidef_matrix);
}

#endif //SPCHOL_H