libMesh
Public Member Functions | Private Member Functions | List of all members
DenseMatrixTest Class Reference
Inheritance diagram for DenseMatrixTest:
[legend]

Public Member Functions

void setUp ()
 
void tearDown ()
 
 LIBMESH_CPPUNIT_TEST_SUITE (DenseMatrixTest)
 
 CPPUNIT_TEST (testOuterProduct)
 
 CPPUNIT_TEST (testSVD)
 
 CPPUNIT_TEST (testEVDreal)
 
 CPPUNIT_TEST (testEVDcomplex)
 
 CPPUNIT_TEST (testComplexSVD)
 
 CPPUNIT_TEST (testSubMatrix)
 
 CPPUNIT_TEST_SUITE_END ()
 

Private Member Functions

void testOuterProduct ()
 
void testSVD ()
 
void testEVD_helper (DenseMatrix< Real > &A, std::vector< Real > true_lambda_real, std::vector< Real > true_lambda_imag, Real tol)
 
void testEVDreal ()
 
void testEVDcomplex ()
 
void testComplexSVD ()
 
void testSubMatrix ()
 

Detailed Description

Definition at line 14 of file dense_matrix_test.C.

Member Function Documentation

◆ CPPUNIT_TEST() [1/6]

DenseMatrixTest::CPPUNIT_TEST ( testOuterProduct  )

◆ CPPUNIT_TEST() [2/6]

DenseMatrixTest::CPPUNIT_TEST ( testSVD  )

◆ CPPUNIT_TEST() [3/6]

DenseMatrixTest::CPPUNIT_TEST ( testEVDreal  )

◆ CPPUNIT_TEST() [4/6]

DenseMatrixTest::CPPUNIT_TEST ( testEVDcomplex  )

◆ CPPUNIT_TEST() [5/6]

DenseMatrixTest::CPPUNIT_TEST ( testComplexSVD  )

◆ CPPUNIT_TEST() [6/6]

DenseMatrixTest::CPPUNIT_TEST ( testSubMatrix  )

◆ CPPUNIT_TEST_SUITE_END()

DenseMatrixTest::CPPUNIT_TEST_SUITE_END ( )

◆ LIBMESH_CPPUNIT_TEST_SUITE()

DenseMatrixTest::LIBMESH_CPPUNIT_TEST_SUITE ( DenseMatrixTest  )

◆ setUp()

void DenseMatrixTest::setUp ( )
inline

Definition at line 17 of file dense_matrix_test.C.

17 {}

◆ tearDown()

void DenseMatrixTest::tearDown ( )
inline

Definition at line 19 of file dense_matrix_test.C.

19 {}

◆ testComplexSVD()

void DenseMatrixTest::testComplexSVD ( )
inlineprivate

Definition at line 339 of file dense_matrix_test.C.

References libMesh::DenseVector< T >::size(), and libMesh::DenseMatrix< T >::svd().

340  {
341 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
342  LOG_UNIT_TEST;
343 
344  DenseMatrix<Complex> A(3,3);
345 
346  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);
347  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);
348  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);
349 
350  DenseVector<Real> sigma;
351  A.svd(sigma);
352 
353  DenseVector<Real> true_sigma(3);
354  true_sigma(0) = 7.54942516052;
355  true_sigma(1) = 3.17479511368e-06;
356  true_sigma(2) = 6.64680908281e-07;
357 
358  for (unsigned i=0; i<sigma.size(); ++i)
359  LIBMESH_ASSERT_FP_EQUAL(sigma(i), true_sigma(i), 1.e-10);
360 #endif
361  }
std::complex< Real > Complex
virtual unsigned int size() const override final
Definition: dense_vector.h:104
Defines a dense vector for use in Finite Element-type computations.
Defines a dense matrix for use in Finite Element-type computations.
Definition: dof_map.h:65

◆ testEVD_helper()

void DenseMatrixTest::testEVD_helper ( DenseMatrix< Real > &  A,
std::vector< Real >  true_lambda_real,
std::vector< Real >  true_lambda_imag,
Real  tol 
)
inlineprivate

Definition at line 99 of file dense_matrix_test.C.

References std::abs(), libMesh::DenseMatrix< T >::evd_left_and_right(), libMesh::DenseVector< T >::get_values(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseVector< T >::size(), std::sqrt(), and libMesh::DenseVector< T >::zero().

103  {
104  // Note: see bottom of this file, we only do this test if PETSc is
105  // available, but this function currently only exists if we're
106  // using real numbers.
107 #ifdef LIBMESH_USE_REAL_NUMBERS
108  // Let's compute the eigenvalues on a copy of A, so that we can
109  // use the original to check the computation.
110  DenseMatrix<Real> A_copy = A;
111 
112  DenseVector<Real> lambda_real, lambda_imag;
113  DenseMatrix<Real> VR; // right eigenvectors
114  DenseMatrix<Real> VL; // left eigenvectors
115  A_copy.evd_left_and_right(lambda_real, lambda_imag, VL, VR);
116 
117  // The matrix is square and of size N x N.
118  const unsigned N = A.m();
119 
120  // Verify left eigen-values.
121  // Test that the right eigenvalues are self-consistent by computing
122  // u_j**H * A = lambda_j * u_j**H
123  // Note that we have to handle real and complex eigenvalues
124  // differently, since complex eigenvectors share their storage.
125  for (unsigned eigenval=0; eigenval<N; ++eigenval)
126  {
127  // Only check real eigenvalues
128  if (std::abs(lambda_imag(eigenval)) < tol*tol)
129  {
130  // remove print libMesh::out << "Checking eigenvalue: " << eigenval << std::endl;
131  DenseVector<Real> lhs(N), rhs(N);
132  for (unsigned i=0; i<N; ++i)
133  {
134  rhs(i) = lambda_real(eigenval) * VL(i, eigenval);
135  for (unsigned j=0; j<N; ++j)
136  lhs(i) += A(j, i) * VL(j, eigenval); // Note: A(j,i)
137  }
138 
139  // Subtract and assert that the norm of the difference is
140  // below some tolerance.
141  lhs -= rhs;
142  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(tol)*tol);
143  }
144  else
145  {
146  // This is a complex eigenvalue, so:
147  // a.) It occurs in a complex-conjugate pair
148  // b.) the real part of the eigenvector is stored is VL(:,eigenval)
149  // c.) the imag part of the eigenvector is stored in VL(:,eigenval+1)
150  //
151  // Equating the real and imaginary parts of Ax=lambda*x leads to two sets
152  // of relations that must hold:
153  // 1.) A^T x_r = lambda_r*x_r + lambda_i*x_i
154  // 2.) A^T x_i = -lambda_i*x_r + lambda_r*x_i
155  // which we can verify.
156 
157  // 1.)
158  DenseVector<Real> lhs(N), rhs(N);
159  for (unsigned i=0; i<N; ++i)
160  {
161  rhs(i) = lambda_real(eigenval) * VL(i, eigenval) + lambda_imag(eigenval) * VL(i, eigenval+1);
162  for (unsigned j=0; j<N; ++j)
163  lhs(i) += A(j, i) * VL(j, eigenval); // Note: A(j,i)
164  }
165 
166  lhs -= rhs;
167  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(tol)*tol);
168 
169  // libMesh::out << "lhs=" << std::endl;
170  // lhs.print_scientific(libMesh::out, /*precision=*/15);
171  //
172  // libMesh::out << "rhs=" << std::endl;
173  // rhs.print_scientific(libMesh::out, /*precision=*/15);
174 
175  // 2.)
176  lhs.zero();
177  rhs.zero();
178  for (unsigned i=0; i<N; ++i)
179  {
180  rhs(i) = -lambda_imag(eigenval) * VL(i, eigenval) + lambda_real(eigenval) * VL(i, eigenval+1);
181  for (unsigned j=0; j<N; ++j)
182  lhs(i) += A(j, i) * VL(j, eigenval+1); // Note: A(j,i)
183  }
184 
185  lhs -= rhs;
186  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(tol)*tol);
187 
188  // libMesh::out << "lhs=" << std::endl;
189  // lhs.print_scientific(libMesh::out, /*precision=*/15);
190  //
191  // libMesh::out << "rhs=" << std::endl;
192  // rhs.print_scientific(libMesh::out, /*precision=*/15);
193 
194  // We'll skip the second member of the complex conjugate
195  // pair. If the first one worked, the second one should
196  // as well...
197  eigenval += 1;
198  }
199  }
200 
201  // Verify right eigen-values.
202  // Test that the right eigenvalues are self-consistent by computing
203  // A * v_j - lambda_j * v_j
204  // Note that we have to handle real and complex eigenvalues
205  // differently, since complex eigenvectors share their storage.
206  for (unsigned eigenval=0; eigenval<N; ++eigenval)
207  {
208  // Only check real eigenvalues
209  if (std::abs(lambda_imag(eigenval)) < tol*tol)
210  {
211  // remove print libMesh::out << "Checking eigenvalue: " << eigenval << std::endl;
212  DenseVector<Real> lhs(N), rhs(N);
213  for (unsigned i=0; i<N; ++i)
214  {
215  rhs(i) = lambda_real(eigenval) * VR(i, eigenval);
216  for (unsigned j=0; j<N; ++j)
217  lhs(i) += A(i, j) * VR(j, eigenval);
218  }
219 
220  lhs -= rhs;
221  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(tol)*tol);
222  }
223  else
224  {
225  // This is a complex eigenvalue, so:
226  // a.) It occurs in a complex-conjugate pair
227  // b.) the real part of the eigenvector is stored is VR(:,eigenval)
228  // c.) the imag part of the eigenvector is stored in VR(:,eigenval+1)
229  //
230  // Equating the real and imaginary parts of Ax=lambda*x leads to two sets
231  // of relations that must hold:
232  // 1.) Ax_r = lambda_r*x_r - lambda_i*x_i
233  // 2.) Ax_i = lambda_i*x_r + lambda_r*x_i
234  // which we can verify.
235 
236  // 1.)
237  DenseVector<Real> lhs(N), rhs(N);
238  for (unsigned i=0; i<N; ++i)
239  {
240  rhs(i) = lambda_real(eigenval) * VR(i, eigenval) - lambda_imag(eigenval) * VR(i, eigenval+1);
241  for (unsigned j=0; j<N; ++j)
242  lhs(i) += A(i, j) * VR(j, eigenval);
243  }
244 
245  lhs -= rhs;
246  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(tol)*tol);
247 
248  // 2.)
249  lhs.zero();
250  rhs.zero();
251  for (unsigned i=0; i<N; ++i)
252  {
253  rhs(i) = lambda_imag(eigenval) * VR(i, eigenval) + lambda_real(eigenval) * VR(i, eigenval+1);
254  for (unsigned j=0; j<N; ++j)
255  lhs(i) += A(i, j) * VR(j, eigenval+1);
256  }
257 
258  lhs -= rhs;
259  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0., /*actual=*/lhs.l2_norm(), std::sqrt(tol)*tol);
260 
261  // We'll skip the second member of the complex conjugate
262  // pair. If the first one worked, the second one should
263  // as well...
264  eigenval += 1;
265  }
266  }
267 
268  // Sort the results from Lapack *individually*.
269  std::sort(lambda_real.get_values().begin(), lambda_real.get_values().end());
270  std::sort(lambda_imag.get_values().begin(), lambda_imag.get_values().end());
271 
272  // Sort the true eigenvalues *individually*.
273  std::sort(true_lambda_real.begin(), true_lambda_real.end());
274  std::sort(true_lambda_imag.begin(), true_lambda_imag.end());
275 
276  // Compare the individually-sorted values.
277  for (unsigned i=0; i<lambda_real.size(); ++i)
278  {
279  // Note: I initially verified the results with TOLERANCE**2,
280  // but that turned out to be just a bit too tight for some of
281  // the test problems. I'm not sure what controls the accuracy
282  // of the eigenvalue computation in LAPACK, there is no way to
283  // set a tolerance in the LAPACKgeev_ interface.
284  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/true_lambda_real[i], /*actual=*/lambda_real(i), std::sqrt(tol)*tol);
285  LIBMESH_ASSERT_FP_EQUAL(/*expected=*/true_lambda_imag[i], /*actual=*/lambda_imag(i), std::sqrt(tol)*tol);
286  }
287 #endif
288  }
std::vector< T > & get_values()
Definition: dense_vector.h:268
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...
unsigned int m() const
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
virtual unsigned int size() const override final
Definition: dense_vector.h:104
Defines a dense vector for use in Finite Element-type computations.

◆ testEVDcomplex()

void DenseMatrixTest::testEVDcomplex ( )
inlineprivate

Definition at line 311 of file dense_matrix_test.C.

References libMesh::Real, and libMesh::TOLERANCE.

312  {
313  LOG_UNIT_TEST;
314 
315  // This test is also from a Matlab example, and has complex eigenvalues.
316  // http://www.mathworks.com/help/matlab/math/eigenvalues.html?s_tid=gn_loc_drop
317  DenseMatrix<Real> A(3, 3);
318  A(0,0) = 0; A(0,1) = -6; A(0,2) = -1;
319  A(1,0) = 6; A(1,1) = 2; A(1,2) = -16;
320  A(2,0) = -5; A(2,1) = 20; A(2,2) = -10;
321 
322  std::vector<Real> true_lambda_real(3);
323  true_lambda_real[0] = -3.070950351248293;
324  true_lambda_real[1] = -2.464524824375853;
325  true_lambda_real[2] = -2.464524824375853;
326  std::vector<Real> true_lambda_imag(3);
327  true_lambda_imag[0] = 0.;
328  true_lambda_imag[1] = 17.60083096447099;
329  true_lambda_imag[2] = -17.60083096447099;
330 
331  // We're good up to double precision (due to the truncated
332  // literals above) or Real precision, whichever is worse
333  const Real tol = std::max(Real(1e-6), TOLERANCE);
334 
335  // call helper function to compute and verify results
336  testEVD_helper(A, true_lambda_real, true_lambda_imag, tol);
337  }
static constexpr Real TOLERANCE
void testEVD_helper(DenseMatrix< Real > &A, std::vector< Real > true_lambda_real, std::vector< Real > true_lambda_imag, Real tol)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ testEVDreal()

void DenseMatrixTest::testEVDreal ( )
inlineprivate

Definition at line 290 of file dense_matrix_test.C.

References libMesh::TOLERANCE.

291  {
292  LOG_UNIT_TEST;
293 
294  // This is an example from Matlab's gallery(3) which is a
295  // non-symmetric 3x3 matrix with eigen values lambda = 1, 2, 3.
296  DenseMatrix<Real> A(3, 3);
297  A(0,0) = -149; A(0,1) = -50; A(0,2) = -154;
298  A(1,0) = 537; A(1,1) = 180; A(1,2) = 546;
299  A(2,0) = -27; A(2,1) = -9; A(2,2) = -25;
300 
301  std::vector<Real> true_lambda_real(3);
302  true_lambda_real[0] = 1.;
303  true_lambda_real[1] = 2.;
304  true_lambda_real[2] = 3.;
305  std::vector<Real> true_lambda_imag(3); // all zero
306 
307  // call helper function to compute and verify results.
308  testEVD_helper(A, true_lambda_real, true_lambda_imag, TOLERANCE);
309  }
static constexpr Real TOLERANCE
void testEVD_helper(DenseMatrix< Real > &A, std::vector< Real > true_lambda_real, std::vector< Real > true_lambda_imag, Real tol)

◆ testOuterProduct()

void DenseMatrixTest::testOuterProduct ( )
inlineprivate

Definition at line 35 of file dense_matrix_test.C.

References libMesh::libmesh_real(), libMesh::DenseMatrix< T >::outer_product(), libMesh::DenseVector< T >::size(), and libMesh::TOLERANCE.

36  {
37  LOG_UNIT_TEST;
38 
39  DenseVector<Real> a = {1.0, 2.0};
40  DenseVector<Real> b = {3.0, 4.0, 5.0};
41 
42  DenseMatrix<Real> a_times_b;
43  a_times_b.outer_product(a, b);
44 
45  DenseMatrix<Real> a_times_b_correct(2, 3,
46  {3., 4., 5.,
47  6., 8., 10.});
48 
49  for (unsigned int i = 0; i < a.size(); ++i)
50  for (unsigned int j = 0; j < b.size(); ++j)
51  LIBMESH_ASSERT_FP_EQUAL(libmesh_real(a_times_b(i,j)), libmesh_real(a_times_b_correct(i,j)), TOLERANCE*TOLERANCE);
52  }
T libmesh_real(T a)
static constexpr Real TOLERANCE
void outer_product(const DenseVector< T > &a, const DenseVector< T > &b)
Computes the outer (dyadic) product of two vectors and stores in (*this).
virtual unsigned int size() const override final
Definition: dense_vector.h:104
Defines a dense vector for use in Finite Element-type computations.

◆ testSubMatrix()

void DenseMatrixTest::testSubMatrix ( )
inlineprivate

Definition at line 363 of file dense_matrix_test.C.

References libMesh::DenseMatrix< T >::sub_matrix().

364  {
365  LOG_UNIT_TEST;
366 
367  DenseMatrix<Number> A(4, 3);
368  A(0,0) = 1.0; A(0,1) = 2.0; A(0,2) = 3.0;
369  A(1,0) = 4.0; A(1,1) = 5.0; A(1,2) = 6.0;
370  A(2,0) = 7.0; A(2,1) = 8.0; A(2,2) = 9.0;
371  A(3,0) =10.0; A(3,1) =11.0; A(3,2) =12.0;
372 
373  DenseMatrix<Number> B(2, 2);
374  B(0,0) = 7.0; B(0,1) = 8.0;
375  B(1,0) =10.0; B(1,1) =11.0;
376 
377  DenseMatrix<Number> C = A.sub_matrix(2, 2, 0, 2);
378  CPPUNIT_ASSERT(B == C);
379  }
Definition: assembly.h:38
DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
Get submatrix with the smallest row and column indices and the submatrix size.
Definition: dense_matrix.h:922

◆ testSVD()

void DenseMatrixTest::testSVD ( )
inlineprivate

Definition at line 54 of file dense_matrix_test.C.

References libMesh::libmesh_real(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::Real, libMesh::DenseVector< T >::size(), libMesh::DenseMatrix< T >::svd(), and libMesh::TOLERANCE.

55  {
56  LOG_UNIT_TEST;
57 
58  DenseMatrix<Number> U, VT;
59  DenseVector<Real> sigma;
60  DenseMatrix<Number> A(3, 2, {
61  1.0, 2.0,
62  3.0, 4.0,
63  5.0, 6.0});
64 
65  A.svd(sigma, U, VT);
66 
67  // Solution for this case is (verified with numpy)
68  DenseMatrix<Number> true_U(3, 2, {
69  -2.298476964000715e-01, 8.834610176985250e-01,
70  -5.247448187602936e-01, 2.407824921325463e-01,
71  -8.196419411205157e-01, -4.018960334334318e-01});
72 
73  DenseMatrix<Number> true_VT(2, 2, {
74  -6.196294838293402e-01, -7.848944532670524e-01,
75  -7.848944532670524e-01, 6.196294838293400e-01});
76 
77  DenseVector<Real> true_sigma = {9.525518091565109e+00, 5.143005806586446e-01};
78 
79  // Tolerance is bounded by the double literals above and by using Real
80  const Real tol = std::max(Real(1e-12), TOLERANCE*TOLERANCE);
81 
82  for (unsigned i=0; i<U.m(); ++i)
83  for (unsigned j=0; j<U.n(); ++j)
84  LIBMESH_ASSERT_FP_EQUAL( libmesh_real(U(i,j)), libmesh_real(true_U(i,j)), tol);
85 
86  for (unsigned i=0; i<VT.m(); ++i)
87  for (unsigned j=0; j<VT.n(); ++j)
88  LIBMESH_ASSERT_FP_EQUAL( libmesh_real(VT(i,j)), libmesh_real(true_VT(i,j)), tol);
89 
90  for (unsigned i=0; i<sigma.size(); ++i)
91  LIBMESH_ASSERT_FP_EQUAL(sigma(i), true_sigma(i), tol);
92  }
T libmesh_real(T a)
static constexpr Real TOLERANCE
unsigned int m() const
void svd(DenseVector< Real > &sigma)
Compute the singular value decomposition of the matrix.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned int size() const override final
Definition: dense_vector.h:104
Defines a dense vector for use in Finite Element-type computations.
unsigned int n() const

The documentation for this class was generated from the following file: