BlockIP
Classes | Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
SparseChol Class Reference

Class for sparse Cholesky factorizations. More...

#include <SparseChol.h>

Classes

struct  SRC_DST_ARC
 Auxiliary struct for sorting network structure. More...
 

Public Member Functions

 SparseChol ()
 Constructor. It calls initialize().
 
 ~SparseChol ()
 Destructor. It calls free_mem().
 
void reset (CHOL_SOLVER chol_solver=(CHOL_SOLVER) NULL, TYPE_MATRIX type_matrix=(TYPE_MATRIX) NULL, TYPE_ORIENTATION type_orientation=ORIENTED)
 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)
 Initialize attributes (set variables to 0, NULL...)
 
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)
 Computes symmetric ordering and symbolic factorization of AA' for a general matrix A.
 
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)
 Computes symmetric ordering and symbolic factorization of AA' for a network matrix A.
 
void symbolic_fact_MMt (int m, int k=1)
 Computes symbolic factorization of AA' for an IDENTITY or IDTY_IDTY matrix A.
 
void symbolic_fact_MMt (int m, double *d1_in, int k=1)
 Computes symbolic factorization of AA' for a DIAG matrix A.
 
void symbolic_fact_MMt (int m, double *d1_in, double *d2_in, int k=1)
 Computes symbolic factorization of AA' for a DIAG_DIAG matrix A.
 
void numeric_fact_MMt (double *Theta, int i_k=0)
 Computes numerical factorization of A*Theta*A'.
 
void numeric_solve_MMt (double *rhs, int i_k=0, WHO_PERMUTES whoperm=(WHO_PERMUTES) NULL)
 Solve systems with matrix A*Theta*A' and right-hand-side rhs.
 
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)
 Computes symmetric ordering and symbolic factorization of A for a general symmetric matrix.
 
void symbolic_fact_M (int m)
 Computes symbolic factorization of DIAGONAL matrix A.
 
void numeric_fact_M (double *a)
 Computes numerical factorization for a positive semidefinite matrix A.
 
void numeric_solve_M (double *rhs, WHO_PERMUTES whoperm=(WHO_PERMUTES) NULL)
 Solve systems with a positive semidefinite symmetric matrix A and right-hand-side rhs.
 
int get_pfa (int i)
 Provides pfa[i] with no check.
 
int get_ipfa (int i)
 Provides ipfa[i] with no check.
 
int get_maxlnz ()
 Provides maxlnz with no check.
 
int get_maxfillin ()
 Provides maxfillin with no check.
 
int get_njka ()
 Provides njka with no check.
 
int get_num_zero_pivots ()
 Provides num_zero_pivots with no check.
 
int get_num_semidef_matrix ()
 Provides num_semidef_matrix with no check.
 

Private Member Functions

void free_mem ()
 Frees memory (uses C free(), not C++ delete[])
 
void symbolic_fact_MMt_sprsblkllt_general (int n, int nz, int *icola, int *inirowa, double *a)
 Computes symmetric ordering and symbolic factorization for AA' using Ng-Peyton package for a general matrix.
 
void symbolic_fact_MMt_sprsblkllt_network (int *src, int *dst)
 Computes symmetric ordering and symbolic factorization for AA' using Ng-Peyton package for a network matrix.
 
void get_indices_a_general (int *icola, int *inirowa, double *a)
 Computes data structure of arrays ia, ja, ka, la, va as detailed in Monma and Morton paper.
 
void get_ipk_ipl_network (int *src, int *dst)
 Computes data structures with network sorted by src and dst, to be used later.
 
void get_pfa_ipfa_network ()
 Computes row permutation of AA' for a network matrix A.
 
void get_pfa_ipfa_general (int *inp_ia, int *inp_ja)
 Computes row permutation of AA'/A for a general matrix.
 
void symbolic_AThetaAt_A ()
 Symbolic factorization of AThetaA' or symmetric A.
 
void get_ilnz_network ()
 Computes variables and indices for later fast filling of lnz for a network matrix.
 
void get_ilnz_ifillin_general (int *inp_ia, int *inp_ja)
 Compute arrays of indices to lnz(): array ilnz(njka), for a general matrix.
 
void numeric_fact_MMt_sprsblkllt_general (double *Theta, int i_k)
 Computes numerical factorization of A*Theta*A' when A is general matrix.
 
void numeric_fact_MMt_sprsblkllt_network (double *Theta, int i_k)
 Computes numerical factorization of A*Theta*A' when A is a network.
 
void numeric_fact_MMt_identity (double *Theta, int i_k)
 Numerical factorization of I*Theta*I'.
 
void numeric_fact_MMt_idty_idty (double *Theta, int i_k)
 Numerical factorization of [I I]*(Theta^+, Theta^-)*[I I]' matrix.
 
void numeric_fact_MMt_diagonal (double *Theta, int i_k)
 Numerical factorization of D*Theta*D' (diagonal matrix)
 
void numeric_fact_MMt_diag_diag (double *Theta, int i_k)
 Numerical factorization of [D1 D2]*(Theta^+, Theta^-)*[D1 D2]' (DIAG_DIAG matrix)
 
void numeric_solve_MMt_sprsblkllt (double *rhs, int i_k, WHO_PERMUTES whoperm)
 Numerical solve for GENERAL and NETWORK matrices (require sprsblkllt)
 
void numeric_solve_MMt_diag (double *rhs, int i_k)
 Numerical solve for IDENTITY, IDTY_IDTY, DIAGONAL, DIAG_DIAG (M*M' is diagonal)
 
void symbolic_fact_M_sprsblkllt_general (int nz, int *icola, int *inirowa)
 Computes symmetric ordering and symbolic factorization of symmetric upper triangular A using Ng-Peyton package.
 
void numeric_fact_M_sprsblkllt_general (double *a)
 Computes numerical factorization of A of symmetric upper triangular A using Ng-Peyton package.
 
void numeric_fact_M_diagonal (double *d1)
 Numerical factorization of D diagonal positive semidefinite matrix.
 
void numeric_solve_M_sprsblkllt_general (double *rhs, WHO_PERMUTES whoperm)
 Numerical solve for symmetric matrices.
 
void numeric_solve_M_diagonal (double *rhs)
 Numerical solve for DIAGONAL matrix.
 

Static Private Member Functions

static bool comp_field1 (const SRC_DST_ARC &i, const SRC_DST_ARC &j)
 Auxiliary routine to compare first field of type SRC_DST_ARC, used in sorting arcs.
 
static bool comp_field2 (const SRC_DST_ARC &i, const SRC_DST_ARC &j)
 Auxiliary routine to compare second field of type SRC_DST_ARC, used in sorting arcs.
 

Private Attributes

CHOL_SOLVER chol_solver
 Cholesky solver to be used.
 
TYPE_MATRIX type_matrix
 Type of matrix (general, network, identity, etc),.
 
TYPE_ORIENTATION type_orientation
 Type of arc orientation for network matrices (by default, oriented).
 
int m
 Dimension of system A*A'/A : m x m matrix.
 
int njka
 Nonzeros in diagonal and subdiagonal of A*A'/A.
 
int maxfillin
 Nonzeros in factorization of A*A'/A due to fill-in.
 
int maxlnz
 Total nonzeros in factorization of A*A'/A (= maxfillin+njka).
 
int k
 Number of matrices A*A' with same topology to be factorized, see below explanation of lnz() and pnlz() (k=1 for A)
 
int num_zero_pivots
 Number of zero pivots found during factorizations.
 
int num_semidef_matrix
 Number of semidefinite matrices found during factorizations.
 
double * d1
 Diagonal of DIAG or first diagonal of DIAG_DIAG [D1 D2] matrices.
 
double * d2
 Second diagonal of DIAG_DIAG [D1 D2] matrices.
 

Vectors related to symmetric row-column permutation:

int * pfa
 pfa[i]= k : original row k is now (after permutation) in position i
 
int * ipfa
 ipfa[k]= i : i is current row (after permutation) of original row k.
 
NUMBERING pfa_numbering
 Numbering style of pfa (used in C and Fortran code).
 
double * permrhs
 Auxliary array to store the permutation of the rhs of the system to be solved.
 

When chol_solver is sprsblkllt and A is 'GENERAL' matrix :

int * ia
 ia(m+1): pointer to first in ja().
 
int * ja
 ja(njka): pair of rows (ia,ja) of A that generate arc in graph of A*A' .
 
int * ka
 ka(njka+1): pointer to first in la() and va().
 
int * la
 la(nla): columns that intervene in arc (i,j) in graph of A*A'.
 
int nla
 Size of la() and va().
 
double * va
 va(nla): premultiplied values associated to column of la().
 

When chol_solver is sprsblkllt and A is 'NETWORK' matrix:

int nnu
 Number of nodes.
 
int nar
 Number of arcs.
 
int * inik
 inik(nnu+1): pointers to first arc sorted by source (compressed form).
 
int * dst2
 dst2(nar): destination nodes of arcs in inik.
 
int * arck_l
 arck_l(nar): arcs of inik, dst2.
 
int * inil
 inil(nnu+1): pointers to first arc sorted by destination (compressed form).
 
int * src2
 src2(nar): source nodes of arcs in inik
 
int * arcl_k
 arcl_k(nar): arcs of inil, src2.
 

Variables needed by sprsblkllt:

int * xadj
 xadj(m+1): indices to adjacency matrix ajcncy().
 
int * ajcncy
 ajcncy(2*(njka-m)): 2*arcs in graph of A*A'/A without diagonal entries (i,i) of A*A'/A.
 
int * workspak
 
int * xlnz
 xlnz(m+1): indices to first column/row in lnz(.).
 
int * xlindx
 xlindx(maxsuper): indices to lindx (compressed form).
 
int * lindx
 lindx(maxsub): compressed form of sub and diagonal of factorization.
 
int maxsub
 Maximum size of some vectors, computed by sprsblkllt.
 
int * colcnt
 colnct(m): elements per column in factorization.
 
int * snode
 snode(m): supernode of each column.
 
int maxsuper
 Maximum size of some vectors, computed by sprsblkllt.
 
int * xsuper
 xsuper(m+1): pointer to first column of supernode.
 
int iwsiz
 Size of integer working space, required by sprsblkllt.
 
int * split
 split(m): splitting of supernodes for exploiting cache memory.
 
int tmpsiz
 Size of temporary working space, required by sprsblkllt.
 
double * lnz
 lnz(k*maxlnz): sub(upper) and diagonal elements of ATheta(k)A' (and its Cholesky factorization, once performed).
 
double ** plnz
 plnz(k): pointers to lnz(i).
 

Variables with indices for fast filling of matrix to factorize:

int * ilnz
 ilnz(njka): indices to lnz(.) of each nonzero element of A*Theta*A'.
 
int * ifillin
 ifillin(maxfillin): fill-in positions in lnz(.) to be cleaned before filling lnz() with A*Theta*A'.
 
int maxiarc
 Size of iarc.
 
int * iarc
 iarc(maxiarc): indices to arcs intervening in each element nonzero of lnz().
 
int * ini_arc
 ini_arc(njka+1): begin of arcs in iarc for each element of lnz().
 

Static Private Attributes

static constexpr double MIN_PIVOT = 0.0
 Threshold pivot value for Cholesky factorizations of "simple" (diag, idty...) matrices (no pivot less than or equal to MIN_PIVOT allowed)
 
static constexpr double POSITIVE_PIVOT = 1.0E+128
 New value for pivots below MIN_PIVOT in Cholesky factorizations of "simple" matrices (diag, idty...)
 

Detailed Description

Class for sparse Cholesky factorizations.

Constructor & Destructor Documentation

SparseChol::SparseChol ( )

Constructor. It calls initialize().

SparseChol::~SparseChol ( )

Destructor. It calls free_mem().

Member Function Documentation

bool SparseChol::comp_field1 ( const SRC_DST_ARC i,
const SRC_DST_ARC j 
)
staticprivate

Auxiliary routine to compare first field of type SRC_DST_ARC, used in sorting arcs.

bool SparseChol::comp_field2 ( const SRC_DST_ARC i,
const SRC_DST_ARC j 
)
staticprivate

Auxiliary routine to compare second field of type SRC_DST_ARC, used in sorting arcs.

void SparseChol::free_mem ( )
private

Frees memory (uses C free(), not C++ delete[])

void SparseChol::get_ilnz_ifillin_general ( int *  inp_ia,
int *  inp_ja 
)
private

Compute arrays of indices to lnz(): array ilnz(njka), for a general matrix.

For each pair (i,j) stated by inp_ia and inp_ja of diagonal and subdiagonal (or supradiagonal, matrix is symmetric) elements of A*Theta*A'/A ilnz(.) points to position in lnz(.) Array ifillin() with indices to fill-in positions in lnz() is also computed (to fastly clean those positions when needed).

Parameters
inpia Vector of beginning of rows of A*A'/A(compressed form)
inp_jaVector of column indices A*A'/A(.
void SparseChol::get_ilnz_network ( )
private

Computes variables and indices for later fast filling of lnz for a network matrix.

Computes variables and indices for later fast filling of lnz(): maxiarc, iarc(maxiarc), ini_arc(njka+1)

void SparseChol::get_indices_a_general ( int *  icola,
int *  inirowa,
double *  a 
)
private

Computes data structure of arrays ia, ja, ka, la, va as detailed in Monma and Morton paper.

It deals with matrices with empy (zero) rows by adding a fictitious diagonal zero entry. The zero diagonal pivot will be later modified by the sparse Cholesky solver to deal with this positive semidefinite matrix.

Parameters
icolaVector of column indices.
inirowaVector of beginning of rows (compressed form)
aVector of matrix elements.
int SparseChol::get_ipfa ( int  i)
inline

Provides ipfa[i] with no check.

Returns
ipfa[i]
Note
The user must guarantee ipfa previously computed in call to symbolic_fact_MMt().
void SparseChol::get_ipk_ipl_network ( int *  src,
int *  dst 
)
private

Computes data structures with network sorted by src and dst, to be used later.

Parameters
srcVector of arc sources.
dstVector of arc destinations.
int SparseChol::get_maxfillin ( )
inline

Provides maxfillin with no check.

return maxfillin

Note
The user must guarantee maxfillin previously computed in call to symbolic_fact_MMt().
int SparseChol::get_maxlnz ( )
inline

Provides maxlnz with no check.

Returns
maxlnz
Note
The user must guarantee maxlnz previously computed in call to symbolic_fact_MMt().
int SparseChol::get_njka ( )
inline

Provides njka with no check.

Returns
njka
Note
The user must guarantee njka previously computed in call to symbolic_fact_MMt().
int SparseChol::get_num_semidef_matrix ( )
inline

Provides num_semidef_matrix with no check.

Returns
num_semidef_matrix
Note
It will 0 if no numeric factorization made.
int SparseChol::get_num_zero_pivots ( )
inline

Provides num_zero_pivots with no check.

Returns
num_zero_pivots
Note
It will 0 if no numeric factorization made.
int SparseChol::get_pfa ( int  i)
inline

Provides pfa[i] with no check.

Returns
pfa[i]
Note
The user must guarantee pfa previously computed in call to symbolic_fact_MMt().
void SparseChol::get_pfa_ipfa_general ( int *  inp_ia,
int *  inp_ja 
)
private

Computes row permutation of AA'/A for a general matrix.

Computes row permutation of AA'/A using minimum degree ordering to the graph associated to AA'/A (AA'/A is symmetric matrix). It calls routines orddmd() of sprskblkllt to compute pfa (direct permutation) and ipfa(inverse permutation): pfa[i]= k : original row k is now (after permutation) in position i ipfa[k]= i : i is current row (after permutation) of original row k Both pfa and ipfa needed because permutation may be nonsymmetric

WARNING: ordmmd() is a fortran routine, thus permutation in pfa/ipfa is numbered from 1, not from 0 as in C. 1-numbering needed for sprsblkllt routines. After calling sfnit() and symfct() then they can be numbered in C 0-starting style

Parameters
inpia Vector of beginning of rows of A*A'/A(compressed form)
inp_jaVector of column indices A*A'/A(.
void SparseChol::get_pfa_ipfa_network ( )
private

Computes row permutation of AA' for a network matrix A.

Computes row permutation of AA' using minimum degree ordering to the graph associated to AA' (A is network matrix without last node). It calls routines orddmd() of sprskblkllt to compute pfa (direct permutation) and ipfa(inverse permutation):

pfa[i]= k : original row k is now (after permutation) in position i
ipfa[k]= i : i is current row (after permutation) of original row k

Both pfa and ipfa needed because permutation may be nonsymmetric

Note
ordmmd() is a fortran routine, thus permutation in pfa/ipfa is numbered from 1, not from 0 as in C. 1-numbering needed for sprsblkllt routines. After calling sfnit() and symfct() then they can be numbered in C 0-starting style.
void SparseChol::initialize ( CHOL_SOLVER  chol_solver = (CHOL_SOLVER)NULL,
TYPE_MATRIX  type_matrix = (TYPE_MATRIX)NULL,
TYPE_ORIENTATION  type_orientation = ORIENTED 
)

Initialize attributes (set variables to 0, NULL...)

Parameters
chol_solverCholesky solver to be used
type_matrixType of matrix (general, network, IDTY, ...)
type_orientation(Only for network matrices) Whether arcs are oriented or nonoriented.
void SparseChol::numeric_fact_M ( double *  a)
inline

Computes numerical factorization for a positive semidefinite matrix A.

A is either diagonal or general symmetric upper triangular (with diagonal) matrix in compressed rowwise order

Parameters
a(nz)Matrix elements (only diagonal and upper triangle, or diagonal)
Note
symbolic_fact_M() must have been previously called.
void SparseChol::numeric_fact_M_diagonal ( double *  d1)
private

Numerical factorization of D diagonal positive semidefinite matrix.

When matrix is DIAG just allocate some vectors and fill some minimum information.

Parameters
d1diagonal terms of D, d1 of dimension m
void SparseChol::numeric_fact_M_sprsblkllt_general ( double *  a)
private

Computes numerical factorization of A of symmetric upper triangular A using Ng-Peyton package.

A is general symmetric (diagonal and upper triangle) matrix in compressed rowwise order We assume the topology (icola,inirowa) of matrix A is the same that was used in the symbolic factorization, otherwise the code may fail.

Parameters
a(nzor njka) Matrix elements
void SparseChol::numeric_fact_MMt ( double *  Theta,
int  i_k = 0 
)
inline

Computes numerical factorization of A*Theta*A'.

Parameters
Theta(n)Diagonal matrix Theta of A*Theta*A', of size number of columns of A.
i_kWhich matrix to factorize, from 0 to k-1.
Note
symbolic_fact_MMt() must have been previously called.
void SparseChol::numeric_fact_MMt_diag_diag ( double *  Theta,
int  i_k 
)
private

Numerical factorization of [D1 D2]*(Theta^+, Theta^-)*[D1 D2]' (DIAG_DIAG matrix)

Numerical factorization of [D1 D2]*(Theta^+, Theta^-)*[D1 D2]': for DIAG_DIAG just stores D1*(Theta^+)*D1 + D2*(Theta^-)*D2 in lnz.

Theta(n=2*m), where Theta(0:m-1) is Theta^+ and Theta(m:n-1) is Theta^-

int i_k: which matrix to factorize, from 0 to k-1

void SparseChol::numeric_fact_MMt_diagonal ( double *  Theta,
int  i_k 
)
private

Numerical factorization of D*Theta*D' (diagonal matrix)

Numerical factorization of D*Theta*D': when D is diagonal just stores D*Theta*D in lnz.

Theta(n): diagonal matrix of size of D int i_k: which matrix to factorize, from 0 to k-1

void SparseChol::numeric_fact_MMt_identity ( double *  Theta,
int  i_k 
)
private

Numerical factorization of I*Theta*I'.

Numerical factorization of I*Theta*I': when I is identity just stores Theta in lnz.

Theta(n): diagonal matrix of size dimension of I int i_k: which matrix to factorize, from 0 to k-1

void SparseChol::numeric_fact_MMt_idty_idty ( double *  Theta,
int  i_k 
)
private

Numerical factorization of [I I]*(Theta^+, Theta^-)*[I I]' matrix.

Numerical factorization of [I I]*(Theta^+, Theta^-)*[I I]': for IDTY_IDTY just stores (Theta^+) + (Theta^-) in lnz.

Theta(n=2*m), where Theta(0:m-1) is Theta^+ and Theta(m:n-1) is Theta^-

int i_k: which matrix to factorize, from 0 to k-1

void SparseChol::numeric_fact_MMt_sprsblkllt_general ( double *  Theta,
int  i_k 
)
private

Computes numerical factorization of A*Theta*A' when A is general matrix.

Parameters
Theta(n)Diagonal matrix Theta of A*Theta*A', of size number of columns of A
i_kWhich matrix to factorize, from 0 to k-1
void SparseChol::numeric_fact_MMt_sprsblkllt_network ( double *  Theta,
int  i_k 
)
private

Computes numerical factorization of A*Theta*A' when A is a network.

Computes numerical factorization of A*Theta*A' when A is either an ORIENTED or NONORIENTED network matrix. Exploits that get_ilnz_network computed indices for fast filling of lnz for the ORIENTED case.

If ORIENTED:

A=N of size m x nar (nar=n), and A*Theta*A'= N*Theta*N', where

Theta(n): diagonal matrix Theta of A*Theta*A', of size number of columns of A

If NONORIENTED:

A=[N -N] of size m x (2*nar) (n=2*nar) and A*Theta*A'= N*(Theta^+ + Theta^-)*N' where

Theta(n=2*nar), where Theta(0:nar-1) is Theta^+ and Theta(nar:n-1) is Theta^-

Parameters
ThetaDiagonal Theta matrix.
i_kWhich matrix to factorize, from 0 to k-1.
void SparseChol::numeric_solve_M ( double *  rhs,
WHO_PERMUTES  whoperm = (WHO_PERMUTES)NULL 
)
inline

Solve systems with a positive semidefinite symmetric matrix A and right-hand-side rhs.

The solution is returned by the same vector rhs (then it is overwritten).

Parameters
rhsOn input vector rhs; on output the solution vector.
whopermIf CHOLESKY, then the reordering by pfa/ipfa must internally be applied to rhs by the routine; if USER_OF_CHOLESKY, then the user already provides rhs in the permutted order given by pfa/ipfa.
Note
symbolic_fact_M() and numeric_fact_M() must have been previously called.
void SparseChol::numeric_solve_M_diagonal ( double *  rhs)
private

Numerical solve for DIAGONAL matrix.

Solves M*sol = rhs, when M is a diagonal positive semidefinite matrix.

Parameters
rhson input rhs vector, on output it contains the solution.

Previous calls to symbolic_fact_M and num_fact_M required for a correct "factorization" of M (i.e., storage of in lnz of the diagonal terms)

void SparseChol::numeric_solve_M_sprsblkllt_general ( double *  rhs,
WHO_PERMUTES  whoperm = (WHO_PERMUTES)NULL 
)
private

Numerical solve for symmetric matrices.

Solves A*sol = rhs. On output, rhs contains the solution. Previous calls to symbolic_fact_M and numeric_fact_M required for a correct factorization of A

whoperm: if CHOLESKY, then the reordering by pfa/ipfa must internally be applied to rhs by the routine; if USER_OF_CHOLESKY, then the user already provides rhs in the permutted order given by pfa/ipfa.

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

Solve systems with matrix A*Theta*A' and right-hand-side rhs.

The solution is returned by the same vector rhs (then it is overwritten).

Parameters
rhsOn input vector rhs; on output the solution vector.
whopermIf CHOLESKY, then the reordering by pfa/ipfa must internally be applied to rhs by the routine; if USER_OF_CHOLESKY, then the user already provides rhs in the permutted order given by pfa/ipfa.
i_kWhich matrix to use, from 0 to k-1.
Note
symbolic_fact_MMt() and numeric_fact_MMt() must have been previously called.
void SparseChol::numeric_solve_MMt_diag ( double *  rhs,
int  i_k 
)
private

Numerical solve for IDENTITY, IDTY_IDTY, DIAGONAL, DIAG_DIAG (M*M' is diagonal)

Solves (M*Theta(i_k)*M')*sol = rhs, when (M*Theta(i_k)*M' is a diagonal matrix. This holds for matrices: IDENTITY, IDTY_IDTY, DIAGONAL, DIAG_DIAG

On output, rhs contains the solution. Previous calls to symbolic_fact_MMt and num_fact required for a correct "factorization" of M*Theta(i_k)*M (i.e., storage of in lnz(i_k) of the diagonal M*Theta(i_k)*M matrix)

void SparseChol::numeric_solve_MMt_sprsblkllt ( double *  rhs,
int  i_k = 0,
WHO_PERMUTES  whoperm = (WHO_PERMUTES)NULL 
)
private

Numerical solve for GENERAL and NETWORK matrices (require sprsblkllt)

Solves (ATheta(i_k)A')*sol = rhs. On output, rhs contains the solution. Previous calls to symbolic_fact_MMt and num_fact required for a correct factorization of ATheta(i_k)A'

whoperm: if CHOLESKY, then the reordering by pfa/ipfa must internally be applied to rhs by the routine; if USER_OF_CHOLESKY, then the user already provides rhs in the permutted order given by pfa/ipfa.

int i_k: which matrix to use, from 0 to k-1

void SparseChol::reset ( CHOL_SOLVER  chol_solver = (CHOL_SOLVER)NULL,
TYPE_MATRIX  type_matrix = (TYPE_MATRIX)NULL,
TYPE_ORIENTATION  type_orientation = ORIENTED 
)

Like calling destructor+constructor.

Parameters
chol_solverCholesky solver to be used
type_matrixType of matrix (general, network, IDTY, ...)
type_orientation(Only for network matrices) Whether arcs are oriented or nonoriented.
void SparseChol::symbolic_AThetaAt_A ( )
private

Symbolic factorization of AThetaA' or symmetric A.

It only requires xadj, ajcncy, which must have been previously created. Uses sfinit() and symfct() of sprsblkllt. Also uses bfinit() to prepare for forthcoming numerical factorizations.

void SparseChol::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 
)
inline

Computes symmetric ordering and symbolic factorization of A for a general symmetric matrix.

A is general symmetric upper triangular (with diagonal) matrix in compressed rowwise order

Parameters
mNumber of rows/columns of A.
nzNumber of nonzeros in A.
icola(nz)Column of each element of A (only diagonal and upper triangle)
inirowa(m+1)Pointers to first by row in icola().
prov_pfa(m),prov_ipfa(m)Ordering provided by user, to be used; the routine then only performs the symbolic factorization using this ordering. If NULL, then no ordering is provided and it will be computed. NULL is the default value if the parameter not provided. No check performed about the correctness of the ordering provided.
prov_pfa_numberingNumbering of provided ordering.
void SparseChol::symbolic_fact_M ( int  m)

Computes symbolic factorization of DIAGONAL matrix A.

When matrix is DIAG just allocate some vectors and fill some minimum information.

Parameters
mDimension of A.
void SparseChol::symbolic_fact_M_sprsblkllt_general ( int  nz,
int *  icola,
int *  inirowa 
)
private

Computes symmetric ordering and symbolic factorization of symmetric upper triangular A using Ng-Peyton package.

A is general symmetric (diagonal and upper triangle) matrix in compressed rowwise order

Parameters
nzNumber of nonzeros in A (diagonal and upper triangle).
icola(nz)Column of each element of A.
inirowa(m+1)Pointers to first by row in icola().
void SparseChol::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 
)
inline

Computes symmetric ordering and symbolic factorization of AA' for a general matrix A.

A is general matrix in compressed rowwise order.

Parameters
mNumber of rows of A.
nNumber of columns of A.
nzNumber of nonzeros in A.
icola(nz)Column of each element of A.
inirowa(m+1)Pointers to first by row in icola().
a(nz)Matrix elements.
kNumber of matrices with same topology to be factorized.
prov_pfa(m),prov_ipfa(m)Ordering provided by user, to be used; the routine then only performs the symbolic factorization using this ordering. If NULL, then no ordering is provided and it will be computed. NULL is the default value if the parameter not provided. No check performed about the correctness of the ordering provided.
prov_pfa_numberingNumbering of provided ordering.
void SparseChol::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 
)
inline

Computes symmetric ordering and symbolic factorization of AA' for a network matrix A.

A is node-arc incidence matrix (last node not considered)

Parameters
nnuNumber of nodes.
narNumber of arcs.
src(nar)Source node for each arc.
dst(nar)Destination node for each arc.
kNumber of matrices with same topology to be factorized.
prov_pfa(m),prov_ipfa(m)Ordering provided by user, to be used; the routine then only performs the symbolic factorization using this ordering. If NULL, then no ordering is provided and it will be computed. NULL is the default value if the parameter not provided. No check performed about the correctness of the ordering provided.
prov_pfa_numberingNumbering of provided ordering.
void SparseChol::symbolic_fact_MMt ( int  m,
int  k = 1 
)

Computes symbolic factorization of AA' for an IDENTITY or IDTY_IDTY matrix A.

When matrix is IDENTITY or IDTY_IDTY, just allocate some vectors and fill some minimum information

Parameters
mDimension of A.
kNumber of matrices with same topology to be factorized.
void SparseChol::symbolic_fact_MMt ( int  m,
double *  d1_in,
int  k = 1 
)

Computes symbolic factorization of AA' for a DIAG matrix A.

When matrix is DIAG just copy d1, allocate some vectors and fill some minimum information.

Parameters
mDimension of A.
d1_in(m)Vector of diagonal elements.
kNumber of matrices with same topology to be factorized.
Note
d1_in(m) is copied, and not free'd.
void SparseChol::symbolic_fact_MMt ( int  m,
double *  d1_in,
double *  d2_in,
int  k = 1 
)

Computes symbolic factorization of AA' for a DIAG_DIAG matrix A.

When matrix is [D1 D2] just copy d1,d2, allocate some vectors and fill some minimum information.

Parameters
mDimension of A.
d1_in(m)Vector of diagonal elements of D1.
d2_in(m)Vector of diagonal elements of D2.
kNumber of matrices with same topology to be factorized.
Note
d1_in(m), d2_in(m) are copied, and not free'd.
void SparseChol::symbolic_fact_MMt_sprsblkllt_general ( int  n,
int  nz,
int *  icola,
int *  inirowa,
double *  a 
)
private

Computes symmetric ordering and symbolic factorization for AA' using Ng-Peyton package for a general matrix.

A is general matrix in compressed rowwise order

Parameters
nNumber of columns of A.
nzNumber of nonzeros in A.
icola(nz)Column of each element of A.
inirowa(m+1)Pointers to first by row in icola().
a(nz)Matrix elements.
void SparseChol::symbolic_fact_MMt_sprsblkllt_network ( int *  src,
int *  dst 
)
private

Computes symmetric ordering and symbolic factorization for AA' using Ng-Peyton package for a network matrix.

Routines for both ORIENTED and NONORIENTED networks For ORIENTED the matrix is A=N For NONORIENTED the matrix is A=[N -N]. A*A' is then N*N'+(-N)(-N')= 2*NN'. Then the the topology of A*A' is the same for ORIENTED and NONORIENTED routines. Most routines are thus the same for the oriented and nonoriented case. The only difference is when computing A*Theta*A':

  • for oriented network: A*Theta*A'=N*Theta*N'
  • for nonoriented network: A*Theta*A'=[N -N]*diag(Theta^+, Theta^-)*[N -N']'= N*(Theta^+)*N'+N*(Theta^-)*N= N*(Theta^+ + Theta^-)*N'.

    Routine get_ilnz_network() computes indices for fast filling of lnz always considering the oriented network N. Therefore, the routine numeric_fact_MMt_sprsblkllt_network has to deal with both the oriented and nonoriented case later.

A is node-arc incidence matrix (last node not considered)

Parameters
src(nar)Source node for each arc.
dst(nar)Destination node for each arc,

Member Data Documentation

int* SparseChol::ajcncy
private

ajcncy(2*(njka-m)): 2*arcs in graph of A*A'/A without diagonal entries (i,i) of A*A'/A.

int* SparseChol::arck_l
private

arck_l(nar): arcs of inik, dst2.

int* SparseChol::arcl_k
private

arcl_k(nar): arcs of inil, src2.

CHOL_SOLVER SparseChol::chol_solver
private

Cholesky solver to be used.

int* SparseChol::colcnt
private

colnct(m): elements per column in factorization.

double* SparseChol::d1
private

Diagonal of DIAG or first diagonal of DIAG_DIAG [D1 D2] matrices.

double* SparseChol::d2
private

Second diagonal of DIAG_DIAG [D1 D2] matrices.

int* SparseChol::dst2
private

dst2(nar): destination nodes of arcs in inik.

int* SparseChol::ia
private

ia(m+1): pointer to first in ja().

When A is 'GENERAL' matrix:
ia,ja,ka,nla,la and va are internal indices for efficiently computing A*Theta*A' and symbolic factorization using the procedure described in "Monma, C. L. and A.J. Morton. 1987. Computational experience with a dual affine variant of Karmarkar's method for linear programming. Operations Research Letters V.6 n.6".

ia(m+1): pointer to first in ja(.)
ja(njka): pair of rows (ia,ja) of A that generate arc in graph of A*A'
ka(njka+1): pointer to first in la(.) and va(.)
nla: size of la() and va()
la(nla): columns that intervene in arc (i,j) in graph of A*A'
va(nla): premultiplied values associated to column of la(.)

int* SparseChol::iarc
private

iarc(maxiarc): indices to arcs intervening in each element nonzero of lnz().

int* SparseChol::ifillin
private

ifillin(maxfillin): fill-in positions in lnz(.) to be cleaned before filling lnz() with A*Theta*A'.

int* SparseChol::ilnz
private

ilnz(njka): indices to lnz(.) of each nonzero element of A*Theta*A'.

Variables with indices for fast filling of lnz().
If A is 'GENERAL' matrix:
ilnz(njka): indices to lnz(.) of each nonzero element of A*Theta*A'/A.
ifillin(maxfillin): fill-in positions in lnz(.) to be cleaned before filling lnz() with A*Theta*A'.

If A is 'NETWORK' matrix:
iarc(maxiarc): indices to arcs intervening in each element nonzero of lnz().
maxiarc: size of iarc.
ini_arc(njka+1): begin of arcs in iarc for each element of lnz().

int* SparseChol::ini_arc
private

ini_arc(njka+1): begin of arcs in iarc for each element of lnz().

int* SparseChol::inik
private

inik(nnu+1): pointers to first arc sorted by source (compressed form).

int* SparseChol::inil
private

inil(nnu+1): pointers to first arc sorted by destination (compressed form).

int* SparseChol::ipfa
private

ipfa[k]= i : i is current row (after permutation) of original row k.

int SparseChol::iwsiz
private

Size of integer working space, required by sprsblkllt.

int* SparseChol::ja
private

ja(njka): pair of rows (ia,ja) of A that generate arc in graph of A*A' .

int SparseChol::k
private

Number of matrices A*A' with same topology to be factorized, see below explanation of lnz() and pnlz() (k=1 for A)

int* SparseChol::ka
private

ka(njka+1): pointer to first in la() and va().

int* SparseChol::la
private

la(nla): columns that intervene in arc (i,j) in graph of A*A'.

int* SparseChol::lindx
private

lindx(maxsub): compressed form of sub and diagonal of factorization.

double* SparseChol::lnz
private

lnz(k*maxlnz): sub(upper) and diagonal elements of ATheta(k)A' (and its Cholesky factorization, once performed).

int SparseChol::m
private

Dimension of system A*A'/A : m x m matrix.

int SparseChol::maxfillin
private

Nonzeros in factorization of A*A'/A due to fill-in.

int SparseChol::maxiarc
private

Size of iarc.

int SparseChol::maxlnz
private

Total nonzeros in factorization of A*A'/A (= maxfillin+njka).

int SparseChol::maxsub
private

Maximum size of some vectors, computed by sprsblkllt.

int SparseChol::maxsuper
private

Maximum size of some vectors, computed by sprsblkllt.

constexpr double SparseChol::MIN_PIVOT = 0.0
staticprivate

Threshold pivot value for Cholesky factorizations of "simple" (diag, idty...) matrices (no pivot less than or equal to MIN_PIVOT allowed)

int SparseChol::nar
private

Number of arcs.

int SparseChol::njka
private

Nonzeros in diagonal and subdiagonal of A*A'/A.

int SparseChol::nla
private

Size of la() and va().

int SparseChol::nnu
private

Number of nodes.

When A is 'NETWORK' matrix:
nnu, nar, *inik, *dst2, *arck_l, *inil, *src2, *arcl_k are internal indices and variables for efficiently computing A*Theta*A' and symbolic factorization.

int SparseChol::num_semidef_matrix
private

Number of semidefinite matrices found during factorizations.

int SparseChol::num_zero_pivots
private

Number of zero pivots found during factorizations.

double* SparseChol::permrhs
private

Auxliary array to store the permutation of the rhs of the system to be solved.

int* SparseChol::pfa
private

pfa[i]= k : original row k is now (after permutation) in position i

Vector pfa of direct row permutation of matrix A*A', and its inverse (ipfa) computed by Cholesky solver.

pfa[i]= k : original row k is now (after permutation) in position i.
ipfa[k]= i : i is current row (after permutation) of original row k.
Both pfa and ipfa needed because permutation may be nonsymmetric.
pfa_numbering for either C-0 or Fortran-1 numbering style.
permrhs(): auxiliary array, stores the permutation of the rhs when solving the system.

NUMBERING SparseChol::pfa_numbering
private

Numbering style of pfa (used in C and Fortran code).

double** SparseChol::plnz
private

plnz(k): pointers to lnz(i).

This is used when k matrices with same structure must be factorized; their share symbolic factorization, the only difference is in lnz(); we need a different lnz() for each of them; i=0..k-1, so first matrix is matrix 0, last matrix is k-1.

constexpr double SparseChol::POSITIVE_PIVOT = 1.0E+128
staticprivate

New value for pivots below MIN_PIVOT in Cholesky factorizations of "simple" matrices (diag, idty...)

int* SparseChol::snode
private

snode(m): supernode of each column.

int* SparseChol::split
private

split(m): splitting of supernodes for exploiting cache memory.

int* SparseChol::src2
private

src2(nar): source nodes of arcs in inik

int SparseChol::tmpsiz
private

Size of temporary working space, required by sprsblkllt.

TYPE_MATRIX SparseChol::type_matrix
private

Type of matrix (general, network, identity, etc),.

TYPE_ORIENTATION SparseChol::type_orientation
private

Type of arc orientation for network matrices (by default, oriented).

double* SparseChol::va
private

va(nla): premultiplied values associated to column of la().

int* SparseChol::workspak
private
int* SparseChol::xadj
private

xadj(m+1): indices to adjacency matrix ajcncy().

Variables needed by sprsblkllt:

xadj(m+1): indices to adjacency matrix ajcncy().
ajcncy(2*(njka-m)): 2*arcs in graph of A*A'/A without diagonal entries (i,i) of A*A'/A.
colnct(m): elements per column in factorization.
snode(m): supernode of each column.
xsuper(m+1): pointer to first column of supernode.
xlindx(maxsuper): indices to lindx (compressed form).
lindx(maxsub): compressed form of sub and diagonal of factorization.
split(m): splitting of supernodes for exploiting cache memory.
xlnz(m+1): indices to first column/row in lnz(.).
lnz(k*maxlnz): sub(upper) and diagonal elements of ATheta(k)A' (and its Cholesky factorization, once performed).
plnz(k): pointer to lnz(i); this is used when k matrices with same structure must be factorized; their share symbolic factorization, thes only difference is in lnz(); we need a different lnz() for each of them; i=0..k-1, so first matrix is matrix 0, last matrix is k-1.

int* SparseChol::xlindx
private

xlindx(maxsuper): indices to lindx (compressed form).

int* SparseChol::xlnz
private

xlnz(m+1): indices to first column/row in lnz(.).

int* SparseChol::xsuper
private

xsuper(m+1): pointer to first column of supernode.


The documentation for this class was generated from the following files: