Stupid Matrix class test

//          Copyright Nicolas Janin 2004 - 20013.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)


// Test case for http://pastebin.com/UE80Dpm1

#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
#include "SMatrix.hpp"

using namespace std;


typedef SMatrix<double> dmat;

// matrix display
template <typename T>
ostream& operator<<(ostream &os, SMatrix<T> &mat) {
	os << "(" << mat.rows() << "," << mat.cols() << ") =\n";
	os.setf(ios::fixed, ios::floatfield);
	os << setprecision(4);
	size_t L = mat.rows(), C = mat.cols();
	for(size_t l = 0; l < L; ++l) {
		os << "[ ";
		size_t c = 0;
		for(; c < C-1; ++c) os << setw(5) << mat(l, c) << '\t';
		os << setw(5) << mat(l, c) << " ]" << endl;
	}
	return os;
}

double f1(double x){ return sqrt(fabs(x)); }
double f2(double x, int i){ return (x * i); }


int main(int argc, char* argv[])
{

	{
		dmat M(5, 5);
		cout << "M" << M << endl;
		M.fill(-1);
		cout << "M" << M << endl;

		cout << "Divers acces en ecriture" << endl;
		for(size_t i = 1; i < M.rows(); ++i)
			for(size_t j = 0; j < M.cols(); ++j) M(i-1,j) = i / (j+1.);
		cout << "M" << M << endl;

		cout << "M.tr()" << M.tr() << endl;

		M.apply(&f1);
		cout << "M.apply(sqrt(fabs(elem))) " << M << endl;
		M.apply(&f2, 2);
		cout << "M x 2 " << M << endl;
	}

	// Test dot product
	{
		double a[] = {1,2,0, 4,3,-1};
		double b[] = {5,1, 2,3, 3,4};
		dmat A(2,3,a);
		dmat AA(A);
		dmat B(3,2,b);
		cout << "A" << A << endl;
		cout << "B" << B << endl;
		A.dot(B);
		cout << "A.dot(B) = \n" << A << endl;
		B.dot(AA);
		cout << "B.dot(A) = \n" << B << endl;
	}

	{
		dmat A(4,4); dmat B(4,4);
		for (int i = 0; i < A.rows(); i++)
			for (int j = 0; j < A.cols(); j++) {
				A(i,j) = i + j;
				B(i,j) = i - j;
			}
		cout << "Check the formula tr(AB) = tr(B).tr(A), with (" << A.rows() << "," << A.cols() << ") matrices" << endl;
		cout << "A" << A;
		cout << "B" << B;
		int t = clock();
		dmat AA(A);
		(A.dot(B)).tr();
		(B.tr()).dot(AA.tr());
		cout << "A" << A;
		cout << "B" << B;
		// check equality
		for (size_t i = 0; i < A.rows(); i++)
			for (size_t j = 0; j < A.cols(); j++)
				assert(A(i,j) == B(i,j));
		cout << "Done in " << (clock() - t)/1000 << " ms" <<endl;
	}

	{
		cout << "Test arithmetic operations with scalars" << endl;
		double a[] = {1,2,3, 4,5,6};
		double b[] = {1,2, 3,4, 5,6};
		dmat A(2,3,a);
		dmat B(3,2,b);

		cout << "A + 3" << A+3 << endl;
		cout << "A - 3" << A-3 << endl;
		cout << "B x 2" << B*2 << endl;
		cout << "B / 2" << B/2 << endl;
	}

	{
		cout << "Test submatrix" << endl;
		double a[] = {1,2,3,4, 1,2,3,4, 1,2,3,4, 1,2,3,4};
		dmat A(4,4,a);
		dmat B(3,3);

		A.submatrix(B,1,1,3,3);
		cout << "B" << B << endl;

		cout << "Extract row and column" << endl;
		dmat row(1,4);
		dmat col(4,1);
		A.row(row, 2);
		A.column(col, 3);
		cout << "A.row(2) " << row << endl;
		cout << "A.col(3) " << col << endl;
	}

	return 0;
}