// $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
};
-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);
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();
<< (int) (lambda.real()+.0001) << '\t'
<< (int) (lambda.imag()+.0001) << std::endl;
}
-
- v1.reinit(4);
- v2.reinit(4);
}
-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
// $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
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);
--- /dev/null
+//--------------------------------------------------------------------
+// $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;
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+//--------------------------------------------------------------------
+// $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);
+}
--- /dev/null
+
+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
--- /dev/null
+
+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
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
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
// $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
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);