cholesky分解java代码,实数矩阵Cholesky分解算法的C++实现

头文件：

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang),

*

* 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),

*

* 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.