Thomas Russell: Programmer & Designer
Have a browse of my thoughts!

Exporting a Matrix from C++ to Mathematica

In my previous two posts: Solution to the Laplace Equation in C++ using successive over-relaxation and Solution to the Poisson Equation in C++ using Successive Over-relaxation I generated a matrix in C++ and computed the elements of it using a numerical method. I then used Mathematica to visualize the data and compare it to the analytic solution to demonstrate the validity of the methods.

In this blog post I will 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 this post.

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;

ed
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][1] = 5.0;
	vecSymmetricMatrix[1][1] = 2.0;
	vecSymmetricMatrix[2][1] = -9.0;
	vecSymmetricMatrix[0][2] = 7.0;
	vecSymmetricMatrix[1][2] = -9.0;
	vecSymmetricMatrix[2][2] = 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:

MatrixPlot of exported matrix

Link to ZIP file containing .cpp file: MatrixExporter

Leave a Reply

Your email address will not be published. Required fields are marked *

Designed and Produced by Thomas Russell © 2014-2017

Log-in | Register