MTL 矩阵逆阵 解线性方程

原创 2005年03月04日 13:42:00

/* 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);

  if (size < MAX_PRINT_SIZE) {
    std::cout << "Ainv * A - I = " << std::endl;
    print_all_matrix(AInvA);
  }

  // find max error
  double max_error = 0.0;
  for(Matrix::iterator i = AInvA.begin(); i != AInvA.end(); ++i)
    for(Matrix::Row::iterator j = (*i).begin(); j != (*i).end(); ++j)
      if(fabs(*j) > fabs(max_error))
 max_error = *j;
       
  std::cout << "max error = " << max_error << std::endl;

  return max_error;
}


void testLUSoln(const Matrix &A, const Vector &b, Vector &x)
{
  // create LU decomposition
  Matrix LU(A.nrows(), A.ncols());
  dense1D<int> pvector(A.nrows());

  copy(A, LU);
  lu_factorize(LU, pvector);
       
  // solve
  //解线形方程
  lu_solve(LU, pvector, b, x);
}

void testLUInv(const Matrix &A, int size)
{
  // invert it
  Matrix AInv(size,size);
       
  // create LU decomposition
  Matrix LU(A.nrows(), A.ncols());
  dense1D<int> pvector(A.nrows());

  copy(A, LU);
  lu_factor(LU, pvector);
       
        //求逆阵
  // solve
  lu_inverse(LU, pvector, AInv);
       

  if(size < MAX_PRINT_SIZE) {
    std::cout << "Ainv = " << std::endl;
    print_all_matrix(AInv);
  }

  // test it
  testMatrixError(A, AInv);

}

int main(int argc, char **argv)
{
  typedef Matrix::size_type sizeT;

  sizeT size = 3;

  if(argc > 1)
    size = atoi(argv[1]);

  std::cout << "inverting matrix of size " << size << std::endl;

  // create a random matrix and invert it.  Then see how close it comes to
  // identity.

  Matrix A(size,size);
  Vector b(size);
  Vector x(size);

  // initialize
  for (sizeT i=0; i<A.nrows(); i++) {
    for (sizeT j=0; j<A.nrows(); j++)
      A(i,j) = (double)(rand() % 200 - 100) / 50.0;
    b[i] = (double)(rand() % 200 - 100) / 50.0;
  }

  if (size < MAX_PRINT_SIZE) {
    std::cout << "A = " << std::endl;
    print_all_matrix(A);
  }

       
  // time LU inv
  std::cout << std::endl
       << " ----------- testing inversion using LU decomposition"
       << std::endl;
  testLUInv(A, size);

  if (size < MAX_PRINT_SIZE) {
    std::cout << "solution = ";
    print_vector(x);
  }
       
  // test LU solution
  mtl::set_value(x, 0.0);
  testLUSoln(A, b, x);

  if(size < MAX_PRINT_SIZE) {
    std::cout << "solution = ";
    print_vector(x);
  }

  if(argc == 1)
    std::cout << std::endl
  << "pass size argument to program to time larger matrices."
  << std::endl;


  return 0;
}

输出:
inverting matrix of size 3
A =
3x3
[
[1.66,-0.28,1.54],
[1.86,0.7,1.72],
[-1.02,-1.58,1.24]
]

//逆阵

 ----------- testing inversion using LU decomposition
Ainv =
3x3
[
[0.978889,-0.56949,-0.42578],
[-1.10862,0.990792,0.00251165],
[-0.607383,0.79401,0.459414]
]
Ainv * A - I =
3x3
[
[4.44089e-16,2.56739e-16,1.01481e-16],
[-5.05491e-16,-2.22045e-16,-2.75514e-16],
[-1.91335e-16,-1.27014e-16,-1.11022e-16]
]
max error = -5.05491e-16
solution = [0,0,0,]
//方程的解向量
solution = [1.00642,-0.49478,-0.980001,]

Delphi7高级应用开发随书源码

  • 2003年04月30日 00:00
  • 676KB
  • 下载

【线性代数公开课MIT Linear Algebra】 第十六课 Ax=b的解、最小二乘法与矩阵

本系列笔记为方便日后自己查阅而写,更多的是个人见解,也算一种学习的复习与总结,望善始善终吧~ 基本子空间与投影矩阵 上一节课我们已经了解了投影矩阵 projection matrix, P=...
  • a352611
  • a352611
  • 2015年11月04日 22:14
  • 644

异想家纯C语言矩阵运算库

Sandeepin最近做的项目中需要在嵌入式芯片里跑一些算法,而这些单片机性能不上不下,它能跑些简单的程序,但又还没到上Linux系统的地步。所以只好用C语言写一些在高级语言里一个函数就解决的算法了,...

选主元的高斯-约当(Gauss-Jordan)消元法解线性方程组和求逆矩阵

选主元的高斯-约当(Gauss-Jordan)消元法在很多地方都会用到,例如求一个矩阵的逆矩阵、解线性方程组(插一句:LM算法求解的一个步骤),等等。它的速度不是最快的,但是它非常稳定(来自网上的定义...

矩阵求逆和解线性方程c源代码

  • 2016年03月26日 13:10
  • 3KB
  • 下载

数值计算-线性方程组求解(2)-追赶法解三对角矩阵-MATLAB实现

问题描述  在用差分法求解二阶常微分方程的边值问题、热传导问题及三次样条插值函数的求解问题中,都会遇到下面形式的阶数较高的三对角方程组AX=fAX=f,即 ⎡⎣⎢⎢⎢⎢⎢⎢⎢b1a2⋱c1b2an−...

MTL矩阵模板库

  • 2014年07月28日 17:27
  • 1.36MB
  • 下载

求解线性方程组之Gauss-Jordan消去法求矩阵的逆

//用具有最大主元素的Gauss-Jordan消去法求矩阵的逆 #include #include #include #include using namespace std; c...

LUP分解求解线性方程组及求逆矩阵 java

具体的分析参考算法导论 public class LU { /**矩阵相乘*/ public static float[][] multiply(float[][]A,float[][]B){...
  • cumtwyc
  • cumtwyc
  • 2015年11月22日 19:07
  • 954
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:MTL 矩阵逆阵 解线性方程
举报原因:
原因补充:

(最多只允许输入30个字)