头文件：

/** 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*//******************************************************************************                               ccholesky.h** Class template of complex matrix Cholesky factorization.** For a conjugate symmetric, positive definite matrix A, this class computes* the Cholesky factorization, i.e. it computes a lower triangular matrix L* such that A = L*L^H. If the matrix is not conjugate symmetric or positive* definite, the class computes only a partial decomposition. This can be* tested with the isSpd() flag.** Zhang Ming, 2010-12, Xi'an Jiaotong University.*****************************************************************************/#ifndef CCHOLESKY_H
#define CCHOLESKY_H#include <matrix.h>namespace splab
{template <typename Type>class CCholesky{public:CCholesky();~CCholesky();bool isSpd() const;void dec( const Matrix<Type> &A );Matrix<Type> getL() const;Vector<Type> solve( const Vector<Type> &b );Matrix<Type> solve( const Matrix<Type> &B );private:bool spd;Matrix<Type> L;};//    class CCholesky#include <ccholesky-impl.h>}
// namespace splab#endif
// CCHOLESKY_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*//******************************************************************************                             ccholesky-impl.h** Implementation for CCholesky (complex Cholesky) class.** Zhang Ming, 2010-12, Xi'an Jiaotong University.*****************************************************************************//*** constructor and destructor*/
template<typename Type>
CCholesky<Type>::CCholesky() : spd(true)
{
}template<typename Type>
CCholesky<Type>::~CCholesky()
{
}/*** return true, if original matrix is symmetric positive-definite.*/
template<typename Type>
inline bool CCholesky<Type>::isSpd() const
{return spd;
}/*** Constructs a lower triangular matrix L, such that L*L^H= 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 <typename Type>
void CCholesky<Type>::dec( const Matrix<Type> &A )
{int m = A.rows();int n = A.cols();spd = (m == n);if( !spd )return;L = Matrix<Type>(A.cols(),A.cols());// main loopfor( int j=0; j<A.rows(); ++j ){Type d = 0;spd = spd && (imag(A[j][j]) == 0);for( int k=0; k<j; ++k ){Type s = 0;for( int i=0; i<k; ++i )s += L[k][i] * conj(L[j][i]);L[j][k] = s = (A[j][k]-s) / L[k][k];d = d + s*conj(s);spd = spd && (A[k][j] == conj(A[j][k]));}d = A[j][j] - d;spd = spd && ( real(d) > 0 );L[j][j] = sqrt( real(d) > 0 ? d : 0 );for( int k=j+1; k<A.rows(); ++k )L[j][k] = 0;}
}/*** return the lower triangular factor, L, such that L*L'=A.*/
template<typename Type>
inline Matrix<Type> CCholesky<Type>::getL() const
{return L;
}/*** Solve a linear system A*x = b, using the previously computed* cholesky factorization of A: L*L'.*/
template <typename Type>
Vector<Type> CCholesky<Type>::solve( const Vector<Type> &b )
{int n = L.rows();if( b.dim() != n )return Vector<Type>();Vector<Type> x = b;// solve L*y = bfor( int k=0; k<n; ++k ){for( int i=0; i<k; ++i )x[k] -= x[i]*L[k][i];x[k] /= L[k][k];}// solve L^H*x = yfor( int k=n-1; k>=0; --k ){for( int i=k+1; i<n; ++i )x[k] -= x[i]*conj(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 <typename Type>
Matrix<Type> CCholesky<Type>::solve( const Matrix<Type> &B )
{int n = L.rows();if( B.rows() != n )return Matrix<Type>();Matrix<Type> X = B;int nx = B.cols();// solve L*Y = Bfor( int j=0; j<nx; ++j )for( int k=0; k<n; ++k ){for( int i=0; i<k; ++i )X[k][j] -= X[i][j]*L[k][i];X[k][j] /= L[k][k];}// solve L^H*x = yfor( int j=0; j<nx; ++j )for( int k=n-1; k>=0; --k ){for( int i=k+1; i<n; ++i )X[k][j] -= X[i][j]*conj(L[i][k]);X[k][j] /= L[k][k];}return X;
}

测试代码：

/******************************************************************************                             ccholesky_test.cpp** Comnplex Cholesky class testing.** Zhang Ming, 2010-12, Xi'an Jiaotong University.*****************************************************************************/#define BOUNDS_CHECK#include <iostream>
#include <iomanip>
#include <ccholesky.h>using namespace std;
using namespace splab;typedef double  Type;
const   int     N = 5;int main()
{cout << setiosflags(ios::fixed) << setprecision(3);Matrix<complex<Type> > A(N,N), L(N,N);Vector<complex<Type> > b(N);for( int i=1; i<N+1; ++i ){for( int j=1; j<N+1; ++j )if( i == j )A(i,i) = complex<Type>(N+i,0);elseif( i < j )A(i,j) = complex<Type>(Type(i),cos(Type(i)));elseA(i,j) = complex<Type>(Type(j),-cos(Type(j)));b(i) = i*(i+1)/2 + i*(N-i);}cout << "The original matrix A : " << A << endl;CCholesky<complex<Type> > 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^H is : " << A - L*trH(L) << endl;Vector<complex<Type> > x = cho.solve(b);cout << "The constant vector b : " << b << endl;cout << "The solution of Ax = b : " << x << endl;cout << "The solution of Ax - b : " << A*x-b << endl;Matrix<complex<Type> > IA = cho.solve(eye(N,complex<Type>(1,0)));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
(6.000,0.000)   (1.000,0.540)   (1.000,0.540)   (1.000,0.540)   (1.000,0.540)(1.000,-0.540)  (7.000,0.000)   (2.000,-0.416)  (2.000,-0.416)  (2.000,-0.416)(1.000,-0.540)  (2.000,0.416)   (8.000,0.000)   (3.000,-0.990)  (3.000,-0.990)(1.000,-0.540)  (2.000,0.416)   (3.000,0.990)   (9.000,0.000)   (4.000,-0.654)(1.000,-0.540)  (2.000,0.416)   (3.000,0.990)   (4.000,0.654)   (10.000,0.000)The lower triangular matrix L is : size: 5 by 5
(2.449,0.000)   (0.000,0.000)   (0.000,0.000)   (0.000,0.000)   (0.000,0.000)(0.408,-0.221)  (2.605,0.000)   (0.000,0.000)   (0.000,0.000)   (0.000,0.000)(0.408,-0.221)  (0.685,0.160)   (2.700,0.000)   (0.000,0.000)   (0.000,0.000)(0.408,-0.221)  (0.685,0.160)   (0.848,0.367)   (2.727,0.000)   (0.000,0.000)(0.408,-0.221)  (0.685,0.160)   (0.848,0.367)   (0.893,0.240)   (2.753,0.000)A - L*L^H 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,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,0.000)
(9.000,0.000)
(12.000,0.000)
(14.000,0.000)
(15.000,0.000)The solution of Ax = b : size: 5 by 1
(0.356,-0.310)
(0.562,0.207)
(0.746,0.284)
(0.843,-0.037)
(0.842,-0.214)The solution of Ax - b : size: 5 by 1
(0.000,-0.000)
(-0.000,0.000)
(-0.000,0.000)
(-0.000,0.000)
(0.000,0.000)The inverse matrix of A : size: 5 by 5
(0.177,0.000)   (-0.018,-0.007) (-0.014,-0.004) (-0.008,-0.006) (-0.005,-0.007)(-0.018,0.007)  (0.164,-0.000)  (-0.024,0.013)  (-0.020,0.003)  (-0.017,-0.002)(-0.014,0.004)  (-0.024,-0.013) (0.160,-0.000)  (-0.032,0.018)  (-0.029,0.008)(-0.008,0.006)  (-0.020,-0.003) (-0.032,-0.018) (0.150,0.000)   (-0.043,0.012)(-0.005,0.007)  (-0.017,0.002)  (-0.029,-0.008) (-0.043,-0.012) (0.132,0.000)The product of  A*inv(A) : size: 5 by 5
(1.000,0.000)   (-0.000,-0.000) (-0.000,-0.000) (-0.000,-0.000) (0.000,-0.000)(0.000,-0.000)  (1.000,0.000)   (-0.000,0.000)  (0.000,-0.000)  (0.000,-0.000)(0.000,-0.000)  (-0.000,-0.000) (1.000,0.000)   (-0.000,0.000)  (-0.000,-0.000)(0.000,-0.000)  (0.000,0.000)   (-0.000,-0.000) (1.000,0.000)   (0.000,0.000)(0.000,-0.000)  (0.000,0.000)   (-0.000,0.000)  (0.000,0.000)   (1.000,-0.000)Process returned 0 (0x0)   execution time : 0.125 s
Press any key to continue.