/*
 *  Portable Agile C++ Classes (PACC)
 *  Copyright (C) 2001-2004 by Marc Parizeau
 *  http://manitou.gel.ulaval.ca/~parizeau/PACC
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *  Contact:
 *  Laboratoire de Vision et Systemes Numeriques
 *  Departement de genie electrique et de genie informatique
 *  Universite Laval, Quebec, Canada, G1K 7P4
 *  http://vision.gel.ulaval.ca
 *
 */

/*!
 * \file   PACC/Math/Matrix.cpp
 * \brief  Method definitions for class Matrix.
 * \author Marc Parizeau and Christian Gagné, Laboratoire de vision et 
 systèmes numériques, Université Laval
 * $Revision: 1.9.2.1 $
 * $Date: 2007/09/10 18:24:08 $
 */

#include "Math/Matrix.hpp"
#include "Math/Vector.hpp"
#include "Util/StringFunc.hpp"
#include <stdexcept>
#include <iomanip>
#include <cmath>

using namespace std;
using namespace PACC;

/*! 
This method also returns a reference to the result.
*/
Matrix& Matrix::add(Matrix& outMatrix, double inScalar) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "add() invalid matrix!");
	outMatrix.setRowsCols(mRows, mCols);	
	for(unsigned int i = 0; i < size(); ++i) outMatrix[i] = (*this)[i] + inScalar;
	return outMatrix;
}

/*! 
This method also returns a reference to the result.
*/
Matrix& Matrix::add(Matrix& outMatrix, const Matrix& inMatrix) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "add() invalid matrix!");
	PACC_AssertM(mRows == inMatrix.mRows && mCols == inMatrix.mCols, "add() matrix mismatch!");	
	outMatrix.setRowsCols(mRows, mCols);	
	for(unsigned int i = 0; i < size(); ++i) outMatrix[i] = (*this)[i] + inMatrix[i];
	return outMatrix;
}

/*!
*/
void Matrix::computeBackSubLU(const vector<unsigned int>& inIndexes, Matrix& ioMatrix) const
{
	unsigned int lII = UINT_MAX;
	for(unsigned int i = 0; i < mRows; ++i) {
		double lSum = ioMatrix(inIndexes[i], 0);
		ioMatrix(inIndexes[i], 0) = ioMatrix(i, 0);
		if(lII != UINT_MAX) {
			for(unsigned int j = lII; j < i; ++j) lSum -= (*this)(i,j) * ioMatrix(j, 0);
		} else if(lSum != 0.0) lII = i;
		ioMatrix(i, 0) = lSum;
	}
	for(unsigned int i = 0; i < mRows; ++i) {
		const unsigned int lR = mRows-i-1;
		double lSum = ioMatrix(lR,0);
		for(unsigned int j = lR+1; j < mCols; ++j) lSum -= (*this)(lR,j) * ioMatrix(j,0);
		ioMatrix(lR,0) = lSum / (*this)(lR,lR);
	}
}

/*!
*/
double Matrix::computeDeterminant(void) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "computeDeterminant() invalid matrix!");
	PACC_AssertM(mRows == mCols, "computeDeterminant() matrix not square!");
	Matrix lTmp = *this;
	vector<unsigned int> lIndexes(mRows);
	int lD;
	lTmp.decomposeLU(lIndexes, lD);
	double lResult = lD;
	for(unsigned int i = 0; i < mRows; ++i) lResult *= lTmp(i,i);
	return lResult;
}

/*!
*/
void Matrix::computeEigens(Vector& outValues, Matrix& outVectors) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "computeEigens() invalid matrix!");
	PACC_AssertM(mRows == mCols, "computeEigens() matrix not square!");
	outValues.resize(mRows);
	outVectors.resize(mRows, mCols);
	
	// Computer eigenvectors/eigenvalues using Triagonal QL method
	Vector lE(mRows);
	tred2(outValues, lE, outVectors);
	tql2(outValues, lE, outVectors);
	
	// Sort by eigenvalues.
	for(unsigned int j = 0; j < outValues.size(); ++j) {
		double lMax=outValues[j];
		unsigned int lMaxArg=j;
		for(unsigned int l = j+1; l<outValues.size(); ++l) {
			if(outValues[l] > lMax) {
				lMax=outValues[l];
				lMaxArg=l;
			}
		}
		if(lMaxArg != j) {
			for(unsigned int r = 0; r < outVectors.mRows; ++r) {
				double lTmp = outVectors(r,j);
				outVectors(r,j) = outVectors(r,lMaxArg);
				outVectors(r,lMaxArg) = lTmp;
			}
			double lTmp = outValues[j];
			outValues[j] = outValues[lMaxArg];
			outValues[lMaxArg] = lTmp;
		}
	}
}

/*!
*/
void Matrix::decomposeLU(vector<unsigned int>& outIndexes, int& outD)
{
	outD = 1;
	vector<double> lScales;
	scaleLU(lScales);
	for(unsigned int j = 0; j < mCols; ++j) {
		for(unsigned int i = 0; i < j; ++i) {
			double lSum = (*this)(i, j);
			for(unsigned int k = 0; k < i; ++k) lSum -= (*this)(i,k) * (*this)(k,j);
			(*this)(i, j) = lSum;
		}
		double lMax = 0;
		unsigned int l = j;
		for(unsigned int i = j; i < mRows; ++i) {
			double lSum = (*this)(i,j);
			for(unsigned int k = 0; k < j; ++k) lSum -= (*this)(i,k) * (*this)(k,j);
			(*this)(i, j) = lSum;
			double lTmp = lScales[i] * fabs(lSum);
			if(lTmp >= lMax) {
				l = i;
				lMax = lTmp;
			}
		}
		if(j != l) {
			for(unsigned int k = 0; k < (*this).mCols; ++k) {
				double lTmp = (*this)(l,k);
				(*this)(l,k) = (*this)(j,k);
				(*this)(j,k) = lTmp;
			}
			outD = -outD;
			lScales[l] = lScales[j];
		}
		outIndexes[j] = l;
		if((*this)(j,j) == 0.0) (*this)(j,j) = 1e-20;
		if(j != (mCols-1)) {
			double lDummy = 1.0 / (*this)(j,j);
			for(unsigned int i = j+1; i < mRows; ++i) (*this)(i,j) *= lDummy;
		}
	}
}

/*!
This method also returns a reference to the result.
*/
Matrix& Matrix::extract(Matrix& outMatrix, unsigned int inRow1, unsigned int inRow2, unsigned int inCol1, unsigned int inCol2) const
{
	PACC_AssertM(inRow1 <= inRow2 && inCol1 <= inCol2 && inRow2 < mRows && inCol2 < mCols, "extract() invalid indexes!");
	if(&outMatrix != this) {
		// output matrix is not self assigning
		outMatrix.setRowsCols(inRow2-inRow1+1, inCol2-inCol1+1);
		for(unsigned int i = inRow1; i <= inRow2; ++i) {
			for(unsigned int j = inCol1; j <= inCol2; ++j) {
				outMatrix(i-inRow1,j-inCol1) = (*this)(i,j);
			}
		}
	} else {
		// use temporary matrix to self assign
		Matrix lMatrix(*this);
		outMatrix.setRowsCols(inRow2-inRow1+1, inCol2-inCol1+1);
		for(unsigned int i = inRow1; i <= inRow2; ++i) {
			for(unsigned int j = inCol1; j <= inCol2; ++j) {
				outMatrix(i-inRow1,j-inCol1) = lMatrix(i,j);
			}
		}
	}
	return outMatrix;
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::extractColumn(Matrix& outMatrix, unsigned int inCol) const
{
	return extract(outMatrix, 0, mRows-1, inCol, inCol);
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::extractRow(Matrix& outMatrix, unsigned int inRow) const
{
	return extract(outMatrix, inRow, inRow, 0, mCols-1);
}

/*!
 */
double Matrix::hypot(double a, double b) const
{
	double r;
	if(abs(a) > abs(b)) {
		r = b/a;
		r = abs(a)*sqrt(1+r*r);
	}
	else if(b != 0) {
		r = a/b;
		r = abs(b)*sqrt(1+r*r);
	}
	else {
		r = 0.0;
	}
	return r;
}

/*!
*/
Matrix Matrix::invert(void) const 
{
	Matrix lMatrix; 
	return invert(lMatrix);
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::invert(Matrix& outMatrix) const
{
	PACC_AssertM(mRows == mCols, "invert() matrix not square!");
	Matrix lTmp = *this;
	vector<unsigned int> lIndexes(mRows);
	int lD;
	lTmp.decomposeLU(lIndexes, lD);
	outMatrix.setIdentity(mRows);
	Matrix lB(mRows, 1);
	for(unsigned int j = 0; j < mCols; ++j) {
		for(unsigned int i = 0; i < mRows; ++i) lB(i,0) = outMatrix(i,j);
		lTmp.computeBackSubLU(lIndexes, lB);
		for(unsigned int i = 0; i < mRows; ++i) outMatrix(i,j) = lB(i,0);
	}
	return outMatrix;
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::multiply(Matrix& outMatrix, double inScalar) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "multiply() invalid matrix!");
	outMatrix.setRowsCols(mRows, mCols);
	for(unsigned int i = 0; i < size(); ++i) outMatrix[i] = (*this)[i] * inScalar;
	return outMatrix;
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::multiply(Matrix& outMatrix, const Matrix& inMatrix) const
{
	PACC_AssertM(mCols == inMatrix.mRows, "multiply() matrix mismatch!");
	if(&outMatrix != this && &outMatrix != &inMatrix) {
		// output matrix is neither left or right matrix (no self assigment)
		outMatrix.setRowsCols(mRows, inMatrix.mCols);
		for(unsigned int i = 0; i < outMatrix.mRows; ++i) {
			for(unsigned int j = 0; j < outMatrix.mCols; ++j) {
				outMatrix(i,j) = 0;
				for(unsigned int k = 0; k < mCols; ++k) {
					outMatrix(i,j) += (*this)(i,k) * inMatrix(k,j);
				}
			}
		}
	} else if(&outMatrix == this && &outMatrix != &inMatrix) {
		// use temporary matrix to self assign with left matrix
		Matrix lMatrix(*this);
		outMatrix.setRowsCols(mRows, inMatrix.mCols);
		for(unsigned int i = 0; i < outMatrix.mRows; ++i) {
			for(unsigned int j = 0; j < outMatrix.mCols; ++j) {
				outMatrix(i,j) = 0;
				for(unsigned int k = 0; k < mCols; ++k) {
					outMatrix(i,j) += lMatrix(i,k) * inMatrix(k,j);
				}
			}
		}
	} else if(&outMatrix != this && &outMatrix == &inMatrix) {
		// use temporary matrix to self assign with right matrix
		Matrix lMatrix(inMatrix);
		outMatrix.setRowsCols(mRows, inMatrix.mCols);
		for(unsigned int i = 0; i < outMatrix.mRows; ++i) {
			for(unsigned int j = 0; j < outMatrix.mCols; ++j) {
				outMatrix(i,j) = 0;
				for(unsigned int k = 0; k < mCols; ++k) {
					outMatrix(i,j) += (*this)(i,k) * lMatrix(k,j);
				}
			}
		}
	} else {
		// use temporary matrix to self assign with both left and right matrices
		Matrix lMatrix(*this);
		outMatrix.setRowsCols(mRows, inMatrix.mCols);
		for(unsigned int i = 0; i < outMatrix.mRows; ++i) {
			for(unsigned int j = 0; j < outMatrix.mCols; ++j) {
				outMatrix(i,j) = 0;
				for(unsigned int k = 0; k < mCols; ++k) {
					outMatrix(i,j) += lMatrix(i,k) * lMatrix(k,j);
				}
			}
		}
	}
	return outMatrix;
}

/*!
Matrix elements must be enumerated in row order and delimited by either commas 
(','), semi-columns (';'), or white space. The recommended style is to seperate 
elements with comas, and rows with semi-columns. For example:
\verbatim
<Matrix name="My Matrix" rows="3" cols="4">1,2,3,4;5,6,7,8;9,10,11,12</Matrix>
\endverbatim
The number of elements must match the product of the "rows" and "cols" attributes.
*/
string Matrix::read(const XML::Iterator& inNode)
{
	if(!inNode) throw runtime_error("Matrix::read() nothing to read!");
	clear();
	for(XML::Iterator lChild = inNode->getFirstChild(); lChild; ++lChild) {
		if(lChild->getType() == XML::eString) {
			istringstream lStream(lChild->getValue());
			Tokenizer lTokenizer(lStream);
			lTokenizer.setDelimiters(" \n\r\t,;", "");
			string lToken;
			while(lTokenizer.getNextToken(lToken)) push_back(String::convertToFloat(lToken));
		}
	}
	mRows = String::convertToInteger(inNode->getAttribute("rows"));
	mCols = String::convertToInteger(inNode->getAttribute("cols"));
	if(vector<double>::size() != mRows*mCols) {
		throwError("Matrix::read() number of elements does not match the rows x cols attributes", inNode);
	}
	string lName = inNode->getAttribute("name");
	if(lName != "") mName = lName;
	return lName;
}

/*!
*/
void Matrix::resize(unsigned int inRows, unsigned int inCols)
{
	Matrix lMat(*this);
	setRowsCols(inRows, inCols);
	for(unsigned int i = 0; i < mRows; ++i) {
		for(unsigned int j = 0; j < mCols; ++j) {
			(*this)(i,j) = (i < lMat.mRows && j < lMat.mCols ? lMat(i,j) : 0.);
		}
	}
}

/*!
*/
void Matrix::scaleLU(vector<double>& outScales) const
{
	outScales.resize(mCols);
	for(unsigned int i = 0; i < mRows; ++i) {
		double lMax = 0.;
		for(unsigned int j = 0; j < mCols; ++j) {
			const double lTmp = fabs((*this)(i,j));     
			if(lTmp > lMax) lMax=lTmp;
		}
		if(lMax == 0.) throw runtime_error("<Matrix::scaleLU> matrix is singular!");
		outScales[i] = 1./lMax;
	}
}

/*!
*/
void Matrix::setIdentity(unsigned int inSize)
{
	setRowsCols(inSize, inSize);
	for(unsigned int j = 0; j < mCols; ++j) {
		for(unsigned int i = 0; i < mRows; ++i) {
			(*this)(i,j) = (i == j ? 1 : 0);
		}
	}
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::subtract(Matrix& outMatrix, double inScalar) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "subtract() invalid matrix!");
	outMatrix.setRowsCols(mRows, mCols);
	for(unsigned int i = 0; i < size(); ++i) outMatrix[i] = (*this)[i] - inScalar;
	return outMatrix;
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::subtract(Matrix& outMatrix, const Matrix& inMatrix) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "subtract() invalid matrix!");
	PACC_AssertM(mRows == inMatrix.mRows && mCols == inMatrix.mCols, "subtract() matrix mismatch!");
	outMatrix.setRowsCols(mRows, mCols);
	for(unsigned int i = 0; i < size(); ++i) outMatrix[i] = (*this)[i] - inMatrix[i];
	return outMatrix;
}

/*!
 *  \param d Real part of eigenvalues computed from the matrix.
 *  \param e Imaginary part of eigenvalues computed from the matrix.
 *  \param V Eigenvectors computed from the matrix.
 *
 *  This method is derived from procedure tql2 of the Java package JAMA,
 *  which is itself derived from the Algol procedures tql2, by
 *  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
 *  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
 *  Fortran subroutine in EISPACK.
 */
void Matrix::tql2(Vector& d, Vector& e, Matrix& V) const
{
	const unsigned int n=mRows;
	for(unsigned int i = 1; i < n; i++) e[i-1] = e[i];
	e[n-1] = 0.0;
	
	double f = 0.0;
	double tst1 = 0.0;
	double eps = std::pow(2.0,-52.0);
	for(unsigned int l = 0; l < n; l++) {
		// Find small subdiagonal element
		tst1 = max(tst1, abs(d[l]) + abs(e[l]));
		unsigned int m=l;
		while((m+1) < n) {
			if(std::abs(e[m]) <= eps*tst1) break;
			m++;
		}
		
		// If m == l, d[l] is an eigenvalue,
		// otherwise, iterate.
		if(m > l) {
			unsigned int iter = 0;
			do {
				iter = iter + 1;  // (Could check iteration count here.)
													// Compute implicit shift
				double g = d[l];
				double p = (d[l+1] - g) / (2.0 * e[l]);
				double r = hypot(p,1.0);
				if(p < 0) r = -r;
				d[l] = e[l] / (p + r);
				d[l+1] = e[l] * (p + r);
				double dl1 = d[l+1];
				double h = g - d[l];
				for(unsigned int i = l+2; i < n; i++) d[i] -= h;
				f = f + h;
				
				// Implicit QL transformation.
				p = d[m];
				double c = 1.0;
				double c2 = c;
				double c3 = c;
				double el1 = e[l+1];
				double s = 0.0;
				double s2 = 0.0;
				for(unsigned int i = m-1; i >= l; i--) {
					c3 = c2;
					c2 = c;
					s2 = s;
					g = c * e[i];
					h = c * p;
					r = hypot(p,e[i]);
					e[i+1] = s * r;
					s = e[i] / r;
					c = p / r;
					p = c * d[i] - s * g;
					d[i+1] = h + s * (c * g + s * d[i]);
					
					// Accumulate transformation.
					for(unsigned int k = 0; k < n; k++) {
						h = V(k,i+1);
						V(k,i+1) = s * V(k,i) + c * h;
						V(k,i) = c * V(k,i) - s * h;
					}
					if(i == 0) break;
				}
				p = -s * s2 * c3 * el1 * e[l] / dl1;
				e[l] = s * p;
				d[l] = c * p;
				
				// Check for convergence.
			} while (std::abs(e[l]) > eps*tst1);
		}
		d[l] = d[l] + f;
		e[l] = 0.0;
	}
}

/*!
 *  \param d Real part of eigenvalues computed from the matrix.
 *  \param e Imaginary part of eigenvalues computed from the matrix.
 *  \param V Eigenvectors computed from the matrix.
 *
 *  This method is derived from procedure tred2 of the Java package JAMA,
 *  which is itself derived from the Algol procedures tred2, by
 *  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
 *  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
 *  Fortran subroutine in EISPACK.
 */
void Matrix::tred2(Vector& d, Vector& e, Matrix& V) const
{
	const unsigned int n=mRows;
	V = *this;
	for(unsigned int j = 0; j < n; ++j) d[j] = V(n-1,j);
	
	// Householder reduction to tridiagonal form.
	for(unsigned int i = n-1; i > 0; --i) {
		
		// Scale to avoid under/overflow.
		double scale = 0.0;
		double h = 0.0;
		for(unsigned int k = 0; k < i; ++k) scale += abs(d[k]);
		if(scale == 0.0) {
			e[i] = d[i-1];
			for(unsigned int j = 0; j < i; ++j) {
				d[j] = V(i-1,j);
				V(i,j) = 0.0;
				V(j,i) = 0.0;
			}
		} else {
			// Generate Householder vector.
			for(unsigned int k=0; k<i; ++k) {
				d[k] /= scale;
				h += d[k] * d[k];
			}
			double f = d[i-1];
			double g = sqrt(h);
			if(f > 0.0) g = -g;
			e[i] = scale * g;
			h = h - f * g;
			d[i-1] = f - g;
			for(unsigned int j = 0; j < i; j++) e[j] = 0.0;
			
			// Apply similarity transformation to remaining columns.
			for(unsigned int j = 0; j < i; j++) {
				f = d[j];
				V(j,i) = f;
				g = e[j] + V(j,j) * f;
				for(unsigned int k = j+1; k <= i-1; k++) {
					g += V(k,j) * d[k];
					e[k] += V(k,j) * f;
				}
				e[j] = g;
			}
			f = 0.0;
			for(unsigned int j = 0; j < i; j++) {
				e[j] /= h;
				f += e[j] * d[j];
			}
			double hh = f / (h + h);
			for(unsigned int j=0; j<i; j++) e[j] -= hh * d[j];
			for(unsigned int j=0; j<i; j++) {
				f = d[j];
				g = e[j];
				for(unsigned int k = j; k <= i-1; k++) V(k,j) -= (f * e[k] + g * d[k]);
				d[j] = V(i-1,j);
				V(i,j) = 0.0;
			}
		}
		d[i] = h;
	}
	
	// Accumulate transformations.
	for(unsigned int i = 0; i < n-1; i++) {
		V(n-1,i) = V(i,i);
		V(i,i) = 1.0;
		double h = d[i+1];
		if(h!=0.0) {
			for(unsigned int k=0; k<=i; k++) d[k] = V(k,i+1) / h;
			for(unsigned int j=0; j<=i; j++) {
				double g = 0.0;
				for(unsigned int k=0; k<=i; k++) g += V(k,i+1) * V(k,j);
				for(unsigned int k=0; k<=i; k++) V(k,j) -= g * d[k];
			}
		}
		for(unsigned int k=0; k<=i; k++) V(k,i+1) = 0.0;
	}
	for(unsigned int j=0; j<n; j++) {
		d[j] = V(n-1,j);
		V(n-1,j) = 0.0;
	}
	V(n-1,n-1) = 1.0;
	e[0] = 0.0;
}

/*!
*/
Matrix Matrix::transpose(void) const 
{
	Matrix lMatrix; 
	return transpose(lMatrix);
}

/*!
This method also returns a reference to the result.
 */
Matrix& Matrix::transpose(Matrix& outMatrix) const
{
	PACC_AssertM(mRows > 0 && mCols > 0, "transpose() invalid matrix!");
	if(&outMatrix != this) {
		// output matrix is not self assigning
		outMatrix.setRowsCols(mCols,mRows);
		// transpose elements
		for(unsigned int i = 0; i< mRows; ++i) {
			for(unsigned int j = 0; j < mCols; ++j) {
				outMatrix(j,i) = (*this)(i,j);
			}
		}
	} else {
		// use temporary matrix to self assign
		Matrix lMatrix(*this);
		outMatrix.setRowsCols(mCols,mRows);
		// transpose elements
		for(unsigned int i = 0; i< mRows; ++i) {
			for(unsigned int j = 0; j < mCols; ++j) {
				outMatrix(j,i) = lMatrix(i,j);
			}
		}
	}		
	return outMatrix;
}

/*!
*/
void Matrix::throwError(const string& inMessage, const XML::Iterator& inNode) const
{
	ostringstream lStream;
	lStream << inMessage << " for markup:\n";
	XML::Streamer lStreamer(lStream);
	inNode->serialize(lStreamer);
	throw runtime_error(lStream.str());
}

/*!
 See Matrix::read for a description of the write format. By default, the precision 
 of the output is set to 15 digits. This value can be changed using method
 Matrix::setPrecision.
*/
void Matrix::write(XML::Streamer& outStream, const string& inTag) const
{
	outStream.openTag(inTag, false);
	if(mName != "") outStream.insertAttribute("name", mName);
	outStream.insertAttribute("rows", mRows);
	outStream.insertAttribute("cols", mCols);
	ostringstream lContent;
	lContent.precision(mPrec);
	for(unsigned int i = 0; i < size(); ++i) {
		if(i != 0 && i % mCols == 0) lContent << ";";
		else if(i != 0) lContent << ",";
		lContent << (*this)[i];
	}
	outStream.insertStringContent(lContent.str());
	outStream.closeTag();
}

/*!
*/
ostream& PACC::operator<<(ostream &outStream, const Matrix& inMatrix)
{
	XML::Streamer lStream(outStream);
	inMatrix.write(lStream);
	return outStream;
}

/*!
This method uses the first data tag of the parse tree to read the matrix. The 
corresponding tree root is then erased. Any read error throws a std::runtime_error.
*/
XML::Document& PACC::operator>>(XML::Document& inDocument, Matrix& outMatrix)
{
	XML::Iterator lNode = inDocument.getFirstDataTag();
	outMatrix.read(lNode);
	inDocument.erase(lNode);
	return inDocument;
}