头文件：
/*
* 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
*/
/*****************************************************************************
* cevd.h
*
* Class template of eigenvalues decomposition for complex matrix.
*
* For a complex matrix A, we have A*V = V*D, where the eigenvalue matrix D
* is diagonal and the eigenvector matrix V is linear independence. That is,
* the diagonal values of D are the eigenvalues and the columns of V represent
* the corresponding eigenvectors of D. If A is Hermitian, then V is a unitary
* matrix, which means A = V*D*V', and eigenvalues are all real numbers.
*
* The matrix V may be badly conditioned, or even singular, so the validity
* of the equation A=V*D*inverse(V) depends upon the condition number of V.
*
* Zhang Ming, 2010-12, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef CEVD_H
#define CEVD_H
#include <evd.h>
namespace splab
{
template <typename Type>
class CEVD
{
public:
CEVD();
~CEVD();
// decomposition
void dec( const Matrix<complex<Type> > &A );
// the eigenvalues are real or complex
bool isHertimian() const;
// get eigenvectors and
Matrix<complex<Type> > getV() const;
Vector<complex<Type> > getD() const;
Vector<Type> getRD() const;
private:
bool hermitian;
// eigenvectors and eigenvalues
Matrix<complex<Type> > V;
Vector<complex<Type> > d;
Vector<Type> rd;
};
// class CEVD
#include <cevd-impl.h>
}
// namespace splab
#endif
// CEVD_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
*/
/*****************************************************************************
* cevd-impl.h
*
* Implementation for CEVD class.
*
* Zhang Ming, 2010-12, Xi'an Jiaotong University.
*****************************************************************************/
/**
* constructor and destructor
*/
template<typename Type>
CEVD<Type>::CEVD() : hermitian(true)
{
}
template<typename Type>
CEVD<Type>::~CEVD()
{
}
/**
* Check for symmetry, then construct the eigenvalue decomposition
*/
template <typename Type>
void CEVD<Type>::dec( const Matrix<complex<Type> > &A )
{
int N = A.cols();
assert( A.rows() == N );
V = Matrix<complex<Type> >(N,N);
Matrix<Type> S(2*N,2*N);
for( int i=0; i<N; ++i )
for( int j=0; j<N; ++j )
{
S[i][j] = A[i][j].real();
S[i][j+N] = -A[i][j].imag();
S[i+N][j] = -S[i][j+N];
S[i+N][j+N] = S[i][j];
}
EVD<Type> eig;
eig.dec(S);
if( eig.isSymmetric() )
{
hermitian = true;
rd = Vector<Type>(N);
Matrix<Type> RV = eig.getV();
Vector<Type> Rd = eig.getD();
for( int i=0; i<N; ++i )
rd[i] = Rd[2*i];
for( int j=0; j<N; ++j )
{
int j2 = 2*j;
for( int i=0; i<N; ++i )
V[i][j] = complex<Type>( RV[i][j2], RV[i+N][j2] );
}
}
else
{
hermitian = false;
d = Vector<complex<Type> >(N);
Matrix<complex<Type> > cV = eig.getCV();
Vector<complex<Type> > cd = eig.getCD();
for( int i=0; i<N; ++i )
d[i] = cd[2*i];
for( int j=0; j<N; ++j )
{
int j2 = 2*j;
for( int i=0; i<N; ++i )
V[i][j] = complex<Type>( cV[i][j2].real()-cV[i+N][j2].imag(),
cV[i][j2].imag()+cV[i+N][j2].real() );
}
}
}
/**
* If the matrix is Hermitian, then return true.
*/
template <typename Type>
bool CEVD<Type>::isHertimian() const
{
return hermitian;
}
/**
* Return the COMPLEX eigenvector matrix
*/
template <typename Type>
inline Matrix<complex<Type> > CEVD<Type>::getV() const
{
return V;
}
/**
* Return the complex eigenvalues vector.
*/
template <typename Type>
inline Vector<complex<Type> > CEVD<Type>::getD() const
{
return d;
}
/**
* Return the real eigenvalues vector.
*/
template <typename Type>
inline Vector<Type> CEVD<Type>::getRD() const
{
return rd;
}测试代码：
/*****************************************************************************
* cevd_test.cpp
*
* CEVD class testing.
*
* Zhang Ming, 2010-12, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <iomanip>
#include <cevd.h>
using namespace std;
using namespace splab;
typedef double Type;
const int N = 4;
int main()
{
Matrix<Type> B(N,N);
B[0][0] = 3.0; B[0][1] = -2.0; B[0][2] = -0.9; B[0][3] = 0.0;
B[1][0] = -2.0; B[1][1] = 4.0; B[1][2] = 1.0; B[1][3] = 0.0;
B[2][0] = 0.0; B[2][1] = 0.0; B[2][2] = -1.0; B[2][3] = 0.0;
B[3][0] = -0.5; B[3][1] = -0.5; B[3][2] = 0.1; B[3][3] = 1.0;
Matrix<complex<Type> > A = complexMatrix( B, elemMult(B,B) );
// A = multTr(A,A);
cout << setiosflags(ios::fixed) << setprecision(2);
cout << "The original complex matrix A : " << A << endl;
CEVD<Type> eig;
eig.dec(A);
if( eig.isHertimian() )
{
Matrix<complex<Type> > V = eig.getV();
Vector<Type> D = eig.getRD();
Matrix<complex<Type> > DM = diag( complexVector(D) );
cout << "The eigenvectors matrix V is : " << V << endl;
cout << "The eigenvalue D is : " << diag(D) << endl;
cout << "The V'*V : " << trMult(V,V) << endl;
cout << "The A*V - V*D : " << A*V - V*DM << endl;
}
else
{
Matrix<complex<Type> > V = eig.getV();
Vector<complex<Type> > D = eig.getD();
Matrix<complex<Type> > DM = diag( D );
cout << "The complex eigenvectors matrix V : " << V << endl;
cout << "The complex eigenvalue D : " << DM << endl;
cout << "The A*V - V*D : " << A*V - V*DM << endl;
}
return 0;
}运行结果：
The original complex matrix A : size: 4 by 4
(3.00,9.00) (-2.00,4.00) (-0.90,0.81) (0.00,0.00)
(-2.00,4.00) (4.00,16.00) (1.00,1.00) (0.00,0.00)
(0.00,0.00) (0.00,0.00) (-1.00,1.00) (0.00,0.00)
(-0.50,0.25) (-0.50,0.25) (0.10,0.01) (1.00,1.00)
The complex eigenvectors matrix V : size: 4 by 4
(-0.54,-0.75) (1.50,-1.10) (0.00,0.00) (-0.19,-0.09)
(-1.49,-0.96) (-0.94,0.24) (-0.00,-0.00) (0.05,0.23)
(0.00,0.00) (0.00,0.00) (-0.00,0.00) (0.54,-1.92)
(0.03,-0.08) (0.06,0.05) (-1.11,1.67) (-0.05,0.15)
The complex eigenvalue D : size: 4 by 4
(2.26,17.55) (0.00,0.00) (0.00,0.00) (0.00,0.00)
(0.00,0.00) (4.74,7.45) (0.00,0.00) (0.00,0.00)
(0.00,0.00) (0.00,0.00) (1.00,1.00) (0.00,0.00)
(0.00,0.00) (0.00,0.00) (0.00,0.00) (-1.00,1.00)
The A*V - V*D : size: 4 by 4
(-0.00,0.00) (0.00,0.00) (0.00,0.00) (0.00,0.00)
(-0.00,0.00) (-0.00,0.00) (0.00,-0.00) (0.00,0.00)
(0.00,-0.00) (-0.00,-0.00) (0.00,-0.00) (-0.00,-0.00)
(-0.00,0.00) (0.00,0.00) (-0.00,0.00) (-0.00,0.00)
Process returned 0 (0x0) execution time : 0.094 s
Press any key to continue.
5052
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
