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 AThetaA'.

void	numeric_solve_MMt (double *rhs, int i_k=0, WHO_PERMUTES whoperm=(WHO_PERMUTES) NULL)
	Solve systems with matrix AThetaA' 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 AThetaA' when A is general matrix.

void	numeric_fact_MMt_sprsblkllt_network (double *Theta, int i_k)
	Computes numerical factorization of AThetaA' when A is a network.

void	numeric_fact_MMt_identity (double *Theta, int i_k)
	Numerical factorization of IThetaI'.

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 DThetaD' (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)): 2arcs in graph of AA'/A without diagonal entries (i,i) of AA'/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 AThetaA'.

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

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

inp	ia Vector of beginning of rows of A*A'/A(compressed form)
inp_ja	Vector 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

icola	Vector of column indices.
inirowa	Vector of beginning of rows (compressed form)
a	Vector 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

src	Vector of arc sources.
dst	Vector 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

inp	ia Vector of beginning of rows of A*A'/A(compressed form)
inp_ja	Vector 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_solver	Cholesky solver to be used
type_matrix	Type 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

d1	diagonal 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(nz	or 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 AThetaA', of size number of columns of A.
i_k	Which 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 AThetaA', of size number of columns of A
i_k	Which 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

Theta	Diagonal Theta matrix.
i_k	Which 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

rhs	On input vector rhs; on output the solution vector.
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.

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

rhs	on 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

rhs	On input vector rhs; on output the solution vector.
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.
i_k	Which 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_solver	Cholesky solver to be used
type_matrix	Type 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

m	Number of rows/columns of A.
nz	Number 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_numbering	Numbering 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

m	Dimension 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

nz	Number 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

m	Number of rows of A.
n	Number of columns of A.
nz	Number 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.
k	Number 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_numbering	Numbering 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

nnu	Number of nodes.
nar	Number of arcs.
src(nar)	Source node for each arc.
dst(nar)	Destination node for each arc.
k	Number 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_numbering	Numbering 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

m	Dimension of A.
k	Number 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

m	Dimension of A.
d1_in(m)	Vector of diagonal elements.
k	Number 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

m	Dimension of A.
d1_in(m)	Vector of diagonal elements of D1.
d2_in(m)	Vector of diagonal elements of D2.
k	Number 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

n	Number of columns of A.
nz	Number 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:

/home/jcastro2/intpoint/BlockIP/BlockIP/src_blockip-v2/SparseChol.h
/home/jcastro2/intpoint/BlockIP/BlockIP/src_blockip-v2/SparseChol.C
/home/jcastro2/intpoint/BlockIP/BlockIP/src_blockip-v2/SparseCholSprsblkllt.C

Classes

Public Member Functions

Private Member Functions

Static Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation