check for FullMatrix::invert

git-svn-id: https://svn.dealii.org/trunk@3275 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/tests/lac/full_matrix.cc b/tests/lac/full_matrix.cc
index 3aa2f3f1fc1a7fc686117557160879b2006880ea..ab77c1752c70c68408fdc79662b5708dfc7cc590 100644 (file)
--- a/tests/lac/full_matrix.cc
+++ b/tests/lac/full_matrix.cc
@@ -15,12 +15,26 @@
 #include <cmath>
 #include <fstream>
 #include <iomanip>
+#include <cstdlib>
 
 #include <base/logstream.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/eigen.h>
 
+// Create a positive definite random matrix
+
+void random_matrix(FullMatrix<double>& A)
+{
+  for (unsigned int i=0; i<A.m();++i)
+    for (unsigned int j=0; j<A.n();++j)
+      {
+       double rnd = rand();
+       rnd /= RAND_MAX;
+       A(i,j) = (i==j) ? A.m()+rnd : rnd;
+      }
+}
+
 int
 main ()
 {
@@ -29,84 +43,113 @@ main ()
   logfile.precision(3);
   deallog.attach(logfile);
   deallog.depth_console(0);
+  srand(3391466);
   
-  FullMatrix<double> A(5,5), C(5,5), D(5,5), H(5,5);
-  D(0,0) = 1.;
-  D(1,1) = 2.;
-  D(2,2) = 3.;
-  D(3,3) = 4.;
-  D(4,4) = 5.;
-
-  A = D;
-  
-  for (unsigned int i=0;i<4;++i)
+  for (unsigned int i=1;i<10;++i)
     {
-                                      // Setup rotation matrix
-      C.clear();
-      C.diagadd(1.);
-      C(i,i) = C(i+1,i+1) = cos(i+1);
-      C(i+1,i) = sin(i+1);
-      C(i,i+1) = -sin(i+1);
+      FullMatrix<double> A(i,i), B(i,i);
+ 
+                                      // Create matrix and its inverse
+      random_matrix(A);
+      B.invert(A);
 
-      C.print_formatted (logfile);
-      deallog << "l1-norm: " << C.l1_norm() << endl;
-      D = C;
-      D.gauss_jordan();
-      D.print_formatted (logfile);
-      deallog << "linfty-norm: " << D.linfty_norm() << endl
-             << "Frobenius-norm: " << D.norm2() << endl;
-
-                                      // Rotate original matrix
-      A.mmult(H,C);
-      C.Tmmult(A,H);
+                                      // Check if unit vectors are recovered
+      deallog << "Inverse(dim=" << i <<"):";
+      for (unsigned int j=0;j<i;++j)
+       {
+         Vector<double> x(i);
+         Vector<double> y(i);
+         Vector<double> z(i);
+         x(j) = 1.;
+         A.vmult(y,x);
+         B.vmult(z,y);
+         z.add(-1.,x);
+         double a = z.l2_norm();
+         if (a > 1.e-12) deallog << a << ' ';
+       }
+      deallog << endl;
     }
 
-  A.print_formatted (logfile);
-
-  Vector<double> u(5);
-  GrowingVectorMemory<Vector<double> > mem;
-  
-  SolverControl control (500,1.e-8, false, true);
-  
-  if (true)
-    {
-      u = 1.;
-      EigenPower<Vector<double> >
-       von_Mises(control, mem, 0.);
-      double eigen = 0.;
-      von_Mises.solve(eigen, A, u);
-      deallog << "Eigenvalue: " << eigen << endl;
-    }
-  if (true)
-    {
-      u = 1.;
-      EigenPower<Vector<double> >
-       von_Mises(control, mem, -4.);
-      double eigen = 0.;
-      von_Mises.solve(eigen, A, u);
-      deallog << "Eigenvalue: " << eigen << endl;
-    }
-  H = A;
-  H.gauss_jordan();
-  H.print_formatted (logfile);
   if (true)
     {
-      u = 1.;
-      EigenPower<Vector<double> >
-       von_Mises(control, mem, 0.);
-      double eigen = 0.;
-      von_Mises.solve(eigen, H, u);
-      deallog << "Eigenvalue: " << eigen << endl;
+      FullMatrix<double> A(5,5), C(5,5), D(5,5), H(5,5);
+      D(0,0) = 1.;
+      D(1,1) = 2.;
+      D(2,2) = 3.;
+      D(3,3) = 4.;
+      D(4,4) = 5.;
+
+      A = D;
+  
+      for (unsigned int i=0;i<4;++i)
+       {
+                                          // Setup rotation matrix
+         C.clear();
+         C.diagadd(1.);
+         C(i,i) = C(i+1,i+1) = cos(i+1);
+         C(i+1,i) = sin(i+1);
+         C(i,i+1) = -sin(i+1);
+         
+         C.print_formatted (logfile);
+         deallog << "l1-norm: " << C.l1_norm() << endl;
+         D = C;
+         D.gauss_jordan();
+         D.print_formatted (logfile);
+         deallog << "linfty-norm: " << D.linfty_norm() << endl
+                 << "Frobenius-norm: " << D.norm2() << endl;
+         
+                                          // Rotate original matrix
+         A.mmult(H,C);
+         C.Tmmult(A,H);
+       }
+      
+      A.print_formatted (logfile);
+      
+      Vector<double> u(5);
+      GrowingVectorMemory<Vector<double> > mem;
+      
+      SolverControl control (500,1.e-8, false, true);
+      
+      if (true)
+       {
+         u = 1.;
+         EigenPower<Vector<double> >
+           von_Mises(control, mem, 0.);
+         double eigen = 0.;
+         von_Mises.solve(eigen, A, u);
+         deallog << "Eigenvalue: " << eigen << endl;
+       }
+      if (true)
+       {
+         u = 1.;
+         EigenPower<Vector<double> >
+           von_Mises(control, mem, -4.);
+         double eigen = 0.;
+         von_Mises.solve(eigen, A, u);
+         deallog << "Eigenvalue: " << eigen << endl;
+       }
+      H = A;
+      H.gauss_jordan();
+      H.print_formatted (logfile);
+      if (true)
+       {
+         u = 1.;
+         EigenPower<Vector<double> >
+           von_Mises(control, mem, 0.);
+         double eigen = 0.;
+         von_Mises.solve(eigen, H, u);
+         deallog << "Eigenvalue: " << eigen << endl;
+       }
+      if (true)
+       {
+         u = 1.;
+         EigenPower<Vector<double> >
+           von_Mises(control, mem, -4.);
+         double eigen = 0.;
+         von_Mises.solve(eigen, H, u);
+         deallog << "Eigenvalue: " << eigen << endl;
+       }
     }
-  if (true)
-    {
-      u = 1.;
-      EigenPower<Vector<double> >
-       von_Mises(control, mem, -4.);
-      double eigen = 0.;
-      von_Mises.solve(eigen, H, u);
-      deallog << "Eigenvalue: " << eigen << endl;
-    }   
 }