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

- 日理万妓 2023-01-15 03:39 126阅读 0赞

头文件：

/\*

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

实现文件：

/\*

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

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 )