From e4769eafd65198e0414faa5c961121a3faf322f0 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 25 Aug 2000 19:31:26 +0000 Subject: [PATCH] check for FullMatrix::invert git-svn-id: https://svn.dealii.org/trunk@3275 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/lac/full_matrix.cc | 183 ++++++++++++++++++++++++--------------- 1 file changed, 113 insertions(+), 70 deletions(-) diff --git a/tests/lac/full_matrix.cc b/tests/lac/full_matrix.cc index 3aa2f3f1fc..ab77c1752c 100644 --- a/tests/lac/full_matrix.cc +++ b/tests/lac/full_matrix.cc @@ -15,12 +15,26 @@ #include #include #include +#include #include #include #include #include +// Create a positive definite random matrix + +void random_matrix(FullMatrix& A) +{ + for (unsigned int i=0; i 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 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 x(i); + Vector y(i); + Vector 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 u(5); - GrowingVectorMemory > mem; - - SolverControl control (500,1.e-8, false, true); - - if (true) - { - u = 1.; - EigenPower > - von_Mises(control, mem, 0.); - double eigen = 0.; - von_Mises.solve(eigen, A, u); - deallog << "Eigenvalue: " << eigen << endl; - } - if (true) - { - u = 1.; - EigenPower > - 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 > - von_Mises(control, mem, 0.); - double eigen = 0.; - von_Mises.solve(eigen, H, u); - deallog << "Eigenvalue: " << eigen << endl; + FullMatrix 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 u(5); + GrowingVectorMemory > mem; + + SolverControl control (500,1.e-8, false, true); + + if (true) + { + u = 1.; + EigenPower > + von_Mises(control, mem, 0.); + double eigen = 0.; + von_Mises.solve(eigen, A, u); + deallog << "Eigenvalue: " << eigen << endl; + } + if (true) + { + u = 1.; + EigenPower > + 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 > + von_Mises(control, mem, 0.); + double eigen = 0.; + von_Mises.solve(eigen, H, u); + deallog << "Eigenvalue: " << eigen << endl; + } + if (true) + { + u = 1.; + EigenPower > + von_Mises(control, mem, -4.); + double eigen = 0.; + von_Mises.solve(eigen, H, u); + deallog << "Eigenvalue: " << eigen << endl; + } } - if (true) - { - u = 1.; - EigenPower > - von_Mises(control, mem, -4.); - double eigen = 0.; - von_Mises.solve(eigen, H, u); - deallog << "Eigenvalue: " << eigen << endl; - } } -- 2.39.5