头文件：
/*
* 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
*/
/*****************************************************************************
* cholesky.h
*
* Class template of Cholesky decomposition.
*
* For a symmetric, positive definite matrix A, this function computes the
* Cholesky factorization, i.e. it computes a lower triangular matrix L such
* that A = L*L’. If the matrix is not symmetric or positive definite, the
* function computes only a partial decomposition. This can be tested with
* the isSpd() flag.
*
* This class also supports factorization of complex matrix by specializing
* some member functions.
*
* Adapted form Template Numerical Toolkit.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi’an Jiaotong University.
*****************************************************************************/
#ifndef CHOLESKY_H
#define CHOLESKY_H
#include
namespace splab
{
template
class Cholesky
{
public:
Cholesky();
~Cholesky();
bool isSpd() const;
void dec( const Matrix &A );
Matrix getL() const;
Vector solve( const Vector &b );
Matrix solve( const Matrix &B );
private:
bool spd;
Matrix L;
};
//class Cholesky
#include
}
// namespace splab
#endif
// CHOLESKY_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
*/
/*****************************************************************************
* cholesky-impl.h
*
* Implementation for Cholesky class.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi’an Jiaotong University.
*****************************************************************************/
/**
* constructor and destructor
*/
template
Cholesky::Cholesky() : spd(true)
{
}
template
Cholesky::~Cholesky()
{
}
/**
* return true, if original matrix is symmetric positive-definite.
*/
template
inline bool Cholesky::isSpd() const
{
return spd;
}
/**
* Constructs a lower triangular matrix L, such that L*L’= A.
* If A is not symmetric positive-definite (SPD), only a partial
* factorization is performed. If isspd() evalutate true then
* the factorizaiton was successful.
*/
template
void Cholesky::dec( const Matrix &A )
{
int m = A.rows();
int n = A.cols();
spd = (m == n);
if( !spd )
return;
L = Matrix(n,n);
// main loop.
for( int j=0; j
{
Type d = 0;
for( int k=0; k
{
Type s = 0;
for( int i=0; i
s = L[k][i]*L[j][i];
L[j][k] = s = (A[j][k]-s) / L[k][k];
d = d s*s;
spd = spd && (A[k][j] == A[j][k]);
}
d = A[j][j] – d;
spd = spd && ( d > 0 );
L[j][j] = sqrt( d > 0 ? d : 0 );
for( int k=j 1; k
L[j][k] = 0;
}
}
/**
* return the lower triangular factor, L, such that L*L’=A.
*/
template
inline Matrix Cholesky::getL() const
{
return L;
}
/**
* Solve a linear system A*x = b, using the previously computed
* cholesky factorization of A: L*L’.
*/
template
Vector Cholesky::solve( const Vector &b )
{
int n = L.rows();
if( b.dim() != n )
return Vector();
Vector x = b;
// solve L*y = b
for( int k=0; k
{
for( int i=0; i
x[k] -= x[i]*L[k][i];
x[k] /= L[k][k];
}
// solve L’*x = y
for( int k=n-1; k>=0; –k )
{
for( int i=k 1; i
x[k] -= x[i]*L[i][k];
x[k] /= L[k][k];
}
return x;
}
/**
* Solve a linear system A*X = B, using the previously computed
* cholesky factorization of A: L*L’.
*/
template
Matrix Cholesky::solve( const Matrix &B )
{
int n = L.rows();
if( B.rows() != n )
return Matrix();
Matrix X = B;
int nx = B.cols();
// solve L*Y = B
for( int j=0; j
for( int k=0; k
{
for( int i=0; i
X[k][j] -= X[i][j]*L[k][i];
X[k][j] /= L[k][k];
}
// solve L’*X = Y
for( int j=0; j
for( int k=n-1; k>=0; –k )
{
for( int i=k 1; i
X[k][j] -= X[i][j]*L[i][k];
X[k][j] /= L[k][k];
}
return X;
}
/**
* Main loop of specialized member function. This macro definition is
* aimed at avoiding code duplication.
*/
#define MAINLOOP
{
int m = A.rows();
int n = A.cols();
spd = (m == n);
if( !spd )
return;
for( int j=0; j
{
spd = spd && (imag(A[j][j]) == 0);
d = 0;
for( int k=0; k
{
s = 0;
for( int i=0; i
s = L[k][i] * conj(L[j][i]);
L[j][k] = s = (A[j][k]-s) / L[k][k];
d = d norm(s);
spd = spd && (A[k][j] == conj(A[j][k]));
}
d = real(A[j][j]) – d;
spd = spd && ( d > 0 );
L[j][j] = sqrt( d > 0 ? d : 0 );
for( int k=j 1; k
L[j][k] = 0;
}
}
/**
* Solving process of specialized member function. This macro definition is
* aimed at avoiding code duplication.
*/
#define SOLVE1
{
for( int k=0; k
{
for( int i=0; i
x[k] -= x[i]*L[k][i];
x[k] /= L[k][k];
}
for( int k=n-1; k>=0; –k )
{
for( int i=k 1; i
x[k] -= x[i]*conj(L[i][k]);
x[k] /= L[k][k];
}
return x;
}
/**
* Solving process of specialized member function. This macro definition is
* aimed at avoiding code duplication.
*/
#define SOLVE2
{
int nx = B.cols();
for( int j=0; j
for( int k=0; k
{
for( int i=0; i
X[k][j] -= X[i][j]*L[k][i];
X[k][j] /= L[k][k];
}
for( int j=0; j
for( int k=n-1; k>=0; –k )
{
for( int i=k 1; i
X[k][j] -= X[i][j]*conj(L[i][k]);
X[k][j] /= L[k][k];
}
return X;
}
/**
* Specializing for “dec” member function.
*/
template <>
void Cholesky >::dec( const Matrix > &A )
{
float d;
complex s;
L = Matrix >(A.cols(),A.cols());
MAINLOOP;
}
template <>
void Cholesky >::dec( const Matrix > &A )
{
double d;
complex s;
L = Matrix >(A.cols(),A.cols());
MAINLOOP;
}
template <>
void Cholesky >::dec( const Matrix > &A )
{
long double d;
complex s;
L = Matrix >(A.cols(),A.cols());
MAINLOOP;
}
/**
* Specializing for “solve” member function.
*/
template <>
Vector > Cholesky >::solve( const Vector > &b )
{
int n = L.rows();
if( b.dim() != n )
return Vector >();
Vector > x = b;
SOLVE1
}
template <>
Vector > Cholesky >::solve( const Vector > &b )
{
int n = L.rows();
if( b.dim() != n )
return Vector >();
Vector > x = b;
SOLVE1
}
template <>
Vector > Cholesky >::solve( const Vector > &b )
{
int n = L.rows();
if( b.dim() != n )
return Vector >();
Vector > x = b;
SOLVE1
}
/**
* Specializing for “solve” member function.
*/
template <>
Matrix > Cholesky >::solve( const Matrix > &B )
{
int n = L.rows();
if( B.rows() != n )
return Matrix >();
Matrix > X = B;
SOLVE2
}
template <>
Matrix > Cholesky >::solve( const Matrix > &B )
{
int n = L.rows();
if( B.rows() != n )
return Matrix >();
Matrix > X = B;
SOLVE2
}
template <>
Matrix > Cholesky >::solve( const Matrix > &B )
{
int n = L.rows();
if( B.rows() != n )
return Matrix >();
Matrix > X = B;
SOLVE2
}
测试代码：
/*****************************************************************************
* cholesky_test.cpp
*
* Cholesky class testing.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi’an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include
#include
#include
using namespace std;
using namespace splab;
typedef double Type;
const int N = 5;
int main()
{
Matrix A(N,N), L(N,N);
Vector b(N);
for( int i=1; i
{
for( int j=1; j
if( i == j )
A(i,i) = i;
else
if( i < j )
A(i,j) = i;
else
A(i,j) = j;
b(i) = i*(i 1)/2.0 i*(N-i);
}
cout << setiosflags(ios::fixed) << setprecision(3);
cout << “The original matrix A : ” << A << endl;
Cholesky cho;
cho.dec(A);
if( !cho.isSpd() )
cout << “Factorization was not complete.” << endl;
else
{
L = cho.getL();
cout << “The lower triangular matrix L is : ” << L << endl;
cout << “A – L*L^T is : ” << A – L*trT(L) << endl;
Vector x = cho.solve(b);
cout << “The constant vector b : ” << b << endl;
cout << “The solution of Ax = b : ” << x << endl;
cout << “The Ax – b : ” << A*x-b << endl;
Matrix IA = cho.solve(eye(N,Type(1)));
cout << “The inverse matrix of A : ” << IA << endl;
cout << “The product of A*inv(A) : ” << A*IA << endl;
}
return 0;
}
运行结果：
The original matrix A : size: 5 by 5
1.000 1.000 1.000 1.000 1.000
1.000 2.000 2.000 2.000 2.000
1.000 2.000 3.000 3.000 3.000
1.000 2.000 3.000 4.000 4.000
1.000 2.000 3.000 4.000 5.000
The lower triangular matrix L is : size: 5 by 5
1.000 0.000 0.000 0.000 0.000
1.000 1.000 0.000 0.000 0.000
1.000 1.000 1.000 0.000 0.000
1.000 1.000 1.000 1.000 0.000
1.000 1.000 1.000 1.000 1.000
A – L*L^T is : size: 5 by 5
0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000 0.000
The constant vector b : size: 5 by 1
5.000
9.000
12.000
14.000
15.000
The solution of Ax = b : size: 5 by 1
1.000
1.000
1.000
1.000
1.000
The Ax – b : size: 5 by 1
0.000
0.000
0.000
0.000
0.000
The inverse matrix of A : size: 5 by 5
2.000 -1.000 0.000 0.000 0.000
-1.000 2.000 -1.000 0.000 0.000
0.000 -1.000 2.000 -1.000 0.000
0.000 0.000 -1.000 2.000 -1.000
0.000 0.000 0.000 -1.000 1.000
The product of A*inv(A) : size: 5 by 5
1.000 0.000 0.000 0.000 0.000
0.000 1.000 0.000 0.000 0.000
0.000 0.000 1.000 0.000 0.000
0.000 0.000 0.000 1.000 0.000
0.000 0.000 0.000 0.000 1.000
Process returned 0 (0x0) execution time : 0.109 s
Press any key to continue.
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。
发布者：全栈程序员-用户IM，转载请注明出处：https://javaforall.cn/219123.html原文链接：https://javaforall.cn
