--- /dev/null
+//----------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2007, 2009 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.
+//
+//----------------------------------------------------------------------
+
+
+// check FullMatrix::cholesky for correct functionality
+
+#include "../tests.h"
+#include "full_matrix_common.h"
+
+
+std::string output_file_name = "full_matrix_56/output";
+
+const double entries2[4] = { 0.814723686393179 , 0.516389376684563,
+ 0.516389376684563, 0.913375856139019};
+
+
+const double entries3[9] = { 1.808621732261680, 0.845866689167942,
+ 0.667106522517665, 0.845866689167942,
+ 1.398873167281503, 0.281706853672865,
+ 0.667106522517665, 0.281706853672865,
+ 0.741757766593798};
+
+
+
+// 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 = std::rand();
+ rnd /= RAND_MAX;
+ A(i,j) = (i==j) ? A.m()+rnd : rnd;
+ }
+}
+
+template <typename number>
+void
+check ()
+ {
+ deallog << std::fixed;
+ deallog << std::setprecision(4);
+// deallog.attach(logfile);
+ deallog.depth_console(0);
+// deallog.threshold_double(1.e-10);
+ std::srand(3391466);
+
+ FullMatrix<double> T2(2,2,entries2), S2(2,2), R2(2,2);
+ FullMatrix<double>T3(3,3,entries3), S3(3,3), R3(3,3);
+
+ // These results for T2, T3, T5 were compared
+ // against Matlab.
+ deallog << "Original Matrix:" << std::endl;
+ deallog << "================" << std::endl;
+
+ display_matrix(T2);
+ deallog << std::endl;
+
+ S2.cholesky(T2);
+ deallog << "Cholesky Decomposition Matrix:" << std::endl;
+ deallog << "==============================" << std::endl;
+ display_matrix(S2);
+ deallog << std::endl;
+
+ R2.Tadd(1,S2); // R = S^t
+ deallog << "Transposed Cholesky Decomposition Matrix:" << std::endl;
+ deallog << "=========================================" << std::endl;
+ display_matrix(R2);
+ deallog << std::endl;
+
+ R2.Tmmult(S2, R2, false); // S = R^tR
+ deallog << "Reconstructed Matrix:" << std::endl;
+ deallog << "=====================" << std::endl;
+ display_matrix(S2);
+ deallog << std::endl;
+ S2.add(-1.0,T2);
+ if (S2.frobenius_norm()>1.0e-10)
+ {
+ deallog << "NOT the same to 1e-10 tolerance" << std::endl;
+ }
+ else
+ {
+ deallog << "the SAME to 1e-10 tolerance" << std::endl;
+ }
+ deallog << std::endl;
+ deallog << std::endl;
+
+ deallog << "Original Matrix:" << std::endl;
+ deallog << "================" << std::endl;
+ display_matrix(T3);
+ deallog << std::endl;
+
+ S3.cholesky(T3);
+ deallog << "Cholesky Decomposition Matrix:" << std::endl;
+ deallog << "==============================" << std::endl;
+ display_matrix(S3);
+ deallog << std::endl;
+
+ R3.Tadd(1,S3); // R = S^t
+ deallog << "Transposed Cholesky Decomposition Matrix:" << std::endl;
+ deallog << "=========================================" << std::endl;
+ display_matrix(R3);
+ deallog << std::endl;
+
+ R3.Tmmult(S3, R3, false); // S = R^tR
+ deallog << "Reconstructed Matrix:" << std::endl;
+ deallog << "=====================" << std::endl;
+ display_matrix(S3);
+ deallog << std::endl;
+ S3.add(-1.0,T3);
+ if (S3.frobenius_norm()>1.0e-10)
+ {
+ deallog << "NOT the same to 1e-10 tolerance" << std::endl;
+ }
+ else
+ {
+ deallog << "the SAME to 1e-10 tolerance" << std::endl;
+ }
+ deallog << std::endl;
+ deallog << std::endl;
+
+ FullMatrix<double> T5(5,5), S5(5,5), R5(5,5);
+ random_matrix(T5);
+ T5.symmetrize();
+
+ deallog << "Original Matrix:" << std::endl;
+ deallog << "================" << std::endl;
+ display_matrix(T5);
+ deallog << std::endl;
+
+ S5.cholesky(T5);
+ deallog << "Cholesky Decomposition Matrix:" << std::endl;
+ deallog << "==============================" << std::endl;
+ display_matrix(S5);
+ deallog << std::endl;
+
+ R5.Tadd(1,S5); // R = S^t
+ deallog << "Transposed Cholesky Decomposition Matrix:" << std::endl;
+ deallog << "=========================================" << std::endl;
+ display_matrix(R5);
+ deallog << std::endl;
+
+ R5.Tmmult(S5, R5, false); // S = R^tR
+ deallog << "Reconstructed Matrix:" << std::endl;
+ deallog << "=====================" << std::endl;
+ display_matrix(S5);
+
+ deallog << std::endl;
+ S5.add(-1.0,T5);
+ if (S5.frobenius_norm()>1.0e-10)
+ {
+ deallog << "NOT the same to 1e-10 tolerance" << std::endl;
+ }
+ else
+ {
+ deallog << "the SAME to 1e-10 tolerance" << std::endl;
+ }
+ deallog << std::endl;
+ deallog << std::endl;
+
+}
+
+
+
--- /dev/null
+
+DEAL::Original Matrix:
+DEAL::================
+DEAL::2x2 matrix
+DEAL::0.8147 0.5164
+DEAL::0.5164 0.9134
+DEAL::
+DEAL::Cholesky Decomposition Matrix:
+DEAL::==============================
+DEAL::2x2 matrix
+DEAL::0.9026 0
+DEAL::0.5721 0.7656
+DEAL::
+DEAL::Transposed Cholesky Decomposition Matrix:
+DEAL::=========================================
+DEAL::2x2 matrix
+DEAL::0.9026 0.5721
+DEAL::0 0.7656
+DEAL::
+DEAL::Reconstructed Matrix:
+DEAL::=====================
+DEAL::2x2 matrix
+DEAL::0.8147 0.5164
+DEAL::0.5164 0.9134
+DEAL::
+DEAL::the SAME to 1e-10 tolerance
+DEAL::
+DEAL::
+DEAL::Original Matrix:
+DEAL::================
+DEAL::3x3 matrix
+DEAL::1.8086 0.8459 0.6671
+DEAL::0.8459 1.3989 0.2817
+DEAL::0.6671 0.2817 0.7418
+DEAL::
+DEAL::Cholesky Decomposition Matrix:
+DEAL::==============================
+DEAL::3x3 matrix
+DEAL::1.3449 0 0
+DEAL::0.6290 1.0016 0
+DEAL::0.4960 -0.0302 0.7034
+DEAL::
+DEAL::Transposed Cholesky Decomposition Matrix:
+DEAL::=========================================
+DEAL::3x3 matrix
+DEAL::1.3449 0.6290 0.4960
+DEAL::0 1.0016 -0.0302
+DEAL::0 0 0.7034
+DEAL::
+DEAL::Reconstructed Matrix:
+DEAL::=====================
+DEAL::3x3 matrix
+DEAL::1.8086 0.8459 0.6671
+DEAL::0.8459 1.3989 0.2817
+DEAL::0.6671 0.2817 0.7418
+DEAL::
+DEAL::the SAME to 1e-10 tolerance
+DEAL::
+DEAL::
+DEAL::Original Matrix:
+DEAL::================
+DEAL::5x5 matrix
+DEAL::5.5862 0.5230 0.1500 0.6561 0.5712
+DEAL::0.5230 5.5897 0.5640 0.5585 0.6221
+DEAL::0.1500 0.5640 5.5445 0.3284 0.6072
+DEAL::0.6561 0.5585 0.3284 5.5875 0.2760
+DEAL::0.5712 0.6221 0.6072 0.2760 5.8955
+DEAL::
+DEAL::Cholesky Decomposition Matrix:
+DEAL::==============================
+DEAL::5x5 matrix
+DEAL::2.3635 0 0 0 0
+DEAL::0.2213 2.3539 0 0 0
+DEAL::0.0634 0.2337 2.3422 0 0
+DEAL::0.2776 0.2112 0.1116 2.3353 0
+DEAL::0.2417 0.2416 0.2286 0.0567 2.3923
+DEAL::
+DEAL::Transposed Cholesky Decomposition Matrix:
+DEAL::=========================================
+DEAL::5x5 matrix
+DEAL::2.3635 0.2213 0.0634 0.2776 0.2417
+DEAL::0 2.3539 0.2337 0.2112 0.2416
+DEAL::0 0 2.3422 0.1116 0.2286
+DEAL::0 0 0 2.3353 0.0567
+DEAL::0 0 0 0 2.3923
+DEAL::
+DEAL::Reconstructed Matrix:
+DEAL::=====================
+DEAL::5x5 matrix
+DEAL::5.5862 0.5230 0.1500 0.6561 0.5712
+DEAL::0.5230 5.5897 0.5640 0.5585 0.6221
+DEAL::0.1500 0.5640 5.5445 0.3284 0.6072
+DEAL::0.6561 0.5585 0.3284 5.5875 0.2760
+DEAL::0.5712 0.6221 0.6072 0.2760 5.8955
+DEAL::
+DEAL::the SAME to 1e-10 tolerance
+DEAL::
+DEAL::
+DEAL::Original Matrix:
+DEAL::================
+DEAL::2x2 matrix
+DEAL::0.8147 0.5164
+DEAL::0.5164 0.9134
+DEAL::
+DEAL::Cholesky Decomposition Matrix:
+DEAL::==============================
+DEAL::2x2 matrix
+DEAL::0.9026 0
+DEAL::0.5721 0.7656
+DEAL::
+DEAL::Transposed Cholesky Decomposition Matrix:
+DEAL::=========================================
+DEAL::2x2 matrix
+DEAL::0.9026 0.5721
+DEAL::0 0.7656
+DEAL::
+DEAL::Reconstructed Matrix:
+DEAL::=====================
+DEAL::2x2 matrix
+DEAL::0.8147 0.5164
+DEAL::0.5164 0.9134
+DEAL::
+DEAL::the SAME to 1e-10 tolerance
+DEAL::
+DEAL::
+DEAL::Original Matrix:
+DEAL::================
+DEAL::3x3 matrix
+DEAL::1.8086 0.8459 0.6671
+DEAL::0.8459 1.3989 0.2817
+DEAL::0.6671 0.2817 0.7418
+DEAL::
+DEAL::Cholesky Decomposition Matrix:
+DEAL::==============================
+DEAL::3x3 matrix
+DEAL::1.3449 0 0
+DEAL::0.6290 1.0016 0
+DEAL::0.4960 -0.0302 0.7034
+DEAL::
+DEAL::Transposed Cholesky Decomposition Matrix:
+DEAL::=========================================
+DEAL::3x3 matrix
+DEAL::1.3449 0.6290 0.4960
+DEAL::0 1.0016 -0.0302
+DEAL::0 0 0.7034
+DEAL::
+DEAL::Reconstructed Matrix:
+DEAL::=====================
+DEAL::3x3 matrix
+DEAL::1.8086 0.8459 0.6671
+DEAL::0.8459 1.3989 0.2817
+DEAL::0.6671 0.2817 0.7418
+DEAL::
+DEAL::the SAME to 1e-10 tolerance
+DEAL::
+DEAL::
+DEAL::Original Matrix:
+DEAL::================
+DEAL::5x5 matrix
+DEAL::5.5862 0.5230 0.1500 0.6561 0.5712
+DEAL::0.5230 5.5897 0.5640 0.5585 0.6221
+DEAL::0.1500 0.5640 5.5445 0.3284 0.6072
+DEAL::0.6561 0.5585 0.3284 5.5875 0.2760
+DEAL::0.5712 0.6221 0.6072 0.2760 5.8955
+DEAL::
+DEAL::Cholesky Decomposition Matrix:
+DEAL::==============================
+DEAL::5x5 matrix
+DEAL::2.3635 0 0 0 0
+DEAL::0.2213 2.3539 0 0 0
+DEAL::0.0634 0.2337 2.3422 0 0
+DEAL::0.2776 0.2112 0.1116 2.3353 0
+DEAL::0.2417 0.2416 0.2286 0.0567 2.3923
+DEAL::
+DEAL::Transposed Cholesky Decomposition Matrix:
+DEAL::=========================================
+DEAL::5x5 matrix
+DEAL::2.3635 0.2213 0.0634 0.2776 0.2417
+DEAL::0 2.3539 0.2337 0.2112 0.2416
+DEAL::0 0 2.3422 0.1116 0.2286
+DEAL::0 0 0 2.3353 0.0567
+DEAL::0 0 0 0 2.3923
+DEAL::
+DEAL::Reconstructed Matrix:
+DEAL::=====================
+DEAL::5x5 matrix
+DEAL::5.5862 0.5230 0.1500 0.6561 0.5712
+DEAL::0.5230 5.5897 0.5640 0.5585 0.6221
+DEAL::0.1500 0.5640 5.5445 0.3284 0.6072
+DEAL::0.6561 0.5585 0.3284 5.5875 0.2760
+DEAL::0.5712 0.6221 0.6072 0.2760 5.8955
+DEAL::
+DEAL::the SAME to 1e-10 tolerance
+DEAL::
+DEAL::
--- /dev/null
+//----------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2007, 2009 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.
+//
+//----------------------------------------------------------------------
+
+
+// check FullMatrix::outer_product for correct functionality
+
+
+#include "../tests.h"
+#include "full_matrix_common.h"
+
+
+std::string output_file_name = "full_matrix_57/output";
+
+const double ints[9] = { 0, -1., 1., -2., 2., -3., 3., -4., 4.};
+
+template <typename number>
+void
+check ()
+ {
+ deallog << std::fixed;
+ deallog << std::setprecision(1);
+ deallog.depth_console(0);
+
+ Vector<double> First4(4), Second4(4), First9(9);
+ FullMatrix<double>F4(4,4), F9(9,9);
+
+
+ deallog << "Original Vector V" << std::endl;
+ deallog << "=================" << std::endl;
+
+ for (unsigned int i = 0; i<First4.size(); i++)
+ {
+ First4(i) = ints[i];
+ deallog << First4(i) << " ";
+ }
+ deallog << std::endl;
+ deallog << std::endl;
+
+ deallog << "Original Vector W" << std::endl;
+ deallog << "=================" << std::endl;
+
+ for (unsigned int i = 0; i<Second4.size(); i++)
+ {
+ Second4(i) = ints[i+1];
+ deallog << Second4(i) << " ";
+ }
+ deallog << std::endl;
+ deallog << std::endl;
+
+ F4.outer_product(First4, Second4);
+
+ deallog << "Outer_Product of V and W" << std::endl;
+ deallog << "========================" << std::endl;
+ display_matrix(F4);
+ deallog << std::endl;
+ deallog << std::endl;
+
+ deallog << std::endl;
+
+ F4.outer_product(First4, First4);
+
+ deallog << "Outer_Product of V and V" << std::endl;
+ deallog << "========================" << std::endl;
+ display_matrix(F4);
+ deallog << std::endl;
+ deallog << std::endl;
+
+
+ F4.outer_product(Second4, Second4);
+
+ deallog << "Outer_Product of W and W" << std::endl;
+ deallog << "========================" << std::endl;
+ display_matrix(F4);
+ deallog << std::endl;
+ deallog << std::endl;
+ deallog << std::endl;
+
+ deallog << "Vector V" << std::endl;
+ deallog << "========" << std::endl;
+
+ for (unsigned int i = 0; i<First9.size(); i++)
+ {
+ First9(i) = ints[i];
+ deallog << First9(i) << " ";
+ }
+ deallog << std::endl;
+ deallog << std::endl;
+
+ F9.outer_product(First9, First9);
+
+ deallog << "Outer_Product of V and V" << std::endl;
+ deallog << "========================" << std::endl;
+ display_matrix(F9);
+ deallog << std::endl;
+ deallog << std::endl;
+
+ deallog << std::endl;
+
+ deallog << "Vector V" << std::endl;
+ deallog << "========" << std::endl;
+
+ for (unsigned int i = 0; i<First9.size(); i++)
+ {
+ First9(i) = ints[8-i];
+ deallog << First9(i) << " ";
+ }
+ deallog << std::endl;
+ deallog << std::endl;
+
+ F9.outer_product(First9, First9);
+
+ deallog << "Outer_Product of V and V" << std::endl;
+ deallog << "========================" << std::endl;
+ display_matrix(F9);
+ deallog << std::endl;
+ deallog << std::endl;
+
+ deallog << std::endl;
+}
+
+
--- /dev/null
+
+DEAL::Original Vector V
+DEAL::=================
+DEAL::0 -1.0 1.0 -2.0
+DEAL::
+DEAL::Original Vector W
+DEAL::=================
+DEAL::-1.0 1.0 -2.0 2.0
+DEAL::
+DEAL::Outer_Product of V and W
+DEAL::========================
+DEAL::4x4 matrix
+DEAL::0 0 0 0
+DEAL::1.0 -1.0 2.0 -2.0
+DEAL::-1.0 1.0 -2.0 2.0
+DEAL::2.0 -2.0 4.0 -4.0
+DEAL::
+DEAL::
+DEAL::
+DEAL::Outer_Product of V and V
+DEAL::========================
+DEAL::4x4 matrix
+DEAL::0 0 0 0
+DEAL::0 1.0 -1.0 2.0
+DEAL::0 -1.0 1.0 -2.0
+DEAL::0 2.0 -2.0 4.0
+DEAL::
+DEAL::
+DEAL::Outer_Product of W and W
+DEAL::========================
+DEAL::4x4 matrix
+DEAL::1.0 -1.0 2.0 -2.0
+DEAL::-1.0 1.0 -2.0 2.0
+DEAL::2.0 -2.0 4.0 -4.0
+DEAL::-2.0 2.0 -4.0 4.0
+DEAL::
+DEAL::
+DEAL::
+DEAL::Vector V
+DEAL::========
+DEAL::0 -1.0 1.0 -2.0 2.0 -3.0 3.0 -4.0 4.0
+DEAL::
+DEAL::Outer_Product of V and V
+DEAL::========================
+DEAL::9x9 matrix
+DEAL::0 0 0 0 0 0 0 0 0
+DEAL::0 1.0 -1.0 2.0 -2.0 3.0 -3.0 4.0 -4.0
+DEAL::0 -1.0 1.0 -2.0 2.0 -3.0 3.0 -4.0 4.0
+DEAL::0 2.0 -2.0 4.0 -4.0 6.0 -6.0 8.0 -8.0
+DEAL::0 -2.0 2.0 -4.0 4.0 -6.0 6.0 -8.0 8.0
+DEAL::0 3.0 -3.0 6.0 -6.0 9.0 -9.0 12.0 -12.0
+DEAL::0 -3.0 3.0 -6.0 6.0 -9.0 9.0 -12.0 12.0
+DEAL::0 4.0 -4.0 8.0 -8.0 12.0 -12.0 16.0 -16.0
+DEAL::0 -4.0 4.0 -8.0 8.0 -12.0 12.0 -16.0 16.0
+DEAL::
+DEAL::
+DEAL::
+DEAL::Vector V
+DEAL::========
+DEAL::4.0 -4.0 3.0 -3.0 2.0 -2.0 1.0 -1.0 0
+DEAL::
+DEAL::Outer_Product of V and V
+DEAL::========================
+DEAL::9x9 matrix
+DEAL::16.0 -16.0 12.0 -12.0 8.0 -8.0 4.0 -4.0 0
+DEAL::-16.0 16.0 -12.0 12.0 -8.0 8.0 -4.0 4.0 0
+DEAL::12.0 -12.0 9.0 -9.0 6.0 -6.0 3.0 -3.0 0
+DEAL::-12.0 12.0 -9.0 9.0 -6.0 6.0 -3.0 3.0 0
+DEAL::8.0 -8.0 6.0 -6.0 4.0 -4.0 2.0 -2.0 0
+DEAL::-8.0 8.0 -6.0 6.0 -4.0 4.0 -2.0 2.0 0
+DEAL::4.0 -4.0 3.0 -3.0 2.0 -2.0 1.0 -1.0 0
+DEAL::-4.0 4.0 -3.0 3.0 -2.0 2.0 -1.0 1.0 0
+DEAL::0 0 0 0 0 0 0 0 0
+DEAL::
+DEAL::
+DEAL::
+DEAL::Original Vector V
+DEAL::=================
+DEAL::0 -1.0 1.0 -2.0
+DEAL::
+DEAL::Original Vector W
+DEAL::=================
+DEAL::-1.0 1.0 -2.0 2.0
+DEAL::
+DEAL::Outer_Product of V and W
+DEAL::========================
+DEAL::4x4 matrix
+DEAL::0 0 0 0
+DEAL::1.0 -1.0 2.0 -2.0
+DEAL::-1.0 1.0 -2.0 2.0
+DEAL::2.0 -2.0 4.0 -4.0
+DEAL::
+DEAL::
+DEAL::
+DEAL::Outer_Product of V and V
+DEAL::========================
+DEAL::4x4 matrix
+DEAL::0 0 0 0
+DEAL::0 1.0 -1.0 2.0
+DEAL::0 -1.0 1.0 -2.0
+DEAL::0 2.0 -2.0 4.0
+DEAL::
+DEAL::
+DEAL::Outer_Product of W and W
+DEAL::========================
+DEAL::4x4 matrix
+DEAL::1.0 -1.0 2.0 -2.0
+DEAL::-1.0 1.0 -2.0 2.0
+DEAL::2.0 -2.0 4.0 -4.0
+DEAL::-2.0 2.0 -4.0 4.0
+DEAL::
+DEAL::
+DEAL::
+DEAL::Vector V
+DEAL::========
+DEAL::0 -1.0 1.0 -2.0 2.0 -3.0 3.0 -4.0 4.0
+DEAL::
+DEAL::Outer_Product of V and V
+DEAL::========================
+DEAL::9x9 matrix
+DEAL::0 0 0 0 0 0 0 0 0
+DEAL::0 1.0 -1.0 2.0 -2.0 3.0 -3.0 4.0 -4.0
+DEAL::0 -1.0 1.0 -2.0 2.0 -3.0 3.0 -4.0 4.0
+DEAL::0 2.0 -2.0 4.0 -4.0 6.0 -6.0 8.0 -8.0
+DEAL::0 -2.0 2.0 -4.0 4.0 -6.0 6.0 -8.0 8.0
+DEAL::0 3.0 -3.0 6.0 -6.0 9.0 -9.0 12.0 -12.0
+DEAL::0 -3.0 3.0 -6.0 6.0 -9.0 9.0 -12.0 12.0
+DEAL::0 4.0 -4.0 8.0 -8.0 12.0 -12.0 16.0 -16.0
+DEAL::0 -4.0 4.0 -8.0 8.0 -12.0 12.0 -16.0 16.0
+DEAL::
+DEAL::
+DEAL::
+DEAL::Vector V
+DEAL::========
+DEAL::4.0 -4.0 3.0 -3.0 2.0 -2.0 1.0 -1.0 0
+DEAL::
+DEAL::Outer_Product of V and V
+DEAL::========================
+DEAL::9x9 matrix
+DEAL::16.0 -16.0 12.0 -12.0 8.0 -8.0 4.0 -4.0 0
+DEAL::-16.0 16.0 -12.0 12.0 -8.0 8.0 -4.0 4.0 0
+DEAL::12.0 -12.0 9.0 -9.0 6.0 -6.0 3.0 -3.0 0
+DEAL::-12.0 12.0 -9.0 9.0 -6.0 6.0 -3.0 3.0 0
+DEAL::8.0 -8.0 6.0 -6.0 4.0 -4.0 2.0 -2.0 0
+DEAL::-8.0 8.0 -6.0 6.0 -4.0 4.0 -2.0 2.0 0
+DEAL::4.0 -4.0 3.0 -3.0 2.0 -2.0 1.0 -1.0 0
+DEAL::-4.0 4.0 -3.0 3.0 -2.0 2.0 -1.0 1.0 0
+DEAL::0 0 0 0 0 0 0 0 0
+DEAL::
+DEAL::
+DEAL::