MTL 矩阵逆阵 解线性方程
2005-03-04 13:42
330 查看
/* 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,]
//整理 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,]
相关文章推荐
- 逆矩阵及其矩阵线性一元方程的简单算法
- 杭电_hdu1757_矩阵解线性方程_快速幂乘
- MTL 解下三角带状矩阵线形方程
- 线性代数复习之003矩阵初等变换与方程
- n阶常系数线性递推方程的矩阵乘方解法
- 矩阵表示的多项式和线性空间的关系
- 梯度、Hessian矩阵、平面方程的法线以及函数导数的含义
- 线性代数笔记(特征问题与矩阵相似)
- 线性代数基础知识(三)—— 矩阵乘法
- 稀疏矩阵线性表示
- 线性方程求法分类
- codeforces 801 E. Vulnerable Kerbals(DAG最长路+解线性同于方程)
- hdu 5015 233 Matrix 线性序列构造矩阵快速幂
- 线性代数复习-矩阵及其运算
- 线性代数标准型矩阵化简技巧
- 数据结构(二):线性表的使用原则以及链表的应用-稀疏矩阵的三元组表示
- POJ2115——C Looooops(扩展欧几里德+求解模线性方程)
- poj 2115 C Looooops 拓展欧几里德 解模线性方程模板
- 线性代数-矩阵及其运算(总结)
- [转]用矩阵相乘解线性递推方程