main.cpp

强大的矩阵模版类

//#include "mtl/mtl.h"                                                                                        
//
//using namespace std;
//using namespace mtl;
//
//template <class Matrix>
//void print_matrix(Matrix& mat, const string& description)
//{
//	std::cout << description;
//	std::cout << '[';
//	for (Matrix::iterator i = mat.begin(); i!=mat.end(); ++i)
//	{
//		for (Matrix::OneD::iterator j = (*i).begin(); j!=(*i).end(); ++j)
//		{
//			std::cout << '\t' << *j;
//		}
//		std::cout << ((i+1 == mat.end()) ? "\t]\n" :  "\n");
//	}
//}
//
//int main(int argc, char* argv[])
//{
//
//	typedef matrix<float, rectangle<>, dense<>, row_major>::type Matrix;
//
//	const Matrix::size_type MAX_ROW = 3, MAX_COL = 3;
//	Matrix  A(MAX_ROW,MAX_COL),B(MAX_ROW,MAX_COL),C(MAX_ROW,MAX_COL);
//
//	// fill Matrix A with the index syntax
//	for (Matrix::size_type i=0; i<MAX_ROW; ++i)
//	{
//		for (Matrix::size_type j=0; j<MAX_COL; ++j)
//		{
//			A(i, j) = Matrix::value_type(rand() % 50);
//		}
//	}
//
//	// fill Matrix B with the iterator syntax
//	for (Matrix::iterator i=B.begin(); i!=B.end(); ++i)
//	{
//		for (Matrix::OneD::iterator j=(*i).begin(); j!=(*i).end(); ++j)
//		{
//			*j = Matrix::value_type(rand() % 50);
//		}
//	}
//
//	print_matrix(A, "A=\n");
//	print_matrix(B, "B=\n");
//
//	// Matrix C = A + B
//	add(A, C);
//	add(B,C);
//	print_matrix(C, "C = A + B  \n");
//
//	// Matrix C = A * B^T,  B^T: transpose of B
//	transpose(B);
//	print_matrix(B, "B^T=\n");
//	zero_matrix(C);        // because mult(A, B, C): C += A*B
//	mult(A,B,C);
//
//	print_matrix(C, "C = A * B^T\n");
//
//	getchar();
//	return 0;
//}


//=======================================================================

#include "mtl/mtl.h"
#include "mtl/lu.h"

using namespace std;
using namespace mtl;

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

{

	typedef matrix<float, rectangle<>, dense<external>, row_major>::type Matrix;

	// dense : data copy from a float array,not generate them with yourself
	const Matrix::size_type MAX_ROW = 3, MAX_COL = 3;

	// solve the equation Ax=b
	// { 4x - y + z = 7
	//    4x - 8y + z= -21
	//   -2x + y + 5z = 15 }
	// A = [ 4 -1 1
	//          4 -8 1
	//          -2 1 5 ]
	// b = [7 - 21 15]^T

	float a[] = {4.0f, -1.0f, 1.0f, 4.0f, -8.0f, 1.0f, -2.0f, 1.0f, 5.0f};
	Matrix A(a, MAX_ROW, MAX_COL);

	typedef matrix<float, rectangle<>, dense<>, row_major>::type LUMatrix;
	LUMatrix LU(A.nrows(), A.ncols());
	mtl::copy(A, LU);

	typedef dense1D<float> Vector;
	Vector pvector(A.nrows());
	lu_factor(LU, pvector);

	Vector b(A.nrows()), x(A.nrows());
	b[0] = 7.0f, b[1] = -21.0f, b[2] = 15.0f;
	lu_solve(LU, pvector, b, x);

	for (Vector::iterator i=x.begin(); i!=x.end(); ++i)
		cout << *i << '\t';

	system("pause");
	return 0;
}

//=======================================================================


//#include "mtl/mtl.h"
//#include "mtl/lu.h"
//
//using namespace std;
//using namespace mtl;
//
//template <class Matrix>
//void print_matrix(Matrix& mat, const string& description)
//{                                                                                             
//	std::cout << description;
//	std::cout << '[';
//
//	for (Matrix::iterator i = mat.begin(); i!=mat.end(); ++i)
//	{
//		for (Matrix::OneD::iterator j = (*i).begin(); j!=(*i).end(); ++j)
//		{
//			std::cout << '\t' << *j;
//		}
//		std::cout << ((i+1 == mat.end()) ? "\t]\n" :  "\n");
//	}
//}
//
//int main(int argc, char* argv[])
//{
//	typedef matrix<double, rectangle<>, dense<external>, row_major>::type Matrix;
//
//	// dense : data copy from a float array,not generate them with yourself
//	const Matrix::size_type MAX_ROW = 3, MAX_COL = 3;
//
//	// inverse matrix A
//	// A = [ 4 -1 1
//	//          4 -8 1
//	//          -2 1 5 ]
//	double a[] = {4.0, -1.0, 1.0, 4.0, -8.0, 1.0, -2.0, 1.0, 5.0};
//	Matrix A(a, MAX_ROW, MAX_COL);
//
//	typedef matrix<double, rectangle<>, dense<>, row_major>::type CMatrix;
//	CMatrix LU(A.nrows(), A.ncols());
//
//	mtl::copy(A, LU);
//
//	typedef dense1D<double> Vector;
//	Vector pvector(A.nrows());
//	lu_factor(LU, pvector);
//
//	CMatrix InvA(A.nrows(), A.ncols()); 
//	lu_inverse(LU, pvector, InvA);
//
//	print_matrix(A, "A = \n");
//	print_matrix(InvA, "A^(-1) = \n");
//
//	system("pause");
//	return 0;
//
//}