Commit 9886f9a7 authored by James Sutherland's avatar James Sutherland
Browse files

Merge branch '46-complex-eigen-values-with-Field-Matrices' into 'master'

Resolve "What to do with complex eigenvalues"

Closes #46

See merge request !36
parents abb34f07 433f92c2
Pipeline #597 failed with stages
in 17 minutes and 39 seconds
......@@ -416,7 +416,17 @@ public:
* are not supported at this time. In debug mode this will throw an error if
* a complex value is computed. See MatVecOps.h for implementation.
*/
MatUnaryOp< EigenDecomposition<FieldT> > eigen_values();
MatUnaryOp< EigenDecomposition<FieldT> > real_eigenvalues();
/** @brief Calculate eigen values of this matrix
*
* @param realVector vector to return real values into
* @param complexVector vector to return complex values into
*
* This operation calculates the eigenvalues of this matrix and returns the
* REAL and COMPLEX value. See MatVecOps.h for implementation.
*/
void eigenvalues(FieldVector<FieldT>& realVector,
FieldVector<FieldT>& complexVector);
/** @brief Direct linear solve using LU decomposition
*
......
This diff is collapsed.
......@@ -157,7 +157,7 @@ class PoolAutoPtr{
/**
* \brief put managed memory back into other pool on destruction
*/
inline ~PoolAutoPtr() { Pool<T>::put( activeIndex, ptr ); }
inline ~PoolAutoPtr() { if(ptr != NULL) Pool<T>::put( activeIndex, ptr ); }
private:
T* ptr;
......
......@@ -66,6 +66,50 @@ void generate_matrix_with_eigen( const typename FieldT::value_type * const eigen
}
}
template<typename FieldT>
void generate_matrix_with_complex_eigen( const typename FieldT::value_type * const eigenValues,
const size_t eigenSize,
FieldMatrix<FieldT>& M,
FieldMatrix<FieldT>& scratchM,
FieldVector<FieldT>& scratchV,
SpatFldPtr<FieldT>& scratch,
const short int deviceIndex )
{
assert( !scratch.isnull() );
assert( M.elements() % 2 == 0 ); //Can only generate matrices of 2n x 2n
assert( eigenSize == (M.elements()/2) ); //Complex comes in +- pairs
assert( M.elements() == scratchM.elements() && scratchM.elements() == scratchV.elements() );
const size_t elements = M.elements();
size_t m = elements;
size_t twiceN = elements;
//Assign eigen values on block diagonal top right and bottom left (block size 2)
m = elements;
twiceN = elements;
while( m-- ){
twiceN = elements;
while( twiceN-- ){
FieldT& res = M.at(m*elements + twiceN);
//Set top right block value to positive eigenvalue
if( m % 2 == 0 && twiceN == m + 1 )
{
res <<= eigenValues[(m/2)];
}
else if( m % 2 == 1 && twiceN == m - 1 )
{
res <<= -eigenValues[(m/2)];
}
else
{
res <<= 0;
}
}
}
}
//==============================================================================
template<typename FieldT>
......@@ -124,7 +168,7 @@ bool matrix_eigen_tests( FieldMatrix<FieldT>& A,
}
}
x = A.eigen_values();
x = A.real_eigenvalues();
//sort result
typename FieldT::value_type buffer[length];
......@@ -170,25 +214,133 @@ bool matrix_eigen_tests( FieldMatrix<FieldT>& A,
}
}
# ifndef NDEBUG
//Test rotation matrices for failure
if( A.elements() > 1 ){
srand(10);
TestHelper noisyDeep( noisyOn );
//Rotation matrix is very likely to have at least one eigen value complex
random_matrix(A, MatrixTypes::ROTATION, deviceIndex, scratch);
bool pass = false;
try{
x = A.eigen_values();
}
catch(...){
pass = true;
}
return noisy.ok();
}
template<typename FieldT>
bool matrix_complex_eigen_tests( FieldMatrix<FieldT>& A,
FieldMatrix<FieldT>& B,
FieldVector<FieldT>& x,
FieldVector<FieldT>& b,
FieldVector<FieldT>& c,
FieldT& compare,
SpatFldPtr<FieldT>& scratch,
const bool noisyOn,
const size_t length,
const MemoryWindow w,
const short int deviceIndex )
{
TestHelper noisy( noisyOn );
{
std::pair<std::pair<MatrixTypes_e, bool>, char const *> validMatrices[] = {
std::make_pair(std::make_pair(MatrixTypes::RANDOM , true), "Random"),
};
const size_t typeSize = array_size(validMatrices);
for( size_t typeI = 0; typeI < typeSize; typeI++ ){
TestHelper noisyType( noisyOn );
MatrixTypes_e matrixType = validMatrices[typeI].first.first;
bool generateEigen = validMatrices[typeI].first.second;
char const * description = validMatrices[typeI].second;
srand(10);
const size_t matricesToTest = 10;
for( size_t testMatrixIndex = 0; testMatrixIndex < matricesToTest; testMatrixIndex++ ){
TestHelper noisyDeep( noisyOn );
//Generate Matrix A and Eigen Values in b
if( generateEigen ){
const size_t eigenSize = length/2;
typename FieldT::value_type eigenVals[eigenSize];
for( size_t eIndex = 0; eIndex < eigenSize; eIndex++){
eigenVals[eIndex] = (static_cast<typename FieldT::value_type>(rand())/RAND_MAX)*MatrixRandomizer<FieldT>::get_magnitude();
# ifndef NDEBUG
for( size_t checkIndex = 0; checkIndex < eIndex; checkIndex++ ){
assert(eigenVals[eIndex] != eigenVals[checkIndex]);
}
# endif /* NDEBUG */
}
generate_matrix_with_complex_eigen(&eigenVals[0], eigenSize, A, B, b, scratch, deviceIndex);
for( size_t eIndex = 0; eIndex < eigenSize; eIndex++ ){
b.at(eIndex) <<= eigenVals[eIndex];
}
for( size_t eIndex = 0; eIndex < eigenSize; eIndex++ ){
b.at(eigenSize + eIndex) <<= -eigenVals[eIndex];
}
}
else{
//DEPRECIATED: only a diagonal matrix of complex values would have
//complex eigen values on the diagonal, which we cannot do using
//only a real matrix
random_matrix(A, matrixType, deviceIndex, scratch);
for( size_t eIndex = 0; eIndex < length; eIndex++ ){
b.at(eIndex) = A.at(eIndex*length+eIndex);
}
}
A.eigenvalues(x, c);
//sort real result
typename FieldT::value_type buffer[length];
x.set_device_as_active(CPU_INDEX);
for( size_t i = 0; i < w.glob_npts(); i++ ){
for( size_t j = 0; j < length; j++ ){
buffer[j] = x.at(j)[i];
}
std::sort(buffer, buffer+length);
for( size_t j = 0; j < length; j++){
x.at(j)[i] = buffer[j];
}
}
x.set_device_as_active(deviceIndex);
//sort complex result
c.set_device_as_active(CPU_INDEX);
for( size_t i = 0; i < w.glob_npts(); i++ ){
for( size_t j = 0; j < length; j++ ){
buffer[j] = c.at(j)[i];
}
std::sort(buffer, buffer+length);
for( size_t j = 0; j < length; j++){
c.at(j)[i] = buffer[j];
}
}
c.set_device_as_active(deviceIndex);
//sort eigen compare
b.set_device_as_active(CPU_INDEX);
for( size_t i = 0; i < w.glob_npts(); i++ ){
for( size_t j = 0; j < length; j++ ){
buffer[j] = b.at(j)[i];
}
std::sort(buffer, buffer+length);
for( size_t j = 0; j < length; j++ ){
b.at(j)[i] = buffer[j];
}
}
b.set_device_as_active(deviceIndex);
noisy( pass, "Expected failure of rotation matrix eigen values" );
//verify result
for( size_t i = 0; i < length; i++ ){
std::ostringstream msg;
msg << "correct sorted eigen val #" << i;
noisyDeep(field_equal(c(i), b(i), 1e-4), msg.str());
}
for( size_t i = 0; i < length; i++ ){
std::ostringstream msg;
msg << "correct real eigen val #" << i;
noisyDeep(field_equal(0, x(i), 1e-4), msg.str());
}
std::ostringstream msg;
msg << "eigen values of matrix #" << testMatrixIndex;
noisyType( noisyDeep.ok(), msg.str() );
}
std::ostringstream msg;
msg << "Eigen matrix of type " << description;
noisy( noisyType.ok(), msg.str() );
}
}
# endif // NDEBUG
return noisy.ok();
}
......@@ -678,7 +830,7 @@ bool drive_test( const IntVec dim, const int length, const bool verbose )
quietDevice( result, "pointwise matrix solve" );
}
// test eigenvalues randomly
// test only real eigenvalues randomly
{
const bool result = matrix_eigen_tests( A,
B,
......@@ -690,7 +842,23 @@ bool drive_test( const IntVec dim, const int length, const bool verbose )
length,
w,
deviceIndex );
quietDevice( result, "eigen values of matrix" );
quietDevice( result, "only real eigen values of matrix" );
}
// test only complex eigenvalues randomly
if( length % 2 == 0 ) {
const bool result = matrix_complex_eigen_tests( A,
B,
x,
b,
y,
compare,
scratch,
noisyOn,
length,
w,
deviceIndex );
quietDevice( result, "only complex eigen values of matrix" );
}
//print results of device run
......@@ -835,7 +1003,7 @@ double drive_perf( const IntVec listSize,
cudaDeviceSynchronize();
#endif
timer.reset();
x = A.eigen_values();
x = A.real_eigenvalues();
#ifdef __CUDACC__
cudaDeviceSynchronize();
#endif
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment