头文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* linequs1.h
*
* Function template for solving deterministic linear equations.
*
* For a n-by-n coefficient matrix A and n-by-1 constant vector b, Gauss
* elimination method can solve the equations Ax=b. But the decomposition
* methods are more commonly used. if A is not SPD, the LU decomposition
* can be use to solve the equations; if A is SPD, the Cholesky decomposition
* is the better choice; if A is a tridiagonal matrix, then the Forward
* Elimination and Backward Substitution maybe the best choice.
*
* These functions can be used by both REAL or COMPLEX linear equations.
*
* Zhang Ming, 2010-07 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
#ifndef DETLINEQUS_H
#define DETLINEQUS_H
#include <vector.h>
#include <matrix.h>
#include <cholesky.h>
#include <lud.h>
namespace splab
{
template<typename Type>
Vector<Type> gaussSolver( const Matrix<Type>&, const Vector<Type>& );
template<typename Type>
void gaussSolver( Matrix<Type>&, Matrix<Type>& );
template<typename Type>
Vector<Type> luSolver( const Matrix<Type>&, const Vector<Type>& );
template<typename Type>
Matrix<Type> luSolver( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Vector<Type> choleskySolver( const Matrix<Type>&, const Vector<Type>& );
template<typename Type>
Matrix<Type> choleskySolver( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Vector<Type> utSolver( const Matrix<Type>&, const Vector<Type>& );
template<typename Type>
Vector<Type> ltSolver( const Matrix<Type>&, const Vector<Type>& );
template<typename Type>
Vector<Type> febsSolver( const Vector<Type>&, const Vector<Type>&,
const Vector<Type>&, const Vector<Type>& );
#include <linequs1-impl.h>
}
// namespace splab
#endif
// DETLINEQUS_H
实现文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* linequs1-impl.h
*
* Implementation for deterministic linear equations.
*
* Zhang Ming, 2010-07 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
/**
* Linear equations solution by Gauss-Jordan elimination, Ax=B.
* A ---> The n-by-n coefficient matrix(Full Rank);
* b ---> The n-by-1 right-hand side vector;
* x ---> The n-by-1 solution vector.
*/
template <typename Type>
Vector<Type> gaussSolver( const Matrix<Type> &A, const Vector<Type> &b )
{
assert( A.rows() == b.size() );
int rows = b.size();
Vector<Type> x( rows );
Matrix<Type> C(A);
Matrix<Type> B( rows, 1 );
for( int i=0; i<rows; ++i )
B[i][0] = b[i];
gaussSolver( C, B );
for( int i=0; i<rows; ++i )
x[i] = B[i][0];
return x;
}
/**
* Linear equations solution by Gauss-Jordan elimination, AX=B.
* A ---> The n-by-n coefficient matrix(Full Rank);
* B ---> The n-by-m right-hand side vectors;
* When solving is completion, A is replaced by its inverse and
* B is replaced by the corresponding set of solution vectors.
* Adapted from Numerical Recipes.
*/
template <typename Type>
void gaussSolver( Matrix<Type> &A, Matrix<Type> &B )
{
assert( A.rows() == B.rows() );
int i, j, k, l, ll, icol, irow,
n=A.rows(),
m=B.cols();
Type big, dum, pivinv;
Vector<int> indxc(n), indxr(n), ipiv(n);
for( j=0; j<n; j++ )
ipiv[j]=0;
for( i=0; i<n; ++i )
{
big = 0.0;
for( j=0; j<n; ++j )
if( ipiv[j] != 1 )
for( k=0; k<n; ++k )
if( ipiv[k] == 0 )
if( abs(A[j][k]) >= abs(big) )
{
big = abs(A[j][k]);
irow = j;
icol = k;
}
++(ipiv[icol]);
if( irow != icol )
{
for( l=0; l<n; ++l )
swap( A[irow][l], A[icol][l] );
for( l=0; l<m; ++l )
swap( B[irow][l], B[icol][l] );
}
indxr[i] = irow;
indxc[i] = icol;
if( abs(A[icol][icol]) == 0.0 )
{
cerr << "Singular Matrix!" << endl;
return;
}
pivinv = Type(1) / A[icol][icol];
A[icol][icol] = Type(1);
for( l=0; l<n; ++l )
A[icol][l] *= pivinv;
for( l=0; l<m; ++l )
B[icol][l] *= pivinv;
for( ll=0; ll<n; ++ll )
if( ll != icol )
{
dum = A[ll][icol];
A[ll][icol] = 0;
for( l=0; l<n; ++l )
A[ll][l] -= A[icol][l]*dum;
for( l=0; l<m; ++l )
B[ll][l] -= B[icol][l]*dum;
}
}
for( l=n-1; l>=0; --l )
if( indxr[l] != indxc[l] )
for( k=0; k<n; k++ )
swap( A[k][indxr[l]], A[k][indxc[l]] );
}
/**
* Linear equations solution by LU decomposition, Ax=b.
* A ---> The n-by-n coefficient matrix(Full Rank);
* b ---> The n-by-1 right-hand side vector;
* x ---> The n-by-1 solution vector.
*/
template <typename Type>
Vector<Type> luSolver( const Matrix<Type> &A, const Vector<Type> &b )
{
assert( A.rows() == b.size() );
LUD<Type> lu;
lu.dec( A );
return lu.solve( b );
}
/**
* Linear equations solution by LU decomposition, AX=B.
* A ---> The n-by-n coefficient matrix(Full Rank);
* B ---> The n-by-m right-hand side vector;
* X ---> The n-by-m solution vectors.
*/
template <typename Type>
Matrix<Type> luSolver( const Matrix<Type> &A, const Matrix<Type> &B )
{
assert( A.rows() == B.rows() );
LUD<Type> lu;
lu.dec( A );
return lu.solve( B );
}
/**
* Linear equations solution by Cholesky decomposition, Ax=b.
* A ---> The n-by-n coefficient matrix(Full Rank);
* b ---> The n-by-1 right-hand side vector;
* x ---> The n-by-1 solution vector.
*/
template <typename Type>
Vector<Type> choleskySolver( const Matrix<Type> &A, const Vector<Type> &b )
{
assert( A.rows() == b.size() );
Cholesky<Type> cho;
cho.dec(A);
if( cho.isSpd() )
return cho.solve( b );
else
{
cerr << "Factorization was not complete!" << endl;
return Vector<Type>(0);
}
}
/**
* Linear equations solution by Cholesky decomposition, AX=B.
* A ---> The n-by-n coefficient matrix(Full Rank);
* B ---> The n-by-m right-hand side vector;
* X ---> The n-by-m solution vectors.
*/
template <typename Type>
Matrix<Type> choleskySolver( const Matrix<Type> &A, const Matrix<Type> &B )
{
assert( A.rows() == B.rows() );
Cholesky<Type> cho;
cho.dec(A);
if( cho.isSpd() )
return cho.solve( B );
else
{
cerr << "Factorization was not complete!" << endl;
return Matrix<Type>(0,0);
}
}
/**
* Solve the upper triangular system U*x = b.
* U ---> The n-by-n upper triangular matrix(Full Rank);
* b ---> The n-by-1 right-hand side vector;
* x ---> The n-by-1 solution vector.
*/
template <typename Type>
Vector<Type> utSolver( const Matrix<Type> &U, const Vector<Type> &b )
{
int n = b.dim();
assert( U.rows() == n );
assert( U.rows() == U.cols() );
Vector<Type> x(b);
for( int k=n; k >= 1; --k )
{
x(k) /= U(k,k);
for( int i=1; i<k; ++i )
x(i) -= x(k) * U(i,k);
}
return x;
}
/**
* Solve the lower triangular system L*x = b.
* L ---> The n-by-n lower triangular matrix(Full Rank);
* b ---> The n-by-1 right-hand side vector;
* x ---> The n-by-1 solution vector.
*/
template <typename Type>
Vector<Type> ltSolver( const Matrix<Type> &L, const Vector<Type> &b )
{
int n = b.dim();
assert( L.rows() == n );
assert( L.rows() == L.cols() );
Vector<Type> x(b);
for( int k=1; k <= n; ++k )
{
x(k) /= L(k,k);
for( int i=k+1; i<= n; ++i )
x(i) -= x(k) * L(i,k);
}
return x;
}
/**
* Tridiagonal equations solution by Forward Elimination and Backward Substitution.
* a ---> The n-1-by-1 main(0th) diagonal vector;
* b ---> The n-by-1 -1th(above main) diagonal vector;
* c ---> The n-1-by-1 +1th(below main) diagonal vector;
* d ---> The right-hand side vector;
* x ---> The n-by-1 solution vector.
*/
template <typename Type>
Vector<Type> febsSolver( const Vector<Type> &a, const Vector<Type> &b,
const Vector<Type> &c, const Vector<Type> &d )
{
int n = b.size();
assert( a.size() == n-1 );
assert( c.size() == n-1 );
assert( d.size() == n );
Type mu = 0;
Vector<Type> x(d),
bb(b);
for( int i=0; i<n-1; ++i )
{
mu = a[i] / bb[i];
bb[i+1] = bb[i+1] - mu*c[i];
x[i+1] = x[i+1] - mu*x[i];
}
x[n-1] = x[n-1] / bb[n-1];
for( int i=n-2; i>=0; --i )
x[i] = ( x[i]-c[i]*x[i+1] ) / bb[i];
// for( int j=1; j<n; ++j )
// {
// mu = a(j) / bb(j);
// bb(j+1) = bb(j+1) - mu*c(j);
// x(j+1) = x(j+1) - mu*x(j);
// }
//
// x(n) = x(n) / bb(n);
// for( int j=n-1; j>0; --j )
// x(j) = ( x(j)-c(j)*x(j+1) ) / bb(j);
return x;
}
测试代码：
/*****************************************************************************
* linequs1_test.cpp
*
* Deterministic Linear Equations testing.
*
* Zhang Ming, 2010-07 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <iomanip>
#include <linequs1.h>
using namespace std;
using namespace splab;
typedef float Type;
const int M = 3;
const int N = 3;
int main()
{
Matrix<Type> A(M,N), B;
Vector<Type> b(N);
// ordinary linear equations
A[0][0] = 1; A[0][1] = 2; A[0][2] = 1;
A[1][0] = 2; A[1][1] = 5; A[1][2] = 4;
A[2][0] = 1; A[2][1] = 1; A[2][2] = 0;
B = eye( N, Type(1.0) );
b[0] = 1; b[1] = 0; b[2] = 1;
Matrix<Type> invA( A );
Matrix<Type> X( B );
cout << setiosflags(ios::fixed) << setprecision(4) << endl;
cout << "The original matrix A is : " << A << endl;
gaussSolver( invA, X );
cout << "The inverse of A is (Gauss Solver) : "
<< invA << endl;
cout << "The inverse matrix of A is (LUD Solver) : "
<< luSolver( A, B ) << endl;
cout << "The constant vector b : " << b << endl;
cout << "The solution of A * x = b is (Gauss Solver) : "
<< gaussSolver( A, b ) << endl;
cout << "The solution of A * x = b is (LUD Solver) : "
<< luSolver( A, b ) << endl << endl;
Matrix<complex<Type> > cA = complexMatrix( A, B );
Vector<complex<Type> > cb = complexVector( b, b );
cout << "The original complex matrix cA is : " << cA << endl;
cout << "The constant complex vector cb : " << cb << endl;
cout << "The solution of cA * cx = cb is (Gauss Solver) : "
<< gaussSolver( cA, cb ) << endl;
cout << "The solution of cA * cx = cb is (LUD Solver) : "
<< luSolver( cA, cb ) << endl;
cout << "The cA*cx - cb is : "<< cA*luSolver(cA,cb) - cb << endl << endl;
// linear equations with symmetric coefficient matrix
for( int i=1; i<N+1; ++i )
{
for( int j=1; j<N+1; ++j )
if( i == j )
A(i,i) = Type(i);
else
if( i < j )
A(i,j) = Type(i);
else
A(i,j) = Type(j);
b(i) = Type( i*(i+1)/2.0 + i*(N-i) );
}
cout << "The original matrix A : " << A << endl;
cout << "The inverse matrix of A is (Cholesky Solver) : "
<< choleskySolver( A, B ) << endl;
cout << "The constant vector b : " << b << endl;
cout << "The solution of Ax = b is (Cholesky Solver) : "
<< choleskySolver( A, b ) << endl << endl;
cA = complexMatrix( A );
cb = complexVector( b, b );
cout << "The original complex matrix A : " << cA << endl;
cout << "The constant complex vector b : " << cb << endl;
cout << "The solution of Ax = b is (Cholesky Solver) : "
<< choleskySolver( cA, cb ) << endl << endl;
// upper and lower triangular system
for( int i=0; i<N; ++i )
for( int j=i; j<N; ++j )
B[i][j] = A[i][j];
cout << "The original matrix B : " << B << endl;
cout << "The constant vector b : " << b << endl;
cout << "The solution of Ax = b is (Upper Triangular Solver) : "
<< utSolver( B, b ) << endl << endl;
cA = complexMatrix( trT(B), -trT(B) );
cout << "The original complex matrix A : " << cA << endl;
cout << "The constant complex vector b : " << cb << endl;
cout << "The solution of Ax = b is (Lower Triangular Solver) : "
<< ltSolver( cA, cb ) << endl << endl;
// Tridiagonal linear equations
Vector<Type> aa( 3, -1 ), bb( 4, 4 ), cc( 3, -2 ), dd( 4 );
dd(1) = 3; dd(2) = 2; dd(3) = 2; dd(4) = 3;
cout << "The elements below main diagonal is : " << aa << endl;
cout << "The elements on main diagonal is : " << bb << endl;
cout << "The elements above main diagonal is : " << cc << endl;
cout << "The elements constant vector is : " << dd << endl;
cout << "Teh solution is (Forward Elimination and Backward Substitution) : "
<< febsSolver( aa, bb, cc, dd ) << endl;
Vector<complex<Type> > caa = complexVector(aa);
Vector<complex<Type> > cbb = complexVector(bb);
Vector<complex<Type> > ccc = complexVector(cc);
Vector<complex<Type> > cdd = complexVector(Vector<Type>(dd.dim()),dd);
cout << "The elements below main diagonal is : " << caa << endl;
cout << "The elements on main diagonal is : " << cbb << endl;
cout << "The elements above main diagonal is : " << ccc << endl;
cout << "The elements constant vector is : " << cdd << endl;
cout << "Teh solution is (Forward Elimination and Backward Substitution) : "
<< febsSolver( caa, cbb, ccc, cdd ) << endl;
return 0;
}
运行结果：
The original matrix A is : size: 3 by 3
1.0000 2.0000 1.0000
2.0000 5.0000 4.0000
1.0000 1.0000 0.0000
The inverse of A is (Gauss Solver) : size: 3 by 3
-4.0000 1.0000 3.0000
4.0000 -1.0000 -2.0000
-3.0000 1.0000 1.0000
The inverse matrix of A is (LUD Solver) : size: 3 by 3
-4.0000 1.0000 3.0000
4.0000 -1.0000 -2.0000
-3.0000 1.0000 1.0000
The constant vector b : size: 3 by 1
1.0000
0.0000
1.0000
The solution of A * x = b is (Gauss Solver) : size: 3 by 1
-1.0000
2.0000
-2.0000
The solution of A * x = b is (LUD Solver) : size: 3 by 1
-1.0000
2.0000
-2.0000
The original complex matrix cA is : size: 3 by 3
(1.0000,1.0000) (2.0000,0.0000) (1.0000,0.0000)
(2.0000,0.0000) (5.0000,1.0000) (4.0000,0.0000)
(1.0000,0.0000) (1.0000,0.0000) (0.0000,1.0000)
The constant complex vector cb : size: 3 by 1
(1.0000,1.0000)
(0.0000,0.0000)
(1.0000,1.0000)
The solution of cA * cx = cb is (Gauss Solver) : size: 3 by 1
(0.4000,-0.8000)
(-0.4000,1.2000)
(0.6000,-1.0000)
The solution of cA * cx = cb is (LUD Solver) : size: 3 by 1
(0.4000,-0.8000)
(-0.4000,1.2000)
(0.6000,-1.0000)
The cA*cx - cb is : size: 3 by 1
(-0.0000,0.0000)
(0.0000,0.0000)
(-0.0000,-0.0000)
The original matrix A : size: 3 by 3
1.0000 1.0000 1.0000
1.0000 2.0000 2.0000
1.0000 2.0000 3.0000
The inverse matrix of A is (Cholesky Solver) : size: 3 by 3
2.0000 -1.0000 0.0000
-1.0000 2.0000 -1.0000
0.0000 -1.0000 1.0000
The constant vector b : size: 3 by 1
3.0000
5.0000
6.0000
The solution of Ax = b is (Cholesky Solver) : size: 3 by 1
1.0000
1.0000
1.0000
The original complex matrix A : size: 3 by 3
(1.0000,0.0000) (1.0000,0.0000) (1.0000,0.0000)
(1.0000,0.0000) (2.0000,0.0000) (2.0000,0.0000)
(1.0000,0.0000) (2.0000,0.0000) (3.0000,0.0000)
The constant complex vector b : size: 3 by 1
(3.0000,3.0000)
(5.0000,5.0000)
(6.0000,6.0000)
The solution of Ax = b is (Cholesky Solver) : size: 3 by 1
(1.0000,1.0000)
(1.0000,1.0000)
(1.0000,1.0000)
The original matrix B : size: 3 by 3
1.0000 1.0000 1.0000
0.0000 2.0000 2.0000
0.0000 0.0000 3.0000
The constant vector b : size: 3 by 1
3.0000
5.0000
6.0000
The solution of Ax = b is (Upper Triangular Solver) : size: 3 by 1
0.5000
0.5000
2.0000
The original complex matrix A : size: 3 by 3
(1.0000,-1.0000) (0.0000,-0.0000) (0.0000,-0.0000)
(1.0000,-1.0000) (2.0000,-2.0000) (0.0000,-0.0000)
(1.0000,-1.0000) (2.0000,-2.0000) (3.0000,-3.0000)
The constant complex vector b : size: 3 by 1
(3.0000,3.0000)
(5.0000,5.0000)
(6.0000,6.0000)
The solution of Ax = b is (Lower Triangular Solver) : size: 3 by 1
(0.0000,3.0000)
(0.0000,1.0000)
(0.0000,0.3333)
The elements below main diagonal is : size: 3 by 1
-1.0000
-1.0000
-1.0000
The elements on main diagonal is : size: 4 by 1
4.0000
4.0000
4.0000
4.0000
The elements above main diagonal is : size: 3 by 1
-2.0000
-2.0000
-2.0000
The elements constant vector is : size: 4 by 1
3.0000
2.0000
2.0000
3.0000
Teh solution is (Forward Elimination and Backward Substitution) : size: 4 by 1
1.5610
1.6220
1.4634
1.1159
The elements below main diagonal is : size: 3 by 1
(-1.0000,0.0000)
(-1.0000,0.0000)
(-1.0000,0.0000)
The elements on main diagonal is : size: 4 by 1
(4.0000,0.0000)
(4.0000,0.0000)
(4.0000,0.0000)
(4.0000,0.0000)
The elements above main diagonal is : size: 3 by 1
(-2.0000,0.0000)
(-2.0000,0.0000)
(-2.0000,0.0000)
The elements constant vector is : size: 4 by 1
(0.0000,3.0000)
(0.0000,2.0000)
(0.0000,2.0000)
(0.0000,3.0000)
Teh solution is (Forward Elimination and Backward Substitution) : size: 4 by 1
(0.0000,1.5610)
(0.0000,1.6220)
(0.0000,1.4634)
(0.0000,1.1159)
Process returned 0 (0x0) execution time : 0.187 s
Press any key to continue.
扫一扫
2529
08-06
12-17
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
