libMesh
dense_matrix_test.C
Go to the documentation of this file.
1 // Ignore unused parameter warnings coming from cppunit headers
2 #include <libmesh/ignore_warnings.h>
3 #include <cppunit/extensions/HelperMacros.h>
4 #include <cppunit/TestCase.h>
5 #include <libmesh/restore_warnings.h>
6 
7 // libmesh includes
8 #include <libmesh/dense_matrix.h>
9 #include <libmesh/dense_vector.h>
10 
11 #ifdef LIBMESH_HAVE_PETSC
12 #include "libmesh/petsc_macro.h"
13 #endif
14 
15 // THE CPPUNIT_TEST_SUITE_END macro expands to code that involves
16 // std::auto_ptr, which in turn produces -Wdeprecated-declarations
17 // warnings. These can be ignored in GCC as long as we wrap the
18 // offending code in appropriate pragmas. We can't get away with a
19 // single ignore_warnings.h inclusion at the beginning of this file,
20 // since the libmesh headers pull in a restore_warnings.h at some
21 // point. We also don't bother restoring warnings at the end of this
22 // file since it's not a header.
23 #include <libmesh/ignore_warnings.h>
24 
25 using namespace libMesh;
26 
27 class DenseMatrixTest : public CppUnit::TestCase
28 {
29 public:
30  void setUp() {}
31 
32  void tearDown() {}
33 
34  CPPUNIT_TEST_SUITE(DenseMatrixTest);
35 
36  CPPUNIT_TEST(testSVD);
37  CPPUNIT_TEST(testEVDreal);
38  CPPUNIT_TEST(testEVDcomplex);
39  CPPUNIT_TEST(testComplexSVD);
40 
41  CPPUNIT_TEST_SUITE_END();
42 
43 
44 private:
45 
46  void testSVD()
47  {
48  DenseMatrix<Number> U, VT;
49  DenseVector<Real> sigma;
51 
52  A.resize(3, 2);
53  A(0,0) = 1.0; A(0,1) = 2.0;
54  A(1,0) = 3.0; A(1,1) = 4.0;
55  A(2,0) = 5.0; A(2,1) = 6.0;
56 
57  A.svd(sigma, U, VT);
58 
59  // Solution for this case is (verified with numpy)
60  DenseMatrix<Number> true_U(3,2), true_VT(2,2);
61  DenseVector<Real> true_sigma(2);
62  true_U(0,0) = -2.298476964000715e-01; true_U(0,1) = 8.834610176985250e-01;
63  true_U(1,0) = -5.247448187602936e-01; true_U(1,1) = 2.407824921325463e-01;
64  true_U(2,0) = -8.196419411205157e-01; true_U(2,1) = -4.018960334334318e-01;
65 
66  true_VT(0,0) = -6.196294838293402e-01; true_VT(0,1) = -7.848944532670524e-01;
67  true_VT(1,0) = -7.848944532670524e-01; true_VT(1,1) = 6.196294838293400e-01;
68 
69  true_sigma(0) = 9.525518091565109e+00;
70  true_sigma(1) = 5.143005806586446e-01;
71 
72  for (unsigned i=0; i<U.m(); ++i)
73  for (unsigned j=0; j<U.n(); ++j)
74  CPPUNIT_ASSERT_DOUBLES_EQUAL( libmesh_real(U(i,j)), libmesh_real(true_U(i,j)), TOLERANCE*TOLERANCE);
75 
76  for (unsigned i=0; i<VT.m(); ++i)
77  for (unsigned j=0; j<VT.n(); ++j)
78  CPPUNIT_ASSERT_DOUBLES_EQUAL( libmesh_real(VT(i,j)), libmesh_real(true_VT(i,j)), TOLERANCE*TOLERANCE);
79 
80  for (unsigned i=0; i<sigma.size(); ++i)
81  CPPUNIT_ASSERT_DOUBLES_EQUAL(sigma(i), true_sigma(i), TOLERANCE*TOLERANCE);
82  }
83 
84  // This function is called by testEVD for different matrices. The
85  // Lapack results are compared to known eigenvalue real and
86  // imaginary parts for the matrix in question, which must also be
87  // passed in by non-const value, since this routine will sort them
88  // in-place.
90  std::vector<Real> true_lambda_real,
91  std::vector<Real> true_lambda_imag)
92  {
93  // Note: see bottom of this file, we only do this test if PETSc is
94  // available, but this function currently only exists if we're
95  // using real numbers.
96 #ifdef LIBMESH_USE_REAL_NUMBERS
97  // Let's compute the eigenvalues on a copy of A, so that we can
98  // use the original to check the computation.
99  DenseMatrix<Real> A_copy = A;
100 
101  DenseVector<Real> lambda_real, lambda_imag;
102  DenseMatrix<Real> VR; // right eigenvectors
103  DenseMatrix<Real> VL; // left eigenvectors
104  A_copy.evd_left_and_right(lambda_real, lambda_imag, VL, VR);
105 
106  // The matrix is square and of size N x N.
107  const unsigned N = A.m();
108 
109  // Verify left eigen-values.
110  // Test that the right eigenvalues are self-consistent by computing
111  // u_j**H * A = lambda_j * u_j**H
112  // Note that we have to handle real and complex eigenvalues
113  // differently, since complex eigenvectors share their storage.
114  for (unsigned eigenval=0; eigenval<N; ++eigenval)
115  {
116  // Only check real eigenvalues
117  if (std::abs(lambda_imag(eigenval)) < TOLERANCE*TOLERANCE)
118  {
119  // remove print libMesh::out << "Checking eigenvalue: " << eigenval << std::endl;
120  DenseVector<Real> lhs(N), rhs(N);
121  for (unsigned i=0; i<N; ++i)
122  {
123  rhs(i) = lambda_real(eigenval) * VL(i, eigenval);
124  for (unsigned j=0; j<N; ++j)
125  lhs(i) += A(j, i) * VL(j, eigenval); // Note: A(j,i)
126  }
127 
128  // Subtract and assert that the norm of the difference is
129  // below some tolerance.
130  lhs -= rhs;
131  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(TOLERANCE)*TOLERANCE);
132  }
133  else
134  {
135  // This is a complex eigenvalue, so:
136  // a.) It occurs in a complex-conjugate pair
137  // b.) the real part of the eigenvector is stored is VL(:,eigenval)
138  // c.) the imag part of the eigenvector is stored in VL(:,eigenval+1)
139  //
140  // Equating the real and imaginary parts of Ax=lambda*x leads to two sets
141  // of relations that must hold:
142  // 1.) A^T x_r = lambda_r*x_r + lambda_i*x_i
143  // 2.) A^T x_i = -lambda_i*x_r + lambda_r*x_i
144  // which we can verify.
145 
146  // 1.)
147  DenseVector<Real> lhs(N), rhs(N);
148  for (unsigned i=0; i<N; ++i)
149  {
150  rhs(i) = lambda_real(eigenval) * VL(i, eigenval) + lambda_imag(eigenval) * VL(i, eigenval+1);
151  for (unsigned j=0; j<N; ++j)
152  lhs(i) += A(j, i) * VL(j, eigenval); // Note: A(j,i)
153  }
154 
155  lhs -= rhs;
156  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(TOLERANCE)*TOLERANCE);
157 
158  // libMesh::out << "lhs=" << std::endl;
159  // lhs.print_scientific(libMesh::out, /*precision=*/15);
160  //
161  // libMesh::out << "rhs=" << std::endl;
162  // rhs.print_scientific(libMesh::out, /*precision=*/15);
163 
164  // 2.)
165  lhs.zero();
166  rhs.zero();
167  for (unsigned i=0; i<N; ++i)
168  {
169  rhs(i) = -lambda_imag(eigenval) * VL(i, eigenval) + lambda_real(eigenval) * VL(i, eigenval+1);
170  for (unsigned j=0; j<N; ++j)
171  lhs(i) += A(j, i) * VL(j, eigenval+1); // Note: A(j,i)
172  }
173 
174  lhs -= rhs;
175  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(TOLERANCE)*TOLERANCE);
176 
177  // libMesh::out << "lhs=" << std::endl;
178  // lhs.print_scientific(libMesh::out, /*precision=*/15);
179  //
180  // libMesh::out << "rhs=" << std::endl;
181  // rhs.print_scientific(libMesh::out, /*precision=*/15);
182 
183  // We'll skip the second member of the complex conjugate
184  // pair. If the first one worked, the second one should
185  // as well...
186  eigenval += 1;
187  }
188  }
189 
190  // Verify right eigen-values.
191  // Test that the right eigenvalues are self-consistent by computing
192  // A * v_j - lambda_j * v_j
193  // Note that we have to handle real and complex eigenvalues
194  // differently, since complex eigenvectors share their storage.
195  for (unsigned eigenval=0; eigenval<N; ++eigenval)
196  {
197  // Only check real eigenvalues
198  if (std::abs(lambda_imag(eigenval)) < TOLERANCE*TOLERANCE)
199  {
200  // remove print libMesh::out << "Checking eigenvalue: " << eigenval << std::endl;
201  DenseVector<Real> lhs(N), rhs(N);
202  for (unsigned i=0; i<N; ++i)
203  {
204  rhs(i) = lambda_real(eigenval) * VR(i, eigenval);
205  for (unsigned j=0; j<N; ++j)
206  lhs(i) += A(i, j) * VR(j, eigenval);
207  }
208 
209  lhs -= rhs;
210  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(TOLERANCE)*TOLERANCE);
211  }
212  else
213  {
214  // This is a complex eigenvalue, so:
215  // a.) It occurs in a complex-conjugate pair
216  // b.) the real part of the eigenvector is stored is VR(:,eigenval)
217  // c.) the imag part of the eigenvector is stored in VR(:,eigenval+1)
218  //
219  // Equating the real and imaginary parts of Ax=lambda*x leads to two sets
220  // of relations that must hold:
221  // 1.) Ax_r = lambda_r*x_r - lambda_i*x_i
222  // 2.) Ax_i = lambda_i*x_r + lambda_r*x_i
223  // which we can verify.
224 
225  // 1.)
226  DenseVector<Real> lhs(N), rhs(N);
227  for (unsigned i=0; i<N; ++i)
228  {
229  rhs(i) = lambda_real(eigenval) * VR(i, eigenval) - lambda_imag(eigenval) * VR(i, eigenval+1);
230  for (unsigned j=0; j<N; ++j)
231  lhs(i) += A(i, j) * VR(j, eigenval);
232  }
233 
234  lhs -= rhs;
235  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(TOLERANCE)*TOLERANCE);
236 
237  // 2.)
238  lhs.zero();
239  rhs.zero();
240  for (unsigned i=0; i<N; ++i)
241  {
242  rhs(i) = lambda_imag(eigenval) * VR(i, eigenval) + lambda_real(eigenval) * VR(i, eigenval+1);
243  for (unsigned j=0; j<N; ++j)
244  lhs(i) += A(i, j) * VR(j, eigenval+1);
245  }
246 
247  lhs -= rhs;
248  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(TOLERANCE)*TOLERANCE);
249 
250  // We'll skip the second member of the complex conjugate
251  // pair. If the first one worked, the second one should
252  // as well...
253  eigenval += 1;
254  }
255  }
256 
257  // Sort the results from Lapack *individually*.
258  std::sort(lambda_real.get_values().begin(), lambda_real.get_values().end());
259  std::sort(lambda_imag.get_values().begin(), lambda_imag.get_values().end());
260 
261  // Sort the true eigenvalues *individually*.
262  std::sort(true_lambda_real.begin(), true_lambda_real.end());
263  std::sort(true_lambda_imag.begin(), true_lambda_imag.end());
264 
265  // Compare the individually-sorted values.
266  for (unsigned i=0; i<lambda_real.size(); ++i)
267  {
268  // Note: I initially verified the results with TOLERANCE**2,
269  // but that turned out to be just a bit too tight for some of
270  // the test problems. I'm not sure what controls the accuracy
271  // of the eigenvalue computation in LAPACK, there is no way to
272  // set a tolerance in the LAPACKgeev_ interface.
273  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/true_lambda_real[i], /*actual=*/lambda_real(i), std::sqrt(TOLERANCE)*TOLERANCE);
274  CPPUNIT_ASSERT_DOUBLES_EQUAL(/*expected=*/true_lambda_imag[i], /*actual=*/lambda_imag(i), std::sqrt(TOLERANCE)*TOLERANCE);
275  }
276 #endif
277  }
278 
279  void testEVDreal()
280  {
281  // This is an example from Matlab's gallery(3) which is a
282  // non-symmetric 3x3 matrix with eigen values lambda = 1, 2, 3.
283  DenseMatrix<Real> A(3, 3);
284  A(0,0) = -149; A(0,1) = -50; A(0,2) = -154;
285  A(1,0) = 537; A(1,1) = 180; A(1,2) = 546;
286  A(2,0) = -27; A(2,1) = -9; A(2,2) = -25;
287 
288  std::vector<Real> true_lambda_real(3);
289  true_lambda_real[0] = 1.;
290  true_lambda_real[1] = 2.;
291  true_lambda_real[2] = 3.;
292  std::vector<Real> true_lambda_imag(3); // all zero
293 
294  // call helper function to compute and verify results.
295  testEVD_helper(A, true_lambda_real, true_lambda_imag);
296  }
297 
299  {
300  // This test is also from a Matlab example, and has complex eigenvalues.
301  // http://www.mathworks.com/help/matlab/math/eigenvalues.html?s_tid=gn_loc_drop
302  DenseMatrix<Real> A(3, 3);
303  A(0,0) = 0; A(0,1) = -6; A(0,2) = -1;
304  A(1,0) = 6; A(1,1) = 2; A(1,2) = -16;
305  A(2,0) = -5; A(2,1) = 20; A(2,2) = -10;
306 
307  std::vector<Real> true_lambda_real(3);
308  true_lambda_real[0] = -3.070950351248293;
309  true_lambda_real[1] = -2.464524824375853;
310  true_lambda_real[2] = -2.464524824375853;
311  std::vector<Real> true_lambda_imag(3);
312  true_lambda_imag[0] = 0.;
313  true_lambda_imag[1] = 17.60083096447099;
314  true_lambda_imag[2] = -17.60083096447099;
315 
316  // call helper function to compute and verify results
317  testEVD_helper(A, true_lambda_real, true_lambda_imag);
318  }
319 
321  {
322 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
323  DenseMatrix<Complex> A(3,3);
324 
325  A(0,0) = Complex(2.18904,4.44523e-18); A(0,1) = Complex(-3.20491,-0.136699); A(0,2) = Complex(0.716316,-0.964802);
326  A(1,0) = Complex(-3.20491,0.136699); A(1,1) = Complex(4.70076,-3.25261e-18); A(1,2) = Complex(-0.98849,1.45727);
327  A(2,0) = Complex(0.716316,0.964802); A(2,1) = Complex(-0.98849,-1.45727); A(2,2) = Complex(0.659629,-4.01155e-18);
328 
329  DenseVector<Real> sigma;
330  A.svd(sigma);
331 
332  DenseVector<Real> true_sigma(3);
333  true_sigma(0) = 7.54942516052;
334  true_sigma(1) = 3.17479511368e-06;
335  true_sigma(2) = 6.64680908281e-07;
336 
337  for (unsigned i=0; i<sigma.size(); ++i)
338  CPPUNIT_ASSERT_DOUBLES_EQUAL(sigma(i), true_sigma(i), 1.e-10);
339 #endif
340  }
341 
342 };
343 
344 // Only run the test if we expect it can actually work!
345 #ifdef LIBMESH_HAVE_PETSC
346 #if !PETSC_VERSION_LESS_THAN(3,1,0)
348 #endif
349 #endif
T libmesh_real(T a)
double abs(double a)
unsigned int n() const
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
std::vector< T > & get_values()
Definition: dense_vector.h:251
void evd_left_and_right(DenseVector< T > &lambda_real, DenseVector< T > &lambda_imag, DenseMatrix< T > &VL, DenseMatrix< T > &VR)
Compute the eigenvalues (both real and imaginary parts) as well as the left and right eigenvectors of...
Definition: dense_matrix.C:838
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
unsigned int m() const
void svd(DenseVector< Real > &sigma)
Compute the singular value decomposition of the matrix.
Definition: dense_matrix.C:777
CPPUNIT_TEST_SUITE_REGISTRATION(DenseMatrixTest)
std::complex< Real > Complex
void testEVD_helper(DenseMatrix< Real > &A, std::vector< Real > true_lambda_real, std::vector< Real > true_lambda_imag)
static PetscErrorCode Mat * A
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
Defines a dense vector for use in Finite Element-type computations.