MTL 矩阵逆阵 解线性方程

/* thanks to Valient Gough for this example program! */

//整理 by RobinKin

#include <mtl/matrix.h>
#include <mtl/mtl.h>
#include <mtl/utils.h>
#include <mtl/lu.h>

using namespace mtl;

// don't print out the matrices once they get to this size...
#define MAX_PRINT_SIZE 5

typedef matrix<double, rectangle<>, dense<>, row_major>::type Matrix;
typedef dense1D<double> Vector;

double testMatrixError(const Matrix &A, const Matrix &AInv)
{
  int size = A.nrows();

  // test it
  Matrix AInvA(size,size);

  // AInvA = AInv * A
  mult(AInv, A, AInvA);

  // I = identity
  typedef matrix<double, diagonal<>, packed<>, row_major>::type IdentMat;
  IdentMat I(size, size, 0, 0);
  mtl::set_value(I, 1.0);

  // AInvA += -I
  add(scaled(I, -1.0), AInvA);