// hnummatrix.h
#pragma once
#ifndef __HNUMMATRIX_H__
#define __HNUMMATRIX_H__
#include "hnumdef.h"
namespace harlinn
{
namespace numerics
{
enum class MatrixStorageType
{
// column-wise, no supernode
NC,
// column-wise, column-permuted, no supernode
// (The consecutive columns of nonzeros, after permutation, may not be stored contiguously.)
NCP,
// row-wize, no supernode
NR,
// column-wise, supernode
SC,
// supernode, column-wise, permuted
SCP,
// row-wise, supernode
SR,
// Fortran style column-wise storage for dense matrix
DN,
// distributed compressed row format
NR_loc
};
enum class MatrixDataType
{
// float
S,
// double
D,
// std::complex<float>
C,
// std::complex<double>
Z
};
enum class MatrixType
{
// general
GE,
// lower triangular, unit diagonal
TRLU,
// upper triangular, unit diagonal
TRUU,
// lower triangular
TRL,
// upper triangular
TRU,
// symmetric, store lower half
SYL,
// symmetric, store upper half
SYU,
// Hermitian, store lower half
HEL,
// Hermitian, store upper half
HEU
};
/* Stype == SLU_NC (Also known as Harwell-Boeing sparse matrix format) */
template<typename T, typename int_t = int>
struct NCFormat
{
static MatrixStorageType StorageType()
{
return MatrixStorageType::NC;
}
// number of nonzeros in the matrix
int_t nnz;
// pointer to array of nonzero values, packed by column
T *nzval;
// pointer to array of row indices of the nonzeros
int_t *rowind;
// pointer to array of beginning of columns in nzval[] and rowind[]
int_t *colptr;
// Note:
// Zero-based indexing is used; colptr[] has ncol+1 entries,
// the last one pointing beyond the last column,
// so that colptr[ncol] = nnz.
};
/* Stype == SLU_NR */
template<typename T, typename int_t = int>
struct NRFormat
{
static MatrixStorageType StorageType()
{
return MatrixStorageType::NR;
}
// number of nonzeros in the matrix
int_t nnz;
// pointer to array of nonzero values, packed by raw
T* nzval;
// pointer to array of columns indices of the nonzeros
int_t *colind;
// pointer to array of beginning of rows in nzval[] and colind[]
int_t *rowptr;
// Note:
// Zero-based indexing is used;
// rowptr[] has nrow+1 entries, the last one pointing
// beyond the last row, so that rowptr[nrow] = nnz.
};
/* Stype == SLU_SC */
template<typename T, typename int_t = int>
struct SCFormat
{
static MatrixStorageType StorageType()
{
return MatrixStorageType::SC;
}
int_t nsuper; /* number of supernodes, minus 1 */
int_t nnz; /* number of nonzeros in the matrix */
T *nzval; /* pointer to array of nonzero values, packed by column */
int_t *nzval_colptr;/* pointer to array of beginning of columns in nzval[] */
int_t *rowind; /* pointer to array of compressed row indices of
rectangular supernodes */
int_t *rowind_colptr;/* pointer to array of beginning of columns in rowind[] */
int_t *col_to_sup; /* col_to_sup[j] is the supernode number to which column
j belongs; mapping from column to supernode number. */
int_t *sup_to_col; /* sup_to_col[s] points to the start of the s-th
supernode; mapping from supernode number to column.
e.g.: col_to_sup: 0 1 2 2 3 3 3 4 4 4 4 4 4 (ncol=12)
sup_to_col: 0 1 2 4 7 12 (nsuper=4) */
/* Note:
Zero-based indexing is used;
nzval_colptr[], rowind_colptr[], col_to_sup and
sup_to_col[] have ncol+1 entries, the last one
pointing beyond the last column.
For col_to_sup[], only the first ncol entries are
defined. For sup_to_col[], only the first nsuper+2
entries are defined. */
};
/* Stype == SLU_SCP */
template<typename T, typename int_t = int>
struct SCPFormat
{
static MatrixStorageType StorageType()
{
return MatrixStorageType::SCP;
}
/* number of supernodes */
int_t nsuper;
/* number of nonzeros in the matrix */
int_t nnz;
/* pointer to array of nonzero values, packed by column */
T* nzval;
/* nzval_colbeg[j] points to beginning of column j in nzval[] */
int_t *nzval_colbeg;
/* nzval_colend[j] points to one past the last element of column j in nzval[] */
int_t *nzval_colend;
/* pointer to array of compressed row indices of rectangular supernodes */
int_t *rowind;
/* rowind_colbeg[j] points to beginning of column j in rowind[] */
int_t *rowind_colbeg;
/* rowind_colend[j] points to one past the last element of column j in rowind[] */
int_t *rowind_colend;
/* col_to_sup[j] is the supernode number to which column j belongs; mapping from column to supernode. */
int_t *col_to_sup;
/* sup_to_colbeg[s] points to the start of the s-th supernode; mapping from supernode to column.*/
int_t *sup_to_colbeg;
/* sup_to_colend[s] points to one past the end of the
s-th supernode; mapping from supernode number to
column.
e.g.: col_to_sup: 0 1 2 2 3 3 3 4 4 4 4 4 4 (ncol=12)
sup_to_colbeg: 0 1 2 4 7 (nsuper=4)
sup_to_colend: 1 2 4 7 12
*/
int_t *sup_to_colend;
/* Note:
Zero-based indexing is used;
nzval_colptr[], rowind_colptr[], col_to_sup and
sup_to_col[] have ncol+1 entries, the last one
pointing beyond the last column.
*/
};
/* Stype == SLU_NCP */
template<typename T, typename int_t = int>
struct NCPFormat
{
static MatrixStorageType StorageType()
{
return MatrixStorageType::NCP;
}
/* number of nonzeros in the matrix */
int_t nnz;
/* pointer to array of nonzero values, packed by column */
T* nzval;
/* pointer to array of row indices of the nonzeros */
/* Note: nzval[]/rowind[] always have the same length */
int_t *rowind;
/* colbeg[j] points to the beginning of column j in nzval[] and rowind[] */
int_t *colbeg;
/* colend[j] points to one past the last element of columnj in nzval[] and rowind[] */
int_t *colend;
/* Note:
Zero-based indexing is used;
The consecutive columns of the nonzeros may not be
contiguous in storage, because the matrix has been
postmultiplied by a column permutation matrix. */
};
/* Stype == SLU_DN */
template<typename T, typename int_t = int>
struct DNFormat
{
static MatrixStorageType StorageType()
{
return MatrixStorageType::DN;
}
/* leading dimension */
int_t lda;
/* array of size lda*ncol to represent a dense matrix */
void* nzval;
};
/* Stype == SLU_NR_loc (Distributed Compressed Row Format) */
template<typename T, typename int_t = int>
struct NRFormatLoc
{
static MatrixStorageType StorageType()
{
return MatrixStorageType::NR_loc;
}
/* number of nonzeros in the local submatrix */
int_t nnz_loc;
/* number of rows local to this processor */
int_t m_loc;
/* global index of the first row */
int_t fst_row;
/* pointer to array of nonzero values, packed by row */
T* nzval;
/* pointer to array of beginning of rows in nzval[] and colind[] */
int_t *rowptr;
/* pointer to array of column indices of the nonzeros */
int_t *colind;
/* Note:
Zero-based indexing is used;
rowptr[] has n_loc + 1 entries, the last one pointing
beyond the last row, so that rowptr[n_loc] = nnz_loc.
*/
};
template<typename T, typename int_t = int>
class Matrix
{
class Data : public T
{
long referenceCount;
MatrixType matrixType;
// number of rows
int_t nrow;
// number of columns
int_t ncol;
public:
Data()
: referenceCount(1),
nrow(0),
ncol(0)
{
T* ptr = this;
memset(ptr,0,sizeof(T));
}
long AddRef()
{
return InterlockedIncrement(&referenceCount);
}
long Release()
{
long result = InterlockedDecrement(&referenceCount);
if(!result)
{
delete this;
}
return result;
}
static MatrixStorageType StorageType()
{
return T::StorageType();
}
};
typedef Data data_t;
mutable data_t* data;
data_t* Data()
{
if(!data)
{
data = new data_t();
}
return data;
}
const data_t* Data() const
{
if(!data)
{
data = new data_t();
}
return data;
}
public:
Matrix()
: data(nullptr)
{}
Matrix(const Matrix& other)
: data(other.data)
{
if(data)
{
data->AddRef();
}
}
Matrix(Matrix&& other)
: data(other.data)
{
other.data = nullptr;
}
~Matrix( )
{
if(data)
{
data->Release();
}
}
Matrix& operator = (const Matrix& other)
{
if(data != other.data)
{
if(data)
{
data->Release();
}
data = other.data;
if(data)
{
data->AddRef();
}
}
return *this;
}
Matrix& operator = (Matrix&& other)
{
if(this != &other)
{
if(data)
{
data->Release();
}
data = other.data;
if(data)
{
data->AddRef();
}
}
return *this;
}
static MatrixStorageType StorageType()
{
return T::StorageType();
}
};
typedef Matrix< NCFormat< float > > MatrixNCS;
typedef Matrix< NCFormat< double > > MatrixNCD;
typedef Matrix< NCFormat< std::complex<float> > > MatrixNCC;
typedef Matrix< NCFormat< std::complex<double> > > MatrixNCZ;
typedef Matrix< NCPFormat< float > > MatrixNCPS;
typedef Matrix< NCPFormat< double > > MatrixNCPD;
typedef Matrix< NCPFormat< std::complex<float> > > MatrixNCPC;
typedef Matrix< NCPFormat< std::complex<double> > > MatrixNCPZ;
typedef Matrix< NRFormat< float > > MatrixNRS;
typedef Matrix< NRFormat< double > > MatrixNRD;
typedef Matrix< NRFormat< std::complex<float> > > MatrixNRC;
typedef Matrix< NRFormat< std::complex<double> > > MatrixNRZ;
typedef Matrix< SCFormat< float > > MatrixSCS;
typedef Matrix< SCFormat< double > > MatrixSCD;
typedef Matrix< SCFormat< std::complex<float> > > MatrixSCC;
typedef Matrix< SCFormat< std::complex<double> > > MatrixSCZ;
typedef Matrix< SCPFormat< float > > MatrixSCPS;
typedef Matrix< SCPFormat< double > > MatrixSCPD;
typedef Matrix< SCPFormat< std::complex<float> > > MatrixSCPC;
typedef Matrix< SCPFormat< std::complex<double> > > MatrixSCPZ;
typedef Matrix< DNFormat< float > > MatrixDNS;
typedef Matrix< DNFormat< double > > MatrixDND;
typedef Matrix< DNFormat< std::complex<float> > > MatrixDNC;
typedef Matrix< DNFormat< std::complex<double> > > MatrixDNZ;
typedef Matrix< NRFormatLoc< float > > MatrixNRLocS;
typedef Matrix< NRFormatLoc< double > > MatrixNRLocD;
typedef Matrix< NRFormatLoc< std::complex<float> > > MatrixNRLocC;
typedef Matrix< NRFormatLoc< std::complex<double> > > MatrixNRLocZ;
};
};
#endif // __HNUMMATRIX_H__
Submit an article or tip