00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__tridiagonal_matrix_h
00013 #define __deal2__tridiagonal_matrix_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/subscriptor.h>
00017 #include <deal.II/lac/lapack_support.h>
00018
00019 #include <vector>
00020 #include <iomanip>
00021
00022 DEAL_II_NAMESPACE_OPEN
00023
00024
00025 template<typename number> class Vector;
00026
00027
00044 template<typename number>
00045 class TridiagonalMatrix
00046 {
00047 public:
00056 TridiagonalMatrix(unsigned int n = 0,
00057 bool symmetric = false);
00058
00065 void reinit(unsigned int n,
00066 bool symmetric = false);
00067
00068
00070
00071
00072
00078 unsigned int m () const;
00079
00085 unsigned int n () const;
00086
00097 bool all_zero () const;
00098
00099
00100
00102
00103
00104
00110 number operator()(unsigned int i, unsigned int j) const;
00111
00127 number& operator()(unsigned int i, unsigned int j);
00128
00130
00131
00132
00146 void vmult (Vector<number> &w,
00147 const Vector<number> &v,
00148 const bool adding=false) const;
00149
00161 void vmult_add (Vector<number> &w,
00162 const Vector<number> &v) const;
00163
00179 void Tvmult (Vector<number> &w,
00180 const Vector<number> &v,
00181 const bool adding=false) const;
00182
00194 void Tvmult_add (Vector<number> &w,
00195 const Vector<number> &v) const;
00196
00204 number matrix_scalar_product (const Vector<number> &u,
00205 const Vector<number> &v) const;
00206
00226 number matrix_norm_square (const Vector<number> &v) const;
00227
00229
00230
00231
00242 number l1_norm () const;
00243
00255 number linfty_norm () const;
00256
00262 number frobenius_norm () const;
00263
00277 number relative_symmetry_norm2 () const;
00279
00280
00281
00291 void compute_eigenvalues();
00297 number eigenvalue(const unsigned i) const;
00299
00300
00301
00305 template <class OUT>
00306 void print (OUT& s,
00307 const unsigned int width=5,
00308 const unsigned int precision=2) const;
00309
00315 std::size_t memory_consumption () const;
00317
00318 private:
00322 std::vector<number> diagonal;
00338 std::vector<number> left;
00348 std::vector<number> right;
00349
00356 bool is_symmetric;
00357
00369 LAPACKSupport::State state;
00370 };
00371
00374
00375 #ifndef DOXYGEN
00376
00377 template<typename number>
00378 unsigned int
00379 TridiagonalMatrix<number>::m() const
00380 {
00381 return diagonal.size();
00382 }
00383
00384
00385
00386 template<typename number>
00387 unsigned int
00388 TridiagonalMatrix<number>::n() const
00389 {
00390 return diagonal.size();
00391 }
00392
00393
00394 template<typename number>
00395 inline
00396 number
00397 TridiagonalMatrix<number>::operator()(unsigned int i, unsigned int j) const
00398 {
00399 Assert(i<n(), ExcIndexRange(i,0,n()));
00400 Assert(j<n(), ExcIndexRange(j,0,n()));
00401 Assert (i<=j+1, ExcIndexRange(i,j-1,j+2));
00402 Assert (j<=i+1, ExcIndexRange(j,i-1,i+2));
00403
00404 if (j==i)
00405 return diagonal[i];
00406 if (j==i-1)
00407 {
00408 if (is_symmetric)
00409 return right[i-1];
00410 else
00411 return left[i];
00412 }
00413
00414 if (j==i+1)
00415 return right[i];
00416
00417 Assert (false, ExcInternalError());
00418 return 0;
00419 }
00420
00421
00422 template<typename number>
00423 inline
00424 number&
00425 TridiagonalMatrix<number>::operator()(unsigned int i, unsigned int j)
00426 {
00427 Assert(i<n(), ExcIndexRange(i,0,n()));
00428 Assert(j<n(), ExcIndexRange(j,0,n()));
00429 Assert (i<=j+1, ExcIndexRange(i,j-1,j+2));
00430 Assert (j<=i+1, ExcIndexRange(j,i-1,i+2));
00431
00432 if (j==i)
00433 return diagonal[i];
00434 if (j==i-1)
00435 {
00436 if (is_symmetric)
00437 return right[i-1];
00438 else
00439 return left[i];
00440 }
00441
00442 if (j==i+1)
00443 return right[i];
00444
00445 Assert (false, ExcInternalError());
00446 return diagonal[0];
00447 }
00448
00449
00450 template <typename number>
00451 template <class OUT>
00452 void
00453 TridiagonalMatrix<number>::print (
00454 OUT& s,
00455 const unsigned int width,
00456 const unsigned int) const
00457 {
00458 for (unsigned int i=0;i<n();++i)
00459 {
00460 if (i>0)
00461 s << std::setw(width) << (*this)(i,i-1);
00462 else
00463 s << std::setw(width) << "";
00464
00465 s << ' ' << (*this)(i,i) << ' ';
00466
00467 if (i<n()-1)
00468 s << std::setw(width) << (*this)(i,i+1);
00469
00470 s << std::endl;
00471 }
00472 }
00473
00474
00475 #endif // DOXYGEN
00476
00477 DEAL_II_NAMESPACE_CLOSE
00478
00479 #endif
00480
00481