sparse.h

#ifndef __SPA__
#define __SPA__
/****************************************/
/* Sparse Matrix algorithms for CLASS   */
/* 15/11 2010                           */
/* Thomas Tram                          */
/****************************************/
#include <memory>
#include <vector>
 
#include "common.h"
 
/* Structures: */
struct sparse_matrix {
  /* Sparse matrix in compressed column form: */
  int ncols; /* Number of columns */
  int nrows; /* Number of rows */
  int maxnz; /* Maximum number of non-zero entries*/
  std::vector<int>
      Ap; /* Ap[0..ncols]. Ap[k+1]-Ap[k] is the number of entries in the k'th column. */
  std::vector<int> Ai;    /* Ai[0..(maxnz-1)]. Contains the row indices of the entries. */
  std::vector<double> Ax; /* Ax[0..(maxnz-1)]. Contains the values of the entries. */
 
  sparse_matrix(int ncols, int nrows, int maxnz);
};
 
typedef struct sparse_matrix sp_mat;
 
struct sparse_numerical {
  /* Sparse LU decomposition along with enough information to do a fast refactorization: */
  int n;                     /*Matrix assumed square, [nxn] */
  std::unique_ptr<sp_mat> L; /*L and U is the factors of the decomposed matrix.*/
  std::unique_ptr<sp_mat> U;
  std::vector<int*>
      xi; /*xi[k] points to a row of xi, which holds the topological ordered indices.*/
  std::vector<int> xi_data;
  std::vector<int> topvec; /*topvec[k] holds the first index in xi[k].*/
  std::vector<int> pinv;   /*Inverse row permutation. */
  std::vector<int> p;      /*Row permutation. */
  std::vector<int> q;      /* Column permutation */
  std::vector<int> wamd;   /* Work array for sp_amd */
  std::vector<double> w;   /* Work array for sp_lu */
 
  sparse_numerical(int n);
};
 
typedef struct sparse_numerical sp_num;
 
/* Routines and macros: */
int reachr(sp_mat* G, sp_mat* B, int k, int* xik, int* pinv);
void dfsr(int j, sp_mat* G, int* top, int* xik, int* pinv);
void sp_splsolve(sp_mat* G, sp_mat* B, int k, int* xik, int top, double* x, int* pinv);
bool sp_ludcmp(sp_num* N, sp_mat* A, double pivtol);
void sp_lusolve(sp_num* N, double* b, double* x);
void sp_refactor(sp_num* N, sp_mat* A);
int column_grouping(sp_mat* G, int* col_g, int* col_wi);
void sp_amd(int* Cp, int* Ci, int n, int cnzmax, int* P, int* W);
int sp_wclear(int mark, int lemax, int* w, int n);
int sp_tdfs(int j, int k, int* head, const int* next, int* post, int* stack);
 
#define SPFLIP(i) (-(i) - 2)
#define SPUNFLIP(i) (((i) < 0) ? SPFLIP(i) : (i))
#define SPMARKED(w, j) (w[j] < 0)
#define SPMARK(w, j)     \
  {                      \
    w[j] = SPFLIP(w[j]); \
  }
 
#endif