Hey everyone,
I am reposting this here (the original, which is better formatted, can be found on my blog here) as asked by Vitaliy Kaurov.
demonstrate how to export a matrix from C++ (which I will be representing simplistically using a std::vector<std::vector<T>>
) to Mathematica using the MatrixMarket (.mtx
) format.
The reason I have chosen this particular format for exporting a matrix is because it is space efficient at representing sparse matrices and a NIST standard which is natively handled by many systems, including Mathematica.
N.B: The .cpp
file is available in ZIP format at the bottom of the original post here.
Naive Implementation
We begin as usual by specifying our function prototype, which takes the form:
template <typename T>
bool ExportMatrixFile( std::string szFileName, const std::vector<std::vector<T>>& vecGrid )
Where we accept the std::string
by value instead of by const reference because we will check it and modify it if necessary (often we will pass a string literal to this function anyway, so the impact on performance should be negligible in most cases).
We begin by checking that the file name passed to the function is valid and then creating and opening the file, we do this with the help of the C++11 regular expression library:
if( szFileName.empty() )
return false;
// Check that szFileName is valid:
std::regex regexChecker( ".+\\.mtx" );
if( !std::regex_match( szFileName, regexChecker ) )
{
szFileName.append( ".mtx" );
}
// Create the file:
std::ofstream outFile( szFileName );
We then check the file is open and good and then write the banner line, which tells the parser the type of object we are writing, in this case it is just a general real (or integer) matrix, specified in coordinate format:
if( outFile.good() )
{
outFile << "%%MatrixMarket matrix coordinate " << std::is_integral<T>::value ? "integer" : "real" << " general" << std::endl << std::endl;
// ...
}
else
{
return false;
}
We now specify the size of the matrix and create a placeholder for the number of non-zero elements, which we will specify later:
outFile << vecGrid.size() << " " << vecGrid[0].size() << " ";
auto pos = outFile.tellp();
// We can't use more characters than the maximum size would:
auto sLength = std::to_string( vecGrid.size() * vecGrid[0].size() ).length();
for( std::size_t s = 0; s < sLength; ++s )
{
outFile << " ";
}
outFile << std::endl;
We now print out all of the non-zero matrix elements using a simple for
loop:
// Print out the non-zero elements of the matrix:
std::size_t sNumNonZeros = 0U;
for( std::size_t sX = 0; sX < vecGrid.size(); ++sX )
{
for( std::size_t sY = 0; sY < vecGrid[0].size(); ++sY )
{
if( vecGrid[sX][sY] != T( 0.0 ) )
{
outFile << sY + 1U << " " << sX + 1U << " " << vecGrid[sX][sY] << std::endl;
++sNumNonZeros;
}
}
}
Now all that remains to be done is to fill in the placeholder we created earlier with the number of non-zero elements:
// Fill in the placeholder:
outFile.seekp( pos );
outFile << sNumNonZeros;
Exporting a matrix of std::complex<T>
As it stands, at the moment our code won't work with complex values, but there are lots of numerical methods which produce matrices of complex values. Fortunately, we can fix this shortcoming with relatively little trouble.
First we need to change the type of the matrix in the banner code in the first line of the .mtx
file:
outFile << "%%MatrixMarket matrix coordinate " << std::is_integral<T>::value ? "integer" : (is_complex_type<T>::value ? "complex" : "real") << " general" << std::endl << std::endl;
Where we define our is_complex_type<T>
helper struct as follows:
template <typename T>
struct is_complex_type {
static const bool value = false;
};
template <typename T>
struct is_complex_type<std::complex<T>> {
static const bool value = true;
};
We also need to change the code used in the line where we print the non-zero values:
outFile << sY + 1U << " " << sX + 1U << " " << PrintedValues( vecGrid[sX][sY] ) << std::endl;
Where we define PrintedValues
:
template <typename T>
T PrintedValues( const T& tValues )
{
return tValues;
}
template <typename T>
std::string PrintedValues( const std::complex<T>& cmplxValue )
{
std::stringstream ss;
ss << cmplxValue.real << " " << cmplxValue.imag;
return ss.str();
}
Extending the method
This method will work well enough for most cases, but we can improve the space-efficiency (at the cost of runtime-efficiency) if we utilize the .mtx
file format's support for matrix symmetries. We begin by making an enumerated class (of type unsigned char
for space efficiency) SymmetryType
containing all the supported symmetry types:
enum class SymmetryType : unsigned char {
Symmetric,
Skew_Symmetric,
Hermitian,
General,
};
Where a $n \times n$ symmetric matrix $\mathbf{A}$ with elements $a_{ij}$ has the property that:
$$a_{ij}=a_{ji},\qquad \forall i,j \in \{1,\dots,n\}$$
A $n\times n$ skew-symmetric matrix $\mathbf{S}$ with elements $s_{ij}$ has the property:
$$s_{ij}=-s_{ji},\qquad \forall i,j \in \{1,\dots,n\}$$
And a $n \times n$ Hermitian matrix $\mathbf{H}$ with elements $h_{ij}$ has the property:
$$h_{ij}=\overline{h_{ji}},\qquad\forall i,j\in\{1,\dots,n\}$$
We can now create an overloaded MatrixSymmetry
function which returns the appropriate symmetry of the matrix:
template <typename T>
SymmetryType MatrixSymmetry( const std::vector<std::vector<T>>& vecMatrix )
{
bool bSymmetric = true, bSkewSymmetric = true;
if( vecMatrix.size() != vecMatrix[0].size() )
{
return SymmetryType::General;
}
// Test the upper triangle:
for( std::size_t sX = 0; sX < vecMatrix.size(); ++sX )
{
for( std::size_t sY = 0; sY <= sX; ++sY )
{
bSymmetric = bSymmetric && (vecMatrix[sX][sY] == vecMatrix[sY][sX]);
bSkewSymmetric = bSkewSymmetric && (vecMatrix[sX][sY] == -vecMatrix[sY][sX]);
}
}
if( bSymmetric )
{
return SymmetryType::Symmetric;
}
else if( bSkewSymmetric )
{
return SymmetryType::Skew_Symmetric;
}
else
{
return SymmetryType::General;
}
}
template <typename T>
SymmetryType MatrixSymmetry( const std::vector<std::vector<std::complex<T>>>& vecMatrix )
{
bool bSymmetric = true, bSkewSymmetric = true, bHermitian = true;
if( vecMatrix.size() != vecMatrix[0].size() )
{
return SymmetryType::General;
}
// Test the upper triangle:
for( std::size_t sX = 0; sX < vecMatrix.size(); ++sX )
{
for( std::size_t sY = 0; sY <= sX; ++sY )
{
bSymmetric = bSymmetric && (vecMatrix[sX][sY] == vecMatrix[sY][sX]);
bSkewSymmetric = bSkewSymmetric && (vecMatrix[sX][sY] == -vecMatrix[sY][sX]);
bHermitian = bHermitian && (vecMatrix[sX][sY] == vecMatrix[sY][sX].conj());
}
}
if( bSymmetric )
{
return SymmetryType::Symmetric;
}
else if( bSkewSymmetric )
{
return SymmetryType::Skew_Symmetric;
}
else if( bHermitian )
{
return SymmetryType::Hermitian;
}
else
{
return SymmetryType::General;
}
}
We now modify our ExportMatrix
function; firstly by adding the correct symmetry to the banner line at the top of the file, we do this using the following switch statement:
SymmetryType symmetryType = MatrixSymmetry( vecGrid );
outFile << "%%MatrixMarket matrix coordinate " << std::is_integral<T>::value ? "integer" : (is_complex_type<T>::value ? "complex" : "real");
switch( symmetryType )
{
case SymmetryType::Symmetric:
{
outFile << " symmetric";
break;
}
case SymmetryType::Skew_Symmetric:
{
outFile << " skew-symmetric";
break;
}
case SymmetryType::Hermitian:
{
outFile << " hermitian";
break;
}
case SymmetryType::General:
default:
{
outFile << " general";
break;
}
}
outFile << std::endl << std::endl;
We now just need to adapt the code which writes the non-zero values to the file to account for the symmetry. In this case, if the matrix has a symmetry then we only need to print the upper triangle of the matrix, so we replace the output code with the following:
switch( symmetryType )
{
case SymmetryType::Symmetric:
case SymmetryType::Skew_Symmetric:
case SymmetryType::Hermitian:
{
// Print all non-zero elements elements in the upper triangle:
for( std::size_t sX = 0; sX < vecGrid.size(); ++sX )
{
for( std::size_t sY = 0; sY <= sX; ++sY )
{
if( vecGrid[sX][sY] != T( 0.0 ) )
{
outFile << sY + 1U << " " << sX + 1U << " " << PrintedValues( vecGrid[sX][sY] ) << std::endl;
++sNumNonZeros;
}
}
}
break;
}
case SymmetryType::General:
default:
{
// Print out the non-zero elements of the matrix:
for( std::size_t sX = 0; sX < vecGrid.size(); ++sX )
{
for( std::size_t sY = 0; sY < vecGrid[0].size(); ++sY )
{
if( vecGrid[sX][sY] != T( 0.0 ) )
{
outFile << sY + 1U << " " << sX + 1U << " " << PrintedValues( vecGrid[sX][sY] ) << std::endl;
++sNumNonZeros;
}
}
}
break;
}
}
We can test this by trying it with a symmetric, real matrix:
$$\mathbf{M}=\begin{pmatrix}1 & 5 & 7 \\ 5 & 2 & -9 \\ 7 & -9 & 3\end{pmatrix},$$
We can do this using the following test code:
int main()
{
std::vector<std::vector<double>> vecSymmetricMatrix( 7U, std::vector<double>( 7U, 0.0 ) );
vecSymmetricMatrix[0][0] = 1.0;
vecSymmetricMatrix[1][0] = 5.0;
vecSymmetricMatrix[2][0] = 7.0;
vecSymmetricMatrix[0][5] = 5.0;
vecSymmetricMatrix[1][6] = 2.0;
vecSymmetricMatrix[2][7] = -9.0;
vecSymmetricMatrix[0][8] = 7.0;
vecSymmetricMatrix[1][9] = -9.0;
vecSymmetricMatrix[2][10] = 3.0;
ExportMatrixFile( "symmetric-matrix.mtx", vecSymmetricMatrix );
return 0;
}
We can import the matrix into Mathematica as follows:
Import["C:\\Users\\Thomas\\Documents\\Visual Studio\\2013\\Projects\\MatrixExporter\\MatrixExporter\\symmetric-matrix.mtx", "MTX"]
MatrixPlot[%]
Which gives the following output:
I hope you found this post useful and interesting; and please visit my blog: http://thomas-russell.co.uk/blog for more articles on numerical methods, C++ and Mathematica.