@ -21,84 +21,87 @@ namespace mrmc {
namespace sparse {
/ / ! A sparse Matrix for DTMCs with a constant number of non - zero entries on the non - diagonal fields .
/*!
The sparse matrix used for calculation on the DTMCs .
Addressing is NOT zero - based ! The valid range for getValue and addNextValue is 1 to rows ( first argument to constructor ) .
* A sparse matrix class with a constant number of non - zero entries on the non - diagonal fields
* and a seperate dense storage for the diagonal elements .
* NOTE : Addressing * is * zero - based , so the valid range for getValue and addNextValue is 0. . ( rows - 1 )
* where rows is the first argument to the constructor .
*/
template < class T >
class StaticSparseMatrix {
public :
/ / ! Status enum
/*!
An E num representing the internal state of the Matrix .
After creating the Matrix using the Constructor , the Object is in state UnInitialized . After calling initialize ( ) , that state changes to Initialized and after all entries have been entered and finalize ( ) has been called , to ReadReady .
Should a critical error occur in any of the former functions , the state will change to Error .
@ see getState ( )
@ see isReadReady ( )
@ see isInitialized ( )
@ see hasError ( )
* An e num representing the internal state of the Matrix .
* After creating the Matrix using the Constructor , the Object is in state UnInitialized . After calling initialize ( ) , that state changes to Initialized and after all entries have been entered and finalize ( ) has been called , to ReadReady .
* Should a critical error occur in any of the former functions , the state will change to Error .
* @ see getState ( )
* @ see isReadReady ( )
* @ see isInitialized ( )
* @ see hasError ( )
*/
enum MatrixStatus {
Error = - 1 ,
UnInitialized = 0 ,
Initialized = 1 ,
ReadReady = 2 ,
Error = - 1 , UnInitialized = 0 , Initialized = 1 , ReadReady = 2 ,
} ;
/ / ! Constructor
/*!
\ param rows Row - Count and therefore column - count of the square matrix
* Constructs a sparse matrix object with the given number of rows .
* @ param rows The number of rows of the matrix
*/
StaticSparseMatrix ( uint_fast64_t rows ) {
/ / Using direct access instead of setState ( ) because of undefined initialization value
/ / setState ( ) stays in Error should Error be the current value
internal_status = MatrixStatus : : UnInitialized ;
current_size = 0 ;
value_storage = NULL ;
diagonal_storage = NULL ;
column_indications = NULL ;
row_indications = NULL ;
row_count = rows ;
non_zero_entry_count = 0 ;
/ / initialize ( rows , non_zero_entries ) ;
}
StaticSparseMatrix ( uint_fast64_t rows )
: internal_status ( MatrixStatus : : UnInitialized ) ,
current_size ( 0 ) , last_row ( 0 ) , value_storage ( nullptr ) ,
diagonal_storage ( nullptr ) , column_indications ( nullptr ) ,
row_indications ( nullptr ) , row_count ( rows ) , non_zero_entry_count ( 0 ) { }
/ / ! Copy Constructor
/*!
Copy Constructor . Creates an exact copy of the source sparse matrix ssm . Modification of either matrix does not affect the other .
@ param ssm A reference to the matrix that sh ould be copied from
* Copy Constructor . Performs a deep copy of the given sparse matrix .
* @ param ssm A reference to the matrix to be copied .
*/
StaticSparseMatrix ( const StaticSparseMatrix < T > & ssm ) : internal_status ( ssm . internal_status ) , current_size ( ssm . current_size ) , row_count ( ssm . row_count ) , non_zero_entry_count ( ssm . non_zero_entry_count )
{
pantheios : : log_DEBUG ( " StaticSparseMatrix::CopyCTor: Using Copy() Ctor. " ) ;
StaticSparseMatrix ( const StaticSparseMatrix < T > & ssm )
: internal_status ( ssm . internal_status ) ,
current_size ( ssm . current_size ) , last_row ( ssm . last_row ) ,
row_count ( ssm . row_count ) ,
non_zero_entry_count ( ssm . non_zero_entry_count ) {
pantheios : : log_DEBUG ( " StaticSparseMatrix::CopyConstructor: Using copy constructor. " ) ;
/ / Check whether copying the matrix is safe .
if ( ! ssm . hasError ( ) ) {
pantheios : : log_ERROR ( " StaticSparseMatrix::CopyCtor: Throwing invalid_argument: Can not Copy from matrix in Error state. " ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::CopyCtor: Throwing invalid_argument: Can not Copy from matrix in e rror state. " ) ;
throw mrmc : : exceptions : : invalid_argument ( ) ;
} else {
/ / Try to prepare the internal storage and throw an error in case
/ / of a failure .
if ( ! prepareInternalStorage ( ) ) {
pantheios : : log_ERROR ( " StaticSparseMatrix::CopyCtor: Throwing bad_alloc: memory allocation failed. " ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::CopyConstruc tor: Throwing bad_alloc: memory allocation failed. " ) ;
throw std : : bad_alloc ( ) ;
} else {
/ / Now that all storages have been prepared , copy over all
/ / elements . Start by copying the elements of type value and
/ / copy them seperately in order to invoke copy their copy
/ / constructor . This may not be necessary , but it is safer to
/ / do so in any case .
for ( uint_fast64_t i = 0 ; i < non_zero_entry_count ; + + i ) {
/ / use T ( ) to force use of the copy constructor for complex T types
value_storage [ i ] = T ( ssm . value_storage [ i ] ) ;
}
for ( uint_fast64_t i = 0 ; i < = row_count ; + + i ) {
/ / same here
/ / use T ( ) to force use of the copy constructor for complex T types
diagonal_storage [ i ] = T ( ssm . diagonal_storage [ i ] ) ;
}
/ / The elements that are not of the value type but rather the
/ / index type may be copied with memcpy .
memcpy ( column_indications , ssm . column_indications , sizeof ( column_indications [ 0 ] ) * non_zero_entry_count ) ;
memcpy ( row_indications , ssm . row_indications , sizeof ( row_indications [ 0 ] ) * ( row_count + 1 ) ) ;
}
}
}
/ / ! Destructor
/*!
* Destructor . Performs deletion of the reserved storage arrays .
*/
~ StaticSparseMatrix ( ) {
setState ( MatrixStatus : : UnInitialized ) ;
if ( value_storage ! = NULL ) {
@ -119,14 +122,16 @@ class StaticSparseMatrix {
}
}
/ / ! Mandatory initialization of the matrix , variant for initialize ( ) , addNextValue ( ) and finalize ( )
/*!
Mandatory initialization of the matrix , must be called before using any other member function .
This version is to be used together with addNextValue ( ) .
For initialization from a Eigen SparseMatrix , use initialize ( Eigen : : SparseMatrix < T > & ) .
\ param non_zero_entries The exact count of entries that will be submitted through addNextValue * excluding * those on the diagonal ( A_ { i , j } with i = j )
* Initializes the sparse matrix with the given number of non - zero entries
* and prepares it for use with addNextValue ( ) and finalize ( ) .
* NOTE : Calling this method before any other member function is mandatory .
* This version is to be used together with addNextValue ( ) . For
* initialization from an Eigen SparseMatrix , use initialize ( Eigen : : SparseMatrix < T > & ) .
* @ param non_zero_entries
*/
void initialize ( uint_fast64_t non_zero_entries ) {
/ / Check whether initializing the matrix is safe .
if ( internal_status ! = MatrixStatus : : UnInitialized ) {
triggerErrorState ( ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::initialize: Throwing invalid state for status flag != 0 (is " , pantheios : : integer ( internal_status ) , " - Already initialized? " ) ;
@ -140,6 +145,8 @@ class StaticSparseMatrix {
pantheios : : log_ERROR ( " StaticSparseMatrix::initialize: Throwing invalid_argument: More non-zero entries than entries in target matrix " ) ;
throw mrmc : : exceptions : : invalid_argument ( " mrmc::StaticSparseMatrix::initialize: More non-zero entries than entries in target matrix " ) ;
} else {
/ / If it is safe , initialize necessary members and prepare the
/ / internal storage .
non_zero_entry_count = non_zero_entries ;
last_row = 0 ;
@ -153,72 +160,94 @@ class StaticSparseMatrix {
}
}
/ / ! Mandatory initialization of the matrix , variant for initialize ( ) , addNextValue ( ) and finalize ( )
/*!
Mandatory initialization of the matrix , must be called before using any other member function .
This version is to be used for initialization from a Eigen SparseMatrix , use initialize ( uint_fast32_t ) for addNextValue .
\ param eigen_sparse_matrix The Eigen Sparse Matrix to be copied / initialized from . MUST BE in compressed form !
* Initializes the sparse matrix with the given Eigen sparse matrix .
* NOTE : Calling this method before any other member function is mandatory .
* This version is only to be used when copying an Eigen sparse matrix . For
* initialization with addNextValue ( ) and finalize ( ) use initialize ( uint_fast32_t )
* instead .
* @ param eigen_sparse_matrix The Eigen sparse matrix to be copied .
* * NOTE * Has to be in compressed form !
*/
template < int _Options , typename _Index >
void initialize ( const Eigen : : SparseMatrix < T , _Options , _Index > & eigen_sparse_matrix ) {
/ / Throw an error in case the matrix is not in compressed format .
if ( ! eigen_sparse_matrix . isCompressed ( ) ) {
triggerErrorState ( ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::initialize: Throwing invalid_argument: eigen_sparse_matrix is not in Compressed form. " ) ;
throw mrmc : : exceptions : : invalid_argument ( " StaticSparseMatrix::initialize: Throwing invalid_argument: eigen_sparse_matrix is not in Compressed form. " ) ;
}
/ / Compute the actual ( i . e . non - diagonal ) number of non - zero entries .
non_zero_entry_count = getEigenSparseMatrixCorrectNonZeroEntryCount ( eigen_sparse_matrix ) ;
last_row = 0 ;
/ / Try to prepare the internal storage and throw an error in case of
/ / failure .
if ( ! prepareInternalStorage ( ) ) {
triggerErrorState ( ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::initialize: Throwing bad_alloc: memory allocation failed " ) ;
pantheios : : log_ERROR (
" StaticSparseMatrix::initialize: Throwing bad_alloc: memory allocation failed " ) ;
throw std : : bad_alloc ( ) ;
} else {
/ / easy case , we can simply copy the data
/ / RowMajor : Easy , ColMajor : Hmm . But how to detect ?
/ / Get necessary pointers to the contents of the Eigen matrix .
const T * valuePtr = eigen_sparse_matrix . valuePtr ( ) ;
const _Index * indexPtr = eigen_sparse_matrix . innerIndexPtr ( ) ;
const _Index * outerPtr = eigen_sparse_matrix . outerIndexPtr ( ) ;
const _Index entryCount = eigen_sparse_matrix . nonZeros ( ) ;
const _Index outerCount = eigen_sparse_matrix . outerSize ( ) ;
/ / If the given matrix is in RowMajor format , copying can simply
/ / be done by adding all values in order .
/ / Direct copying is , however , prevented because we have to
/ / separate the diagonal entries from others .
if ( isEigenRowMajor ( eigen_sparse_matrix ) ) {
/ / Easy case , all data can be copied with some minor changes .
/ / We have to separate diagonal entries from others
for ( _Index row = 1 ; row < = outerCount ; + + row ) {
for ( _Index col = outerPtr [ row - 1 ] ; col < outerPtr [ row ] ; + + col ) {
addNextValue ( row , indexPtr [ col ] + 1 , valuePtr [ col ] ) ;
/ / Because of the RowMajor format outerSize evaluates to the
/ / number of rows .
const _Index rowCount = eigen_sparse_matrix . outerSize ( ) ;
for ( _Index row = 0 ; row < rowCount ; + + row ) {
for ( _Index col = outerPtr [ row ] ; col < outerPtr [ row + 1 ] ;
+ + col ) {
addNextValue ( row , indexPtr [ col ] , valuePtr [ col ] ) ;
}
}
} else {
/ / temp copies , anyone ?
const _Index eigen_col_count = eigen_sparse_matrix . cols ( ) ;
const _Index eigen_row_count = eigen_sparse_matrix . rows ( ) ;
/ / initialise all column - start positions to known lower boundarys
_Index * positions = new _Index [ eigen_col_count ] ( ) ;
for ( _Index i = 0 ; i < eigen_col_count ; + + i ) {
const _Index entryCount = eigen_sparse_matrix . nonZeros ( ) ;
/ / Because of the ColMajor format outerSize evaluates to the
/ / number of columns .
const _Index colCount = eigen_sparse_matrix . outerSize ( ) ;
/ / Create an array to remember which elements have to still
/ / be searched in each column and initialize it with the starting
/ / index for every column .
_Index * positions = new _Index [ colCount ] ( ) ;
for ( _Index i = 0 ; i < colCount ; + + i ) {
positions [ i ] = outerPtr [ i ] ;
}
/ / Now copy the elements . As the matrix is in ColMajor format ,
/ / we need to iterate over the columns to find the next non - zero
/ / entry .
int i = 0 ;
int currentRow = 0 ;
int currentColumn = 0 ;
while ( i < entryCount ) {
if ( ( positions [ currentColumn ] < outerPtr [ currentColumn + 1 ] ) & & ( indexPtr [ positions [ currentColumn ] ] = = currentRow ) ) {
addNextValue ( currentRow + 1 , currentColumn + 1 , valuePtr [ positions [ currentColumn ] ] ) ;
/ / one more found
+ + i ;
/ / mark this position as " used "
+ + positions [ currentColumn ] ;
}
/ / advance to next column
+ + currentColumn ;
if ( currentColumn = = eigen_col_count ) {
/ / If the current element belongs the the current column ,
/ / add it in case it is also in the current row .
if ( ( positions [ currentColumn ] < outerPtr [ currentColumn + 1 ] )
& & ( indexPtr [ positions [ currentColumn ] ] = = currentRow ) ) {
addNextValue ( currentRow , currentColumn ,
valuePtr [ positions [ currentColumn ] ] ) ;
/ / Remember that we found one more non - zero element .
i + + ;
/ / Mark this position as " used " .
positions [ currentColumn ] + + ;
}
/ / Now we can advance to the next column and also row ,
/ / in case we just iterated through the last column .
currentColumn + + ;
if ( currentColumn = = colCount ) {
currentColumn = 0 ;
+ + currentRow ;
currentRow + + ;
}
}
delete [ ] positions ;
@ -227,234 +256,340 @@ class StaticSparseMatrix {
}
}
/ / ! Linear Setter for matrix entry A_ { row , col } to value
/*!
Linear Setter function for matrix entry A_ { row , col } to value . Must be called consecutively for each element in a row in ascending order of columns AND in ascending order of rows .
Diagonal entries may be set at any time .
* Sets the matrix element at the given row and column to the given value .
* NOTE : This is a linear setter . It must be called consecutively for each element ,
* row by row * and * column by column . Only diagonal entries may be set at any time .
* @ param row The row in which the matrix element is to be set .
* @ param col The column in which the matrix element is to be set .
* @ param value The value that is to be set .
*/
void addNextValue ( const uint_fast64_t row , const uint_fast64_t col , const T & value ) {
if ( ( row > row_count ) | | ( col > row_count ) | | ( row = = 0 ) | | ( col = = 0 ) ) {
/ / Check whether the given row and column positions are valid and throw
/ / error otherwise .
if ( ( row > row_count ) | | ( col > row_count ) ) {
triggerErrorState ( ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::addNextValue: Throwing out_of_range: row or col not in 1 .. rows (is " , pantheios : : integer ( row ) , " x " , pantheios : : integer ( col ) , " , max is " , pantheios : : integer ( row_count ) , " x " , pantheios : : integer ( row_count ) , " ). " ) ;
throw mrmc : : exceptions : : out_of_range ( " StaticSparseMatrix::addNextValue: row or col not in 1 .. rows " ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::addNextValue: Throwing out_of_range: row or col not in 0 .. rows (is " ,
pantheios : : integer ( row ) , " x " , pantheios : : integer ( col ) , " , max is " ,
pantheios : : integer ( row_count ) , " x " , pantheios : : integer ( row_count ) , " ). " ) ;
throw mrmc : : exceptions : : out_of_range ( " StaticSparseMatrix::addNextValue: row or col not in 0 .. rows " ) ;
}
if ( row = = col ) {
if ( row = = col ) { / / Set a diagonal element .
diagonal_storage [ row ] = value ;
} else {
} else { / / Set a non - diagonal element .
/ / If we switched to another row , we have to adjust the missing
/ / entries in the row_indications array .
if ( row ! = last_row ) {
for ( uint_fast64_t i = last_row ; i < row ; + + i ) {
for ( uint_fast64_t i = last_row + 1 ; i < = row ; + + i ) {
row_indications [ i ] = current_size ;
}
last_row = row ;
}
/ / Finally , set the element and increase the current size .
value_storage [ current_size ] = value ;
column_indications [ current_size ] = col ;
/ / Increment counter for stored elements
current_size + + ;
}
}
/*
* Finalizes the sparse matrix to indicate that initialization has been
* completed and the matrix may now be used .
*/
void finalize ( ) {
/ / Check whether it ' s safe to finalize the matrix and throw error
/ / otherwise .
if ( ! isInitialized ( ) ) {
triggerErrorState ( ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::finalize: Throwing invalid state for internal state not Initialized (is " , pantheios : : integer ( internal_status ) , " - Already finalized? " ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::finalize: Throwing invalid state for internal state not Initialized (is " ,
pantheios : : integer ( internal_status ) , " - Already finalized? " ) ;
throw mrmc : : exceptions : : invalid_state ( " StaticSparseMatrix::finalize: Invalid state for internal state not Initialized - Already finalized? " ) ;
} else if ( current_size ! = non_zero_entry_count ) {
triggerErrorState ( ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::finalize: Throwing invalid_state: Wrong call count for addNextValue " ) ;
throw mrmc : : exceptions : : invalid_state ( " StaticSparseMatrix::finalize: Wrong call count for addNextValue " ) ;
} else {
/ / Fill in the missing entries in the row_indications array .
/ / ( Can happen because of empty rows at the end . )
if ( last_row ! = row_count ) {
for ( uint_fast64_t i = last_row ; i < row_count ; + + i ) {
for ( uint_fast64_t i = last_row + 1 ; i < = row_count ; + + i ) {
row_indications [ i ] = current_size ;
}
}
row_indications [ row_count ] = non_zero_entry_count ;
/ / Set a sentinel element at the last position of the row_indications
/ / array . This eases iteration work , as now the indices of row i
/ / are always between row_indications [ i ] and row_indications [ i + 1 ] ,
/ / also for the first and last row .
row_indications [ row_count + 1 ] = non_zero_entry_count ;
setState ( MatrixStatus : : ReadReady ) ;
}
}
/ / ! Getter for saving matrix entry A_ { row , col } to target
/*!
Getter function for the matrix . This function does not check the internal status for errors for performance reasons .
\ param row 1 - based index of the requested row
\ param col 1 - based index of the requested column
\ param target pointer to where the result will be stored
\ return True iff the value was set , false otherwise . On false , 0 will be written to * target .
* Gets the matrix element at the given row and column to the given value .
* NOTE : This function does not check the internal status for errors for performance reasons .
* @ param row The row in which the element is to be read .
* @ param col The column in which the element is to be read .
* @ param target A pointer to the memory location where the read content is
* to be put .
* @ return True iff the value is set in the matrix , false otherwise .
* On false , 0 will be written to * target .
*/
inline bool getValue ( uint_fast64_t row , uint_fast64_t col , T * const target ) {
/ / Check for illegal access indices .
if ( ( row > row_count ) | | ( col > row_count ) ) {
pantheios : : log_ERROR ( " StaticSparseMatrix::getValue: row or col not in 0 .. rows (is " , pantheios : : integer ( row ) , " x " ,
pantheios : : integer ( col ) , " , max is " , pantheios : : integer ( row_count ) , " x " , pantheios : : integer ( row_count ) , " ). " ) ;
throw mrmc : : exceptions : : out_of_range ( " StaticSparseMatrix::getValue: row or col not in 0 .. rows " ) ;
return false ;
}
/ / Read elements on the diagonal directly .
if ( row = = col ) {
/ / storage is row_count + 1 large for direct access without the - 1
* target = diagonal_storage [ row ] ;
return true ;
}
if ( ( row > row_count ) | | ( col > row_count ) | | ( row = = 0 ) | | ( col = = 0 ) ) {
pantheios : : log_ERROR ( " StaticSparseMatrix::getValue: row or col not in 1 .. rows (is " , pantheios : : integer ( row ) , " x " , pantheios : : integer ( col ) , " , max is " , pantheios : : integer ( row_count ) , " x " , pantheios : : integer ( row_count ) , " ). " ) ;
throw mrmc : : exceptions : : out_of_range ( " StaticSparseMatrix::getValue: row or col not in 1 .. rows " ) ;
return false ;
}
uint_fast64_t row_start = row_indications [ row - 1 ] ;
uint_fast64_t row_end = row_indications [ row ] ;
/ / In case the element is not on the diagonal , we have to iterate
/ / over the accessed row to find the element .
uint_fast64_t row_start = row_indications [ row ] ;
uint_fast64_t row_end = row_indications [ row + 1 ] ;
while ( row_start < row_end ) {
/ / If the lement is found , write the content to the specified
/ / position and return true .
if ( column_indications [ row_start ] = = col ) {
* target = value_storage [ row_start ] ;
return true ;
}
/ / If the column of the current element is already larger than the
/ / requested column , the requested element cannot be contained
/ / in the matrix and we may therefore stop searching .
if ( column_indications [ row_start ] > col ) {
break ;
}
row_start + + ;
}
/ / Set 0 as the content and return false in case the element was not found .
* target = 0 ;
return false ;
}
/*!
* Returns the number of rows of the matrix .
*/
uint_fast64_t getRowCount ( ) const {
return row_count ;
}
/*!
* Returns a pointer to the value storage of the matrix . This storage does
* * not * include elements on the diagonal .
* @ return A pointer to the value storage of the matrix .
*/
T * getStoragePointer ( ) const {
return value_storage ;
}
/*!
* Returns a pointer to the storage of elements on the diagonal .
* @ return A pointer to the storage of elements on the diagonal .
*/
T * getDiagonalStoragePointer ( ) const {
return diagonal_storage ;
}
/*!
* Returns a pointer to the array that stores the start indices of non - zero
* entries in the value storage for each row .
* @ return A pointer to the array that stores the start indices of non - zero
* entries in the value storage for each row .
*/
uint_fast64_t * getRowIndicationsPointer ( ) const {
return row_indications ;
}
/*!
* Returns a pointer to an array that stores the column of each non - zero
* element that is not on the diagonal .
* @ return A pointer to an array that stores the column of each non - zero
* element that is not on the diagonal .
*/
uint_fast64_t * getColumnIndicationsPointer ( ) const {
return column_indications ;
}
/*!
* Checks whether the internal status of the matrix makes it ready for
* reading access .
* @ return True iff the internal status of the matrix makes it ready for
* reading access .
*/
bool isReadReady ( ) {
return ( internal_status = = MatrixStatus : : ReadReady ) ;
}
/*!
* Checks whether the matrix was initialized previously . The matrix may
* still require to be finalized , even if this check returns true .
* @ return True iff the matrix was initialized previously .
*/
bool isInitialized ( ) {
return ( internal_status = = MatrixStatus : : Initialized | | internal_status = = MatrixStatus : : ReadReady ) ;
}
/*!
* Returns the internal state of the matrix .
* @ return The internal state of the matrix .
*/
MatrixStatus getState ( ) {
return internal_status ;
}
/*!
* Checks whether the internal state of the matrix signals an error .
* @ return True iff the internal state of the matrix signals an error .
*/
bool hasError ( ) {
return ( internal_status = = MatrixStatus : : Error ) ;
}
/ / ! Converts this matrix to an equivalent sparse matrix in Eigens format .
/*!
Exports this sparse matrix to Eigens SparseMatrix format .
Required this matrix to be in the ReadReady state .
@ return The Eigen SparseMatrix
* Exports this sparse matrix to Eigens sparse matrix format .
* NOTE : this requires this matrix to be in the ReadReady state .
* @ return The sparse matrix in Eigen format .
*/
Eigen : : SparseMatrix < T , Eigen : : RowMajor , int_fast32_t > * toEigenSparseMatrix ( ) {
int_fast32_t eigenRows = static_cast < int_fast32_t > ( row_count ) ;
Eigen : : SparseMatrix < T , Eigen : : RowMajor , int_fast32_t > * mat = new Eigen : : SparseMatrix < T , Eigen : : RowMajor , int_fast32_t > ( eigenRows , eigenRows ) ;
/ / Check whether it is safe to export this matrix .
if ( ! isReadReady ( ) ) {
triggerErrorState ( ) ;
pantheios : : log_ERROR ( " StaticSparseMatrix::toEigenSparseMatrix: Throwing invalid state for internal state not ReadReady (is " , pantheios : : integer ( internal_status ) , " ). " ) ;
throw mrmc : : exceptions : : invalid_state ( " StaticSparseMatrix::toEigenSparseMatrix: Invalid state for internal state not ReadReady. " ) ;
} else {
/ / Create a
int_fast32_t eigenRows = static_cast < int_fast32_t > ( row_count ) ;
Eigen : : SparseMatrix < T , Eigen : : RowMajor , int_fast32_t > * mat = new Eigen : : SparseMatrix < T , Eigen : : RowMajor , int_fast32_t > ( eigenRows , eigenRows ) ;
/ / There are two ways of converting this matrix to Eigen ' s format .
/ / 1. Compute a list of triplets ( row , column , value ) for all
/ / non - zero elements and pass it to Eigen to create a sparse matrix .
/ / 2. Tell Eigen to reserve the average number of non - zero elements
/ / per row in advance and insert the values via a call to Eigen ' s
/ / insert method then . As the given reservation number is only an
/ / estimate , the actual number may be different and Eigen may have
/ / to shift a lot .
/ / In most cases , the second alternative is faster ( about 1 / 2 of the
/ / first , but if there are " heavy " rows that are several times larger
/ / than an average row , the other solution might be faster .
/ / The desired conversion method may be set by an appropriate define .
# ifdef MRMC_USE_TRIPLETCONVERT
/ / FIXME : Wouldn ' t it be more efficient to add the elements in
/ / order including the diagonal elements ? Otherwise , Eigen has to
/ / perform some sorting .
/ / Prepare the triplet storage .
typedef Eigen : : Triplet < T > IntTriplet ;
std : : vector < IntTriplet > tripletList ;
tripletList . reserve ( non_zero_entry_count + row_count ) ;
/ / First , iterate over all elements that are not on the diagonal
/ / and add the corresponding triplet .
uint_fast64_t row_start ;
uint_fast64_t row_end ;
for ( uint_fast64_t row = 1 ; row < = row_count ; + + row ) {
row_start = row_indications [ row - 1 ] ;
row_end = row_indications [ row ] ;
for ( uint_fast64_t row = 0 ; row < = row_count ; + + row ) {
row_start = row_indications [ row ] ;
row_end = row_indications [ row + 1 ] ;
while ( row_start < row_end ) {
tripletList . push_back ( IntTriplet ( row - 1 , column_indications [ row_start ] - 1 , value_storage [ row_start ] ) ) ;
tripletList . push_back ( IntTriplet ( row , column_indications [ row_start ] , value_storage [ row_start ] ) ) ;
+ + row_start ;
}
}
for ( uint_fast64_t i = 1 ; i < = row_count ; + + i ) {
tripletList . push_back ( IntTriplet ( i - 1 , i - 1 , diagonal_storage [ i ] ) ) ;
/ / Then add the elements on the diagonal .
for ( uint_fast64_t i = 0 ; i < = row_count ; + + i ) {
tripletList . push_back ( IntTriplet ( i , i , diagonal_storage [ i ] ) ) ;
}
/ / Let Eigen create a matrix from the given list of triplets .
mat - > setFromTriplets ( tripletList . begin ( ) , tripletList . end ( ) ) ;
# else / / NOT MRMC_USE_TRIPLETCONVERT
/ / In most cases , this is faster ( about 1 / 2 of the above ) . But if there are " heavy " rows that are several times larger than an average row , the other solution might be faster .
mat - > reserve ( Eigen : : VectorXi : : Constant ( static_cast < int_fast32_t > ( mat - > outerSize ( ) ) , static_cast < int_fast32_t > ( ( non_zero_entry_count + row_count ) / mat - > outerSize ( ) ) ) ) ;
/ / Reserve the average number of non - zero elements per row for each
/ / row .
mat - > reserve ( Eigen : : VectorXi : : Constant ( eigenRows , static_cast < int_fast32_t > ( ( non_zero_entry_count + row_count ) / eigenRows ) ) ) ;
/ / Iterate over the all non - zero elements in this matrix and add
/ / them to the matrix individually .
uint_fast64_t row_start ;
uint_fast64_t row_end ;
for ( uint_fast64_t row = 1 ; row < = row_count ; + + row ) {
row_start = row_indications [ row - 1 ] ;
row_end = row_indications [ row ] ;
for ( uint_fast64_t row = 0 ; row < = row_count ; + + row ) {
row_start = row_indications [ row ] ;
row_end = row_indications [ row + 1 ] ;
/ / insert the diagonal entry
mat - > insert ( row - 1 , row - 1 ) = diagonal_storage [ row ] ;
/ / Insert the element on the diagonal .
mat - > insert ( row , row ) = diagonal_storage [ row ] ;
/ / Insert the elements that are not on the diagonal
while ( row_start < row_end ) {
/ / tripletList . push_back ( IntTriplet ( row - 1 , column_indications [ row_start ] - 1 , value_storage [ row_start ] ) ) ;
mat - > insert ( row - 1 , column_indications [ row_start ] - 1 ) = value_storage [ row_start ] ;
mat - > insert ( row , column_indications [ row_start ] ) = value_storage [ row_start ] ;
+ + row_start ;
}
}
# endif / / MRMC_USE_TRIPLETCONVERT
/ / Make the matrix compressed , i . e . remove remaining zero - entries .
mat - > makeCompressed ( ) ;
}
return mat ;
}
/ / ! Returns the exact count of explicit entries in the sparse matrix .
/ / This point can never be reached as both if - branches end in a return
/ / statement .
return nullptr ;
}
/*!
Retuns the exact count of explicit entries in the sparse matrix . While it is called " nonZero " count , the fields may of course be 0 ( the 0 - value of T ) .
@ returns explicit entry count in the matrix
* Returns the number of non - zero entries that are not on the diagonal .
* @ returns The number of non - zero entries that are not on the diagonal .
*/
uint_fast64_t getNonZeroEntryCount ( ) const {
return non_zero_entry_count ;
}
/ / ! Converts a state into a final state .
/*!
This function allows for states to be made final . This means that all entries in row " state " will be changed to 0 but for the one on the diagonal , which will be set to 1 to create a loop on itself .
@ param state The number of the state to be converted . Must be in 1 < = state < = rows .
@ returns Whether the conversion was successful .
* This function makes the given state absorbing . This means that all
* entries in its row will be changed to 0 and the value 1 will be written
* to the element on the diagonal .
* @ param state The state to be made absorbing .
* @ returns True iff the operation was successful .
*/
bool makeStateFinal ( const uint_fast64_t state ) {
if ( ( state > row_count ) | | ( state = = 0 ) ) {
pantheios : : log_ERROR ( " StaticSparseMatrix::makeStateFinal: state not in 1 .. rows (is " , pantheios : : integer ( state ) , " , max is " , pantheios : : integer ( row_count ) , " ). " ) ;
throw mrmc : : exceptions : : out_of_range ( " StaticSparseMatrix::makeStateFinal: state not in 1 .. rows " ) ;
bool makeStateAbsorbing ( const uint_fast64_t state ) {
/ / Check whether the accessed state exists .
if ( state > row_count ) {
pantheios : : log_ERROR ( " StaticSparseMatrix::makeStateFinal: state not in 0 .. rows (is " , pantheios : : integer ( state ) , " , max is " , pantheios : : integer ( row_count ) , " ). " ) ;
throw mrmc : : exceptions : : out_of_range ( " StaticSparseMatrix::makeStateFinal: state not in 0 .. rows " ) ;
return false ;
}
uint_fast64_t row_start = row_indications [ state - 1 ] ;
uint_fast64_t row_end = row_indications [ state ] ;
/ / Iterate over the elements in the row that are not on the diagonal
/ / and set them to zero .
uint_fast64_t row_start = row_indications [ state ] ;
uint_fast64_t row_end = row_indications [ state + 1 ] ;
while ( row_start < row_end ) {
value_storage [ row_start ] = mrmc : : misc : : constGetZero ( value_storage ) ;
row_start + + ;
}
/ / Set the element on the diagonal to one .
diagonal_storage [ state ] = mrmc : : misc : : constGetOne ( diagonal_storage ) ;
return true ;
}
@ -473,71 +608,128 @@ class StaticSparseMatrix {
}
private :
uint_fast64_t current_size ;
/*!
* The number of rows of the matrix .
*/
uint_fast64_t row_count ;
/*!
* The number of non - zero elements that are not on the diagonal .
*/
uint_fast64_t non_zero_entry_count ;
uint_fast64_t last_row ;
/*! Array containing all non-zero values, apart from the diagonal entries */
/*!
* Stores all non - zero values that are not on the diagonal .
*/
T * value_storage ;
/*! Array containing all diagonal values */
/*!
* Stores all elements on the diagonal , even the ones that are zero .
*/
T * diagonal_storage ;
/*! Array containing the column number of the corresponding value_storage entry */
/*!
* Stores the column for each non - zero element that is not on the diagonal .
*/
uint_fast64_t * column_indications ;
/*! Array containing the row boundaries of valueStorage */
/*!
* Array containing the boundaries ( indices ) in the value_storage array
* for each row . All elements of value_storage with indices between the
* i - th and the ( i + 1 ) - st element of this array belong to row i .
*/
uint_fast64_t * row_indications ;
/*! Internal status enum, 0 for constructed, 1 for initialized and 2 for finalized, -1 on errors */
/*!
* The internal status of the matrix .
*/
MatrixStatus internal_status ;
/*! Sets the internal status to Error */
/*!
* Stores the current number of non - zero elements that have been added to
* the matrix . Used for correctly inserting elements in the matrix .
*/
uint_fast64_t current_size ;
/*!
* Stores the row in which the last element was inserted . Used for correctly
* inserting elements in the matrix .
*/
uint_fast64_t last_row ;
/*!
* Sets the internal status to signal an error .
*/
void triggerErrorState ( ) {
setState ( MatrixStatus : : Error ) ;
}
/*!
Sets the internal status to new_state iff the current state is not the Error state .
@ param new_state the new state to be switched to
* Sets the internal status to the given state if the current state is not
* the error state .
* @ param new_state The new state to be switched to .
*/
void setState ( const MatrixStatus new_state ) {
internal_status = ( internal_status = = MatrixStatus : : Error ) ? internal_status : new_state ;
}
/*!
Prepares the internal CSR storage .
Requires non_zero_entry_count and row_count to be set .
@ return t rue on success , false otherwise ( allocation failed ) .
* Prepares the internal CSR storage . For this , it requires
* non_zero_entry_count and row_count to be set correctly .
* @ return T rue on success , false otherwise ( allocation failed ) .
*/
bool prepareInternalStorage ( ) {
/ / Set up the arrays for the elements that are not on the diagonal .
value_storage = new ( std : : nothrow ) T [ non_zero_entry_count ] ( ) ;
column_indications = new ( std : : nothrow ) uint_fast64_t [ non_zero_entry_count ] ( ) ;
/ / Set up the row_indications array and reserve one element more than
/ / there are rows in order to put a sentinel element at the end ,
/ / which eases iteration process .
row_indications = new ( std : : nothrow ) uint_fast64_t [ row_count + 1 ] ( ) ;
/ / row_count + 1 so that access with 1 - based indices can be direct without the overhead of a - 1 each time
diagonal_storage = new ( std : : nothrow ) T [ row_count + 1 ] ( ) ;
/ / Set up the array for the elements on the diagonal .
diagonal_storage = new ( std : : nothrow ) T [ row_count ] ( ) ;
return ( ( value_storage ! = NULL ) & & ( column_indications ! = NULL ) & & ( row_indications ! = NULL ) & & ( diagonal_storage ! = NULL ) ) ;
/ / Return whether all the allocations could be made without error .
return ( ( value_storage ! = NULL ) & & ( column_indications ! = NULL )
& & ( row_indications ! = NULL ) & & ( diagonal_storage ! = NULL ) ) ;
}
/ / !
/*!
* Helper function to determine whether the given Eigen matrix is in RowMajor
* format . Always returns true , but is overloaded , so the compiler will
* only call it in case the Eigen matrix is in RowMajor format .
* @ return True .
*/
template < typename _Scalar , typename _Index >
bool isEigenRowMajor ( Eigen : : SparseMatrix < _Scalar , Eigen : : RowMajor , _Index > ) {
return true ;
}
/*!
* Helper function to determine whether the given Eigen matrix is in RowMajor
* format . Always returns false , but is overloaded , so the compiler will
* only call it in case the Eigen matrix is in ColMajor format .
* @ return False .
*/
template < typename _Scalar , typename _Index >
bool isEigenRowMajor ( Eigen : : SparseMatrix < _Scalar , Eigen : : ColMajor , _Index > ) {
bool isEigenRowMajor (
Eigen : : SparseMatrix < _Scalar , Eigen : : ColMajor , _Index > ) {
return false ;
}
/*!
* Helper function to determine the number of non - zero elements that are
* not on the diagonal of the given Eigen matrix .
* @ param eigen_sparse_matrix The Eigen matrix to analyze .
* @ return The number of non - zero elements that are not on the diagonal of
* the given Eigen matrix .
*/
template < typename _Scalar , int _Options , typename _Index >
_Index getEigenSparseMatrixCorrectNonZeroEntryCount ( const Eigen : : SparseMatrix < _Scalar , _Options , _Index > & eigen_sparse_matrix ) {
_Index getEigenSparseMatrixCorrectNonZeroEntryCount (
const Eigen : : SparseMatrix < _Scalar , _Options , _Index > & eigen_sparse_matrix ) {
const _Index * indexPtr = eigen_sparse_matrix . innerIndexPtr ( ) ;
const _Index * outerPtr = eigen_sparse_matrix . outerIndexPtr ( ) ;
@ -546,8 +738,9 @@ class StaticSparseMatrix {
uint_fast64_t diag_non_zeros = 0 ;
/ / for RowMajor , row is the current Row and col the column
/ / for ColMajor , row is the current Col and col the row
/ / For RowMajor , row is the current Row and col the column and for
/ / ColMajor , row is the current column and col the row , but this is
/ / not important as we are only looking for elements on the diagonal .
_Index innerStart = 0 ;
_Index innerEnd = 0 ;
_Index innerMid = 0 ;
@ -555,7 +748,7 @@ class StaticSparseMatrix {
innerStart = outerPtr [ row ] ;
innerEnd = outerPtr [ row + 1 ] - 1 ;
/ / Now with super fancy binary search , deferred equality detection
/ / Now use binary search ( but defer equality detection ) .
while ( innerStart < innerEnd ) {
innerMid = innerStart + ( ( innerEnd - innerStart ) / 2 ) ;
@ -566,8 +759,8 @@ class StaticSparseMatrix {
}
}
/ / Check whether we have found an element on the diagonal .
if ( ( innerStart = = innerEnd ) & & ( indexPtr [ innerStart ] = = row ) ) {
/ / found a diagonal entry
+ + diag_non_zeros ;
}
}