]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
some restructuring and test for SVD
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 7 Nov 2010 23:11:47 +0000 (23:11 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 7 Nov 2010 23:11:47 +0000 (23:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@22629 0785d39b-7218-0410-832d-ea1e28bc413d

tests/lapack/full_matrix_01.cc
tests/lapack/full_matrix_01/cmp/generic
tests/lapack/full_matrix_02.cc
tests/lapack/full_matrix_03.cc [new file with mode: 0644]
tests/lapack/full_matrix_03/cmp/generic [new file with mode: 0644]
tests/lapack/full_matrix_04.cc [new file with mode: 0644]
tests/lapack/full_matrix_04/cmp/generic [new file with mode: 0644]
tests/lapack/full_matrix_04/full_matrix_03/cmp/generic [new file with mode: 0644]
tests/lapack/solver_cg/cmp/generic
tests/lapack/tridiagonal_ev1.cc

index 8798bfd59944e1eeae53a32709cd69e3d53bca77..56467c6559fd801d78e04e3890a5b0b16ac3f227 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$ 
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -45,30 +45,23 @@ const double rect[] =
 };
 
 
-int main()
+void test_rect(unsigned int m, unsigned int n, const double* values)
 {
-  std::ofstream logfile("full_matrix_01/output");
-  logfile.precision(3);
-  deallog.attach(logfile);
-  deallog.depth_console(0);
-  deallog.threshold_double(1.e-10);
-
-  FullMatrix<double> A(3,4,rect);
-  LAPACKFullMatrix<double> LA(3,4);
+  std::ostringstream prefix;
+  prefix << m << 'x' << n;
+  deallog.push(prefix.str());
+  
+  FullMatrix<double> A(m,n,values);
+  LAPACKFullMatrix<double> LA(m,n);
   LA = A;
   
-  Vector<double> u(4);
-  Vector<double> v1(3);
-  Vector<double> v2(3);
+  Vector<double> u(n);
+  Vector<double> v1(m);
+  Vector<double> v2(m);
   
   for (unsigned int i=0;i<u.size();++i)
     u(i) = i*i;
   
-                                  // Test rectangular vmult. All
-                                  // results compare with same
-                                  // operation for FullMatrix.
-  deallog.push("Rect");
-
   deallog << "operator= (const FullMatrix<number>&) ok" << std::endl;
   
   A.vmult(v1,u);
@@ -97,10 +90,24 @@ int main()
     deallog << "Tvmult_add ok" << std::endl;
   
   deallog.pop();
+}
+
+int main()
+{
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test_rect(3,4,rect);
+  test_rect(4,3,rect);
+  test_rect(4,4,symm);
 
                                   // Test symmetric system
-  A.reinit(4,4);
-  LA.reinit(4,4);
+  FullMatrix<double> A(4,4,symm);
+  LAPACKFullMatrix<double> LA(4,4);
   A.fill(symm);
   LA = A;
   LA.compute_eigenvalues();
@@ -111,7 +118,4 @@ int main()
              << (int) (lambda.real()+.0001) << '\t'
              << (int) (lambda.imag()+.0001) << std::endl;
     }
-  
-  v1.reinit(4);
-  v2.reinit(4);  
 }
index c8eede55a5f60d61d619855f7682720d1172aff7..90b8cf7560c75ab74eb16498d795361cf4a090dc 100644 (file)
@@ -1,9 +1,19 @@
 
-DEAL:Rect::operator= (const FullMatrix<number>&) ok
-DEAL:Rect::vmult ok
-DEAL:Rect::vmult_add ok
-DEAL:Rect::Tvmult ok
-DEAL:Rect::Tvmult_add ok
+DEAL:3x4::operator= (const FullMatrix<number>&) ok
+DEAL:3x4::vmult ok
+DEAL:3x4::vmult_add ok
+DEAL:3x4::Tvmult ok
+DEAL:3x4::Tvmult_add ok
+DEAL:4x3::operator= (const FullMatrix<number>&) ok
+DEAL:4x3::vmult ok
+DEAL:4x3::vmult_add ok
+DEAL:4x3::Tvmult ok
+DEAL:4x3::Tvmult_add ok
+DEAL:4x4::operator= (const FullMatrix<number>&) ok
+DEAL:4x4::vmult ok
+DEAL:4x4::vmult_add ok
+DEAL:4x4::Tvmult ok
+DEAL:4x4::Tvmult_add ok
 DEAL::Eigenvalues 5    0
 DEAL::Eigenvalues 1    0
 DEAL::Eigenvalues 5    0
index cd450766c1b7b3b61268e6a19d05b19ee06d42dd..7777f271e905d7e0dd543b833f533b98bec169f6 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$ 
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -75,7 +75,8 @@ test(const unsigned int size, const bool nonsymmetric)
 
 int main()
 {
-  std::ofstream logfile("full_matrix_02/output");
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
   logfile.precision(3);
   deallog.attach(logfile);
   deallog.depth_console(0);
diff --git a/tests/lapack/full_matrix_03.cc b/tests/lapack/full_matrix_03.cc
new file mode 100644 (file)
index 0000000..4d6498b
--- /dev/null
@@ -0,0 +1,136 @@
+//--------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006, 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//--------------------------------------------------------------------
+
+// Tests SVD of LAPACKFullMatrix by comparing to vmult of FullMatrix
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <lac/lapack_full_matrix.h>
+#include <lac/full_matrix.h>
+#include <lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+/*
+ * Eigenvalues and -vectors of this system are
+ * lambda = 1     v = (1, 1, 1, 1)
+ * lambda = 5     v = (1,-1, 0, 0)
+ * lambda = 5     v = (0, 1,-1, 0)
+ * lambda = 5     v = (0, 0, 1,-1)
+ */
+const double symm[] =
+{
+      4., -1., -1., -1.,
+      -1., 4., -1., -1.,
+      -1., -1., 4., -1.,
+      -1., -1., -1., 4.
+};
+
+const double rect[] =
+{
+      4., 3., 2., 1.,
+      5., 8., 1., -2.,
+      11., 13., -4., -5
+};
+
+
+void test_rect(unsigned int m, unsigned int n, const double* values)
+{
+  std::ostringstream prefix;
+  prefix << m << 'x' << n;
+  deallog.push(prefix.str());
+  
+  FullMatrix<double> A(m,n,values);
+  LAPACKFullMatrix<double> LA(m,n);
+  LA = A;
+  LA.compute_svd();
+  
+  deallog << "Singular values";
+  for (unsigned int i=0;i<LA.n_rows();++i)
+    deallog << ' ' << LA.singular_value(i);
+  deallog << std::endl;
+  
+  Vector<double> u(n);
+  Vector<double> v1(m);
+  Vector<double> v2(m);
+  
+  for (unsigned int i=0;i<u.size();++i)
+    u(i) = i*i;
+  
+                                  // Test rectangular vmult. All
+                                  // results compare with same
+                                  // operation for FullMatrix.
+  
+  A.vmult(v1,u);
+  LA.vmult(v2,u);
+  v1 -= v2;
+  if (v1.l2_norm() < 1.e-12)
+    deallog << "vmult ok" << std::endl;
+  else
+    deallog << "vmult error " << v1.l2_norm() << std::endl;
+  v1 = v2;
+  
+  A.vmult_add(v1,u);
+  LA.vmult_add(v2,u);
+  v1 -= v2;
+  if (v1.l2_norm() < 1.e-12)
+    deallog << "vmult_add ok" << std::endl;
+  else
+    deallog << "vmult_add error " << v1.l2_norm() << std::endl;
+  
+  LA.Tvmult(u, v2);
+  u *= -1;
+  A.Tvmult_add(u, v2);
+  if (u.l2_norm() < 1.e-12)
+    deallog << "Tvmult ok" << std::endl;
+  else
+    deallog << "Tvmult error " << u.l2_norm() << std::endl;
+  
+  A.Tvmult(u, v2);
+  u *= -1;
+  LA.Tvmult_add(u, v2);
+  if (u.l2_norm() < 1.e-12)
+    deallog << "Tvmult_add ok" << std::endl;
+  else
+    deallog << "Tvmult_add error " << u.l2_norm() << std::endl;
+  
+  deallog.pop();
+}
+
+int main()
+{
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test_rect(4,4,symm);
+  test_rect(4,3,rect);
+  test_rect(3,4,rect);
+  
+                                  // Test symmetric system
+  FullMatrix<double> A(4,4,symm);
+  LAPACKFullMatrix<double> LA(4,4);
+  A.fill(symm);
+  LA = A;
+  LA.compute_eigenvalues();
+  for (unsigned int i=0;i<A.m();++i)
+    {
+      std::complex<double> lambda = LA.eigenvalue(i);
+      deallog << "Eigenvalues "
+             << (int) (lambda.real()+.0001) << '\t'
+             << (int) (lambda.imag()+.0001) << std::endl;
+    }
+}
diff --git a/tests/lapack/full_matrix_03/cmp/generic b/tests/lapack/full_matrix_03/cmp/generic
new file mode 100644 (file)
index 0000000..57ea1fd
--- /dev/null
@@ -0,0 +1,20 @@
+
+DEAL:4x4::Singular values 5.00000 5.00000 5.00000 1.00000
+DEAL:4x4::vmult ok
+DEAL:4x4::vmult_add ok
+DEAL:4x4::Tvmult ok
+DEAL:4x4::Tvmult_add ok
+DEAL:4x3::Singular values 16.0315 12.6623 6.13650 0
+DEAL:4x3::vmult ok
+DEAL:4x3::vmult_add ok
+DEAL:4x3::Tvmult ok
+DEAL:4x3::Tvmult_add ok
+DEAL:3x4::Singular values 20.7967 4.30713 1.98624
+DEAL:3x4::vmult ok
+DEAL:3x4::vmult_add ok
+DEAL:3x4::Tvmult ok
+DEAL:3x4::Tvmult_add ok
+DEAL::Eigenvalues 5    0
+DEAL::Eigenvalues 1    0
+DEAL::Eigenvalues 5    0
+DEAL::Eigenvalues 5    0
diff --git a/tests/lapack/full_matrix_04.cc b/tests/lapack/full_matrix_04.cc
new file mode 100644 (file)
index 0000000..0595e6b
--- /dev/null
@@ -0,0 +1,124 @@
+//--------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2006, 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//--------------------------------------------------------------------
+
+// Tests inverse SVD of LAPACKFullMatrix by comparing to vmult of FullMatrix
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <lac/lapack_full_matrix.h>
+#include <lac/full_matrix.h>
+#include <lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+/*
+ * Eigenvalues and -vectors of this system are
+ * lambda = 1     v = (1, 1, 1, 1)
+ * lambda = 5     v = (1,-1, 0, 0)
+ * lambda = 5     v = (0, 1,-1, 0)
+ * lambda = 5     v = (0, 0, 1,-1)
+ */
+const double symm[] =
+{
+      4., -1., -1., -1.,
+      -1., 4., -1., -1.,
+      -1., -1., 4., -1.,
+      -1., -1., -1., 4.
+};
+
+const double rect[] =
+{
+      4., 3., 2., 1.,
+      5., 8., 1., -2.,
+      11., 13., -4., -5
+};
+
+
+void test_rect(unsigned int m, unsigned int n, const double* values)
+{
+  std::ostringstream prefix;
+  prefix << m << 'x' << n;
+  deallog.push(prefix.str());
+  
+  FullMatrix<double> A(m,n,values);
+  LAPACKFullMatrix<double> LA(m,n);
+  LA = A;
+  LA.compute_inverse_svd();
+  
+  deallog << "Singular values";
+  for (unsigned int i=0;i<LA.n_rows();++i)
+    deallog << ' ' << LA.singular_value(i);
+  deallog << std::endl;
+  
+  Vector<double> u1(n);
+  Vector<double> u2(n);
+  Vector<double> v1(m);
+  Vector<double> v2(m);
+  
+  for (unsigned int i=0;i<u1.size();++i)
+    u1(i) = i*i;
+  for (unsigned int i=0;i<v1.size();++i)
+    v1(i) = i*i;
+
+                                  // Test if LA is a left inverse of A
+  A.vmult(v2,u1);
+  LA.vmult(u2,v2);
+  u2 -= u1;
+  if (u2.l2_norm() < 1.e-12)
+    deallog << "vmult ok" << std::endl;
+  else
+    deallog << "vmult error " << u2.l2_norm() << std::endl;
+  
+                                  // Test if LA is a right inverse of A
+  LA.vmult(u2,v1);
+  A.vmult(v2,u2);
+  v2 -= v1;
+  if (v2.l2_norm() < 1.e-12)
+    deallog << "vmult ok" << std::endl;
+  else
+    deallog << "vmult error " << v2.l2_norm() << std::endl;  
+
+                                  // Test if LA^T is a left inverse
+                                  // of A^T
+  A.Tvmult(u2,v1);
+  LA.Tvmult(v2, u2);
+  v2 -= v1;
+  if (v2.l2_norm() < 1.e-12)
+    deallog << "Tvmult ok" << std::endl;
+  else
+    deallog << "Tvmult error " << v2.l2_norm() << std::endl;
+  
+  LA.Tvmult(v2,u1);
+  A.Tvmult(u2,v2);
+  u2 -= u1;
+  if (u2.l2_norm() < 1.e-12)
+    deallog << "Tvmult ok" << std::endl;
+  else
+    deallog << "Tvmult error " << u2.l2_norm() << std::endl;
+  
+  deallog.pop();
+}
+
+int main()
+{
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test_rect(4,4,symm);
+  test_rect(4,3,rect);
+  test_rect(3,4,rect);
+}
diff --git a/tests/lapack/full_matrix_04/cmp/generic b/tests/lapack/full_matrix_04/cmp/generic
new file mode 100644 (file)
index 0000000..68b9a51
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL:4x4::Singular values 0.200000 0.200000 0.200000 1.00000
+DEAL:4x4::vmult ok
+DEAL:4x4::vmult ok
+DEAL:4x4::Tvmult ok
+DEAL:4x4::Tvmult ok
+DEAL:4x3::Singular values 0.0623770 0.0789748 0.162959 0
+DEAL:4x3::vmult ok
+DEAL:4x3::vmult error 1.75406
+DEAL:4x3::Tvmult error 1.75406
+DEAL:4x3::Tvmult ok
+DEAL:3x4::Singular values 0.0480845 0.232173 0.503464
+DEAL:3x4::vmult error 5.58692
+DEAL:3x4::vmult ok
+DEAL:3x4::Tvmult ok
+DEAL:3x4::Tvmult error 5.58692
diff --git a/tests/lapack/full_matrix_04/full_matrix_03/cmp/generic b/tests/lapack/full_matrix_04/full_matrix_03/cmp/generic
new file mode 100644 (file)
index 0000000..57ea1fd
--- /dev/null
@@ -0,0 +1,20 @@
+
+DEAL:4x4::Singular values 5.00000 5.00000 5.00000 1.00000
+DEAL:4x4::vmult ok
+DEAL:4x4::vmult_add ok
+DEAL:4x4::Tvmult ok
+DEAL:4x4::Tvmult_add ok
+DEAL:4x3::Singular values 16.0315 12.6623 6.13650 0
+DEAL:4x3::vmult ok
+DEAL:4x3::vmult_add ok
+DEAL:4x3::Tvmult ok
+DEAL:4x3::Tvmult_add ok
+DEAL:3x4::Singular values 20.7967 4.30713 1.98624
+DEAL:3x4::vmult ok
+DEAL:3x4::vmult_add ok
+DEAL:3x4::Tvmult ok
+DEAL:3x4::Tvmult_add ok
+DEAL::Eigenvalues 5    0
+DEAL::Eigenvalues 1    0
+DEAL::Eigenvalues 5    0
+DEAL::Eigenvalues 5    0
index 343ef8235c32fbcf3764b6b586832aca1ab0ddef..29c895336109bb4a697083887e6e11b4d8544c74 100644 (file)
@@ -16,7 +16,7 @@ DEAL:ssor:cg::Starting value 3.000
 DEAL:ssor:cg::Condition number estimate: 1.384
 DEAL:ssor:cg::Condition number estimate: 1.421
 DEAL:ssor:cg::Convergence step 4 value 6.018e-05
-DEAL:ssor:cg:: 0.5858 0.6810 0.8324
+DEAL:ssor:cg:: 0.7030 0.8172 0.9988
 DEAL::Size 12 Unknowns 121
 DEAL:no-fail:cg::Starting value 11.00
 DEAL:no-fail:cg::Condition number estimate: 11.01
@@ -70,4 +70,4 @@ DEAL:ssor:cg::Condition number estimate: 5.706
 DEAL:ssor:cg::Condition number estimate: 5.715
 DEAL:ssor:cg::Condition number estimate: 5.718
 DEAL:ssor:cg::Convergence step 8 value 0.0005816
-DEAL:ssor:cg:: 0.1457 0.3791 0.4599 0.5650 0.6652 0.7548 0.8330
+DEAL:ssor:cg:: 0.1748 0.4549 0.5519 0.6780 0.7983 0.9058 0.9996
index 30663558047e5b5e7ec5019157cd71ed4b10a4a7..e231b301936263a11cb8f60230080e6ad7c8cfbb 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$ 
 //
-//    Copyright (C) 2006 by the deal.II authors
+//    Copyright (C) 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -42,8 +42,9 @@ void test_laplacian(unsigned int n)
 
 
 int main()
-{
-  std::ofstream logfile("tridiagonal_ev1/output");
+{ 
+  const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+  std::ofstream logfile(logname.c_str());
   deallog.attach(logfile);
   deallog.depth_console(0);
   deallog.threshold_double(1.e-10);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.