#include <fstream>
#include <iostream>
+#include <iomanip>
#include <cmath>
+
+// Symmetric matrix with constant diagonals [-1,2,-1]
template<typename number>
void
-check_vmult()
+matrix1(TridiagonalMatrix<number>& M)
{
- TridiagonalMatrix<number> M(4);
- Vector<number> u(4);
- Vector<number> v(4);
+ const unsigned int n = M.n();
+ for (unsigned int i=0;i<n;++i)
+ {
+ M(i,i) = 2.;
+ if (i>0) M(i,i-1) = -1.;
+ if (i<n-1) M(i,i+1) = -1.;
+ }
+}
- for (unsigned int i=0;i<4;++i)
+
+// Nonsymmetric matrix with constant diagonals [-1,2,-3]
+template<typename number>
+void
+matrix2(TridiagonalMatrix<number>& M)
+{
+ const unsigned int n = M.n();
+ for (unsigned int i=0;i<n;++i)
{
- u(i) = i+1;
- M(i,i) = i+1;
- if (i>0)
- M(i,i-1) = 0.-i;
- if (i<3)
- M(i,i+1) = 4.-i;
+ M(i,i) = 2.;
+ if (i>0) M(i,i-1) = -1.;
+ if (i<n-1) M(i,i+1) = -3.;
}
+}
+
+
+// Matrix with increasing diagonals and symmetric rows
+template<typename number>
+void
+matrix3(TridiagonalMatrix<number>& M)
+{
+ const unsigned int n = M.n();
+ for (unsigned int i=0;i<n;++i)
+ {
+ M(i,i) = 2.+i;
+ if (i>0) M(i,i-1) = -1.-i;
+ if (i<n-1) M(i,i+1) = -1.-i;
+ }
+}
+
+
+// Symmetric matrix with increasing diagonals
+template<typename number>
+void
+matrix4(TridiagonalMatrix<number>& M)
+{
+ const unsigned int n = M.n();
+ for (unsigned int i=0;i<n;++i)
+ {
+ M(i,i) = 2.+i;
+ if (i>0) M(i,i-1) = 0.-i;
+ if (i<n-1) M(i,i+1) = -1.-i;
+ }
+}
+
+
+// Nonsymmetric matrix with increasing diagonals
+template<typename number>
+void
+matrix5(TridiagonalMatrix<number>& M)
+{
+ const unsigned int n = M.n();
+ for (unsigned int i=0;i<n;++i)
+ {
+ M(i,i) = 2.+i;
+ if (i>0) M(i,i-1) = -1.-i;
+ if (i<n-1) M(i,i+1) = 5.-2.*i;
+ }
+}
+
+
+template<typename number>
+void
+check_vmult(TridiagonalMatrix<number>& M)
+{
+ const unsigned int n = M.n();
+ Vector<number> u(n);
+ Vector<number> v(n);
+
+ for (unsigned int i=0;i<n;++i)
+ u(i) = 1 << i;
+
+ for (unsigned int i=0;i<u.size();++i)
+ deallog << ' ' << std::setw(6) << u(i);
+ deallog << std::endl;
+
M.vmult(v,u);
for (unsigned int i=0;i<v.size();++i)
- deallog << ' ' << v(i);
+ deallog << ' ' << std::setw(6) << v(i);
deallog << std::endl;
-
+
M.vmult_add(v,u);
for (unsigned int i=0;i<v.size();++i)
- deallog << ' ' << v(i);
+ deallog << ' ' << std::setw(6) << v(i);
deallog << std::endl;
-
+
M.Tvmult(v,u);
for (unsigned int i=0;i<v.size();++i)
- deallog << ' ' << v(i);
+ deallog << ' ' << std::setw(6) << v(i);
deallog << std::endl;
-
+
M.Tvmult_add(v,u);
for (unsigned int i=0;i<v.size();++i)
- deallog << ' ' << v(i);
+ deallog << ' ' << std::setw(6) << v(i);
deallog << std::endl;
}
+template<typename number>
+void
+check(unsigned int size)
+{
+ TridiagonalMatrix<number>M(size);
+ matrix1(M);
+ deallog << "Matrix [-1,2,-1]" << std::endl;
+ M.print(deallog, 7, 0);
+ check_vmult(M);
+ matrix2(M);
+ deallog << "Matrix [-1,2,-3]" << std::endl;
+ M.print(deallog, 7, 0);
+ check_vmult(M);
+ matrix3(M);
+ deallog << "Matrix increasing symmetric rows" << std::endl;
+ M.print(deallog, 7, 0);
+ check_vmult(M);
+ matrix4(M);
+ deallog << "Matrix increasing symmetric" << std::endl;
+ M.print(deallog, 7, 0);
+ check_vmult(M);
+ matrix5(M);
+ deallog << "Matrix increasing nonsymmetric" << std::endl;
+ M.print(deallog, 7, 0);
+ check_vmult(M);
+}
+
int main()
{
std::ofstream logfile("tridiagonal_matrix/output");
deallog.depth_console(0);
deallog.threshold_double(1.e-10);
- check_vmult<double>();
- check_vmult<float>();
+ check<double>(5);
+ check<float>(4);
}
-DEAL:: 9. 12. 13. 7.
-DEAL:: 18. 24. 26. 14.
-DEAL:: -1. 2. 3. 22.
-DEAL:: -2. 4. 6. 44.
-DEAL:: 9. 12. 13. 7.
-DEAL:: 18. 24. 26. 14.
-DEAL:: -1. 2. 3. 22.
-DEAL:: -2. 4. 6. 44.
+DEAL::Matrix [-1,2,-1]
+DEAL:: 2. -1.
+DEAL::-1. 2. -1.
+DEAL::-1. 2. -1.
+DEAL::-1. 2. -1.
+DEAL::-1. 2.
+DEAL:: 1. 2. 4. 8. 16.
+DEAL:: 0 -1. -2. -4. 24.
+DEAL:: 0 -2. -4. -8. 48.
+DEAL:: 0 -1. -2. -4. 24.
+DEAL:: 0 -2. -4. -8. 48.
+DEAL::Matrix [-1,2,-3]
+DEAL:: 2. -3.
+DEAL::-1. 2. -3.
+DEAL::-1. 2. -3.
+DEAL::-1. 2. -3.
+DEAL::-1. 2.
+DEAL:: 1. 2. 4. 8. 16.
+DEAL:: -4. -9. -18. -36. 24.
+DEAL:: -8. -18. -36. -72. 48.
+DEAL:: 0 -3. -6. -12. 8.
+DEAL:: 0 -6. -12. -24. 16.
+DEAL::Matrix increasing symmetric rows
+DEAL:: 2. -1.
+DEAL::-2. 3. -2.
+DEAL::-3. 4. -3.
+DEAL::-4. 5. -4.
+DEAL::-5. 6.
+DEAL:: 1. 2. 4. 8. 16.
+DEAL:: 0 -4. -14. -40. 56.
+DEAL:: 0 -8. -28. -80. 112.
+DEAL:: -2. -7. -20. -52. 64.
+DEAL:: -4. -14. -40. -104. 128.
+DEAL::Matrix increasing symmetric
+DEAL:: 2. -1.
+DEAL::-1. 3. -2.
+DEAL::-2. 4. -3.
+DEAL::-3. 5. -4.
+DEAL::-4. 6.
+DEAL:: 1. 2. 4. 8. 16.
+DEAL:: 0 -3. -12. -36. 64.
+DEAL:: 0 -6. -24. -72. 128.
+DEAL:: 0 -3. -12. -36. 64.
+DEAL:: 0 -6. -24. -72. 128.
+DEAL::Matrix increasing nonsymmetric
+DEAL:: 2. 5.
+DEAL::-2. 3. 3.
+DEAL::-3. 4. 1.
+DEAL::-4. 5. -1.
+DEAL::-5. 6.
+DEAL:: 1. 2. 4. 8. 16.
+DEAL:: 12. 16. 18. 8. 56.
+DEAL:: 24. 32. 36. 16. 112.
+DEAL:: -2. -1. -10. -36. 88.
+DEAL:: -4. -2. -20. -72. 176.
+DEAL::Matrix [-1,2,-1]
+DEAL:: 2. -1.
+DEAL::-1. 2. -1.
+DEAL::-1. 2. -1.
+DEAL::-1. 2.
+DEAL:: 1. 2. 4. 8.
+DEAL:: 0. -1. -2. 12.
+DEAL:: 0. -2. -4. 24.
+DEAL:: 0. -1. -2. 12.
+DEAL:: 0. -2. -4. 24.
+DEAL::Matrix [-1,2,-3]
+DEAL:: 2. -3.
+DEAL::-1. 2. -3.
+DEAL::-1. 2. -3.
+DEAL::-1. 2.
+DEAL:: 1. 2. 4. 8.
+DEAL:: -4. -9. -18. 12.
+DEAL:: -8. -18. -36. 24.
+DEAL:: 0. -3. -6. 4.
+DEAL:: 0. -6. -12. 8.
+DEAL::Matrix increasing symmetric rows
+DEAL:: 2. -1.
+DEAL::-2. 3. -2.
+DEAL::-3. 4. -3.
+DEAL::-4. 5.
+DEAL:: 1. 2. 4. 8.
+DEAL:: 0. -4. -14. 24.
+DEAL:: 0. -8. -28. 48.
+DEAL:: -2. -7. -20. 28.
+DEAL:: -4. -14. -40. 56.
+DEAL::Matrix increasing symmetric
+DEAL:: 2. -1.
+DEAL::-1. 3. -2.
+DEAL::-2. 4. -3.
+DEAL::-3. 5.
+DEAL:: 1. 2. 4. 8.
+DEAL:: 0. -3. -12. 28.
+DEAL:: 0. -6. -24. 56.
+DEAL:: 0. -3. -12. 28.
+DEAL:: 0. -6. -24. 56.
+DEAL::Matrix increasing nonsymmetric
+DEAL:: 2. 5.
+DEAL::-2. 3. 3.
+DEAL::-3. 4. 1.
+DEAL::-4. 5.
+DEAL:: 1. 2. 4. 8.
+DEAL:: 12. 16. 18. 24.
+DEAL:: 24. 32. 36. 48.
+DEAL:: -2. -1. -10. 44.
+DEAL:: -4. -2. -20. 88.