00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__precondition_block_base_h
00013 #define __deal2__precondition_block_base_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/subscriptor.h>
00019 #include <deal.II/base/smartpointer.h>
00020 #include <deal.II/base/memory_consumption.h>
00021 #include <deal.II/lac/householder.h>
00022 #include <deal.II/lac/lapack_full_matrix.h>
00023
00024 #include <vector>
00025
00026 DEAL_II_NAMESPACE_OPEN
00027
00028 template <typename number> class FullMatrix;
00029 template <typename number> class Vector;
00030
00050 template <typename number>
00051 class PreconditionBlockBase
00052 {
00053 public:
00059 enum Inversion
00060 {
00066 gauss_jordan,
00071 householder,
00076 svd
00077 };
00078
00083 PreconditionBlockBase(bool store_diagonals = false,
00084 Inversion method = gauss_jordan);
00085
00089 ~PreconditionBlockBase();
00090
00098 void clear();
00099
00108 void reinit(unsigned int nblocks, unsigned int blocksize, bool compress,
00109 Inversion method = gauss_jordan);
00110
00115 void inverses_computed(bool are_they);
00116
00130 void set_same_diagonal ();
00131
00136 bool same_diagonal () const;
00137
00143 bool store_diagonals() const;
00144
00149 bool inverses_ready () const;
00150
00154 bool empty () const;
00155
00159 unsigned int size() const;
00160
00167 number el(unsigned int i, unsigned int j) const;
00168
00173 template <typename number2>
00174 void inverse_vmult(unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const;
00175
00180 template <typename number2>
00181 void inverse_Tvmult(unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const;
00182
00187 FullMatrix<number>& inverse (unsigned int i);
00188
00193 Householder<number>& inverse_householder (unsigned int i);
00194
00199 LAPACKFullMatrix<number>& inverse_svd (unsigned int i);
00200
00205 const FullMatrix<number>& inverse (unsigned int i) const;
00206
00211 const Householder<number>& inverse_householder (unsigned int i) const;
00212
00217 const LAPACKFullMatrix<number>& inverse_svd (unsigned int i) const;
00218
00223 FullMatrix<number>& diagonal (unsigned int i);
00224
00229 const FullMatrix<number>& diagonal (unsigned int i) const;
00230
00240 void log_statistics () const;
00241
00247 std::size_t memory_consumption () const;
00248
00255 DeclException0 (ExcDiagonalsNotStored);
00256
00264 DeclException0 (ExcInverseNotAvailable);
00265
00266 protected:
00270 Inversion inversion;
00271
00272 private:
00278 unsigned int n_diagonal_blocks;
00279
00292 std::vector<FullMatrix<number> > var_inverse_full;
00293
00306 std::vector<Householder<number> > var_inverse_householder;
00307
00320 std::vector<LAPACKFullMatrix<number> > var_inverse_svd;
00321
00327 std::vector<FullMatrix<number> > var_diagonal;
00328
00329
00334 bool var_store_diagonals;
00335
00340 bool var_same_diagonal;
00341
00347 bool var_inverses_ready;
00348 };
00349
00350
00351
00352 template <typename number>
00353 inline
00354 PreconditionBlockBase<number>::PreconditionBlockBase(
00355 bool store, Inversion method)
00356 :
00357 inversion(method),
00358 n_diagonal_blocks(0),
00359 var_store_diagonals(store),
00360 var_same_diagonal(false),
00361 var_inverses_ready(false)
00362 {}
00363
00364
00365 template <typename number>
00366 inline
00367 void
00368 PreconditionBlockBase<number>::clear()
00369 {
00370 if (var_inverse_full.size()!=0)
00371 var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
00372 if (var_inverse_householder.size()!=0)
00373 var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end());
00374 if (var_inverse_svd.size()!=0)
00375 var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
00376 if (var_diagonal.size()!=0)
00377 var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
00378 var_same_diagonal = false;
00379 var_inverses_ready = false;
00380 n_diagonal_blocks = 0;
00381 }
00382
00383 template <typename number>
00384 inline
00385 void
00386 PreconditionBlockBase<number>::reinit(unsigned int n, unsigned int b, bool compress,
00387 Inversion method)
00388 {
00389 inversion = method;
00390 var_same_diagonal = compress;
00391 var_inverses_ready = false;
00392 n_diagonal_blocks = n;
00393
00394 if (compress)
00395 {
00396 switch(inversion)
00397 {
00398 case gauss_jordan:
00399 var_inverse_full.resize(1);
00400 var_inverse_full[0].reinit(b,b);
00401 break;
00402 case householder:
00403 var_inverse_householder.resize(1);
00404 break;
00405 case svd:
00406 var_inverse_svd.resize(1);
00407 var_inverse_svd[0].reinit(b,b);
00408 break;
00409 default:
00410 Assert(false, ExcNotImplemented());
00411 }
00412
00413 if (store_diagonals())
00414 {
00415 var_diagonal.resize(1);
00416 var_diagonal[0].reinit(b,b);
00417 }
00418 }
00419 else
00420 {
00421
00422
00423
00424
00425
00426
00427
00428
00429 if (store_diagonals())
00430 {
00431 std::vector<FullMatrix<number> >
00432 tmp(n, FullMatrix<number>(b));
00433 var_diagonal.swap (tmp);
00434 }
00435
00436 switch(inversion)
00437 {
00438 case gauss_jordan:
00439 {
00440 std::vector<FullMatrix<number> >
00441 tmp(n, FullMatrix<number>(b));
00442 var_inverse_full.swap (tmp);
00443 break;
00444 }
00445 case householder:
00446 var_inverse_householder.resize(n);
00447 break;
00448 case svd:
00449 {
00450 std::vector<LAPACKFullMatrix<number> >
00451 tmp(n, LAPACKFullMatrix<number>(b));
00452 var_inverse_svd.swap (tmp);
00453 break;
00454 }
00455 default:
00456 Assert(false, ExcNotImplemented());
00457 }
00458 }
00459 }
00460
00461
00462 template <typename number>
00463 inline
00464 unsigned int
00465 PreconditionBlockBase<number>::size() const
00466 {
00467 return n_diagonal_blocks;
00468 }
00469
00470 template <typename number>
00471 template <typename number2>
00472 inline
00473 void
00474 PreconditionBlockBase<number>::inverse_vmult(
00475 unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const
00476 {
00477 const unsigned int ii = same_diagonal() ? 0U : i;
00478
00479 switch (inversion)
00480 {
00481 case gauss_jordan:
00482 AssertIndexRange (ii, var_inverse_full.size());
00483 var_inverse_full[ii].vmult(dst, src);
00484 break;
00485 case householder:
00486 AssertIndexRange (ii, var_inverse_householder.size());
00487 var_inverse_householder[ii].vmult(dst, src);
00488 break;
00489 case svd:
00490 AssertIndexRange (ii, var_inverse_svd.size());
00491 var_inverse_svd[ii].vmult(dst, src);
00492 break;
00493 default:
00494 Assert(false, ExcNotImplemented());
00495 }
00496 }
00497
00498
00499 template <typename number>
00500 template <typename number2>
00501 inline
00502 void
00503 PreconditionBlockBase<number>::inverse_Tvmult(
00504 unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const
00505 {
00506 const unsigned int ii = same_diagonal() ? 0U : i;
00507
00508 switch (inversion)
00509 {
00510 case gauss_jordan:
00511 AssertIndexRange (ii, var_inverse_full.size());
00512 var_inverse_full[ii].Tvmult(dst, src);
00513 break;
00514 case householder:
00515 AssertIndexRange (ii, var_inverse_householder.size());
00516 var_inverse_householder[ii].Tvmult(dst, src);
00517 break;
00518 case svd:
00519 AssertIndexRange (ii, var_inverse_svd.size());
00520 var_inverse_svd[ii].Tvmult(dst, src);
00521 break;
00522 default:
00523 Assert(false, ExcNotImplemented());
00524 }
00525 }
00526
00527
00528 template <typename number>
00529 inline
00530 const FullMatrix<number>&
00531 PreconditionBlockBase<number>::inverse(unsigned int i) const
00532 {
00533 if (same_diagonal())
00534 return var_inverse_full[0];
00535
00536 Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
00537 return var_inverse_full[i];
00538 }
00539
00540
00541 template <typename number>
00542 inline
00543 const Householder<number>&
00544 PreconditionBlockBase<number>::inverse_householder(unsigned int i) const
00545 {
00546 if (same_diagonal())
00547 return var_inverse_householder[0];
00548
00549 AssertIndexRange (i, var_inverse_householder.size());
00550 return var_inverse_householder[i];
00551 }
00552
00553
00554 template <typename number>
00555 inline
00556 const LAPACKFullMatrix<number>&
00557 PreconditionBlockBase<number>::inverse_svd(unsigned int i) const
00558 {
00559 if (same_diagonal())
00560 return var_inverse_svd[0];
00561
00562 AssertIndexRange (i, var_inverse_svd.size());
00563 return var_inverse_svd[i];
00564 }
00565
00566
00567 template <typename number>
00568 inline
00569 const FullMatrix<number>&
00570 PreconditionBlockBase<number>::diagonal(unsigned int i) const
00571 {
00572 Assert(store_diagonals(), ExcDiagonalsNotStored());
00573
00574 if (same_diagonal())
00575 return var_diagonal[0];
00576
00577 Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
00578 return var_diagonal[i];
00579 }
00580
00581
00582 template <typename number>
00583 inline
00584 FullMatrix<number>&
00585 PreconditionBlockBase<number>::inverse(unsigned int i)
00586 {
00587 Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());
00588
00589 if (same_diagonal())
00590 return var_inverse_full[0];
00591
00592 Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
00593 return var_inverse_full[i];
00594 }
00595
00596
00597 template <typename number>
00598 inline
00599 Householder<number>&
00600 PreconditionBlockBase<number>::inverse_householder(unsigned int i)
00601 {
00602 Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());
00603
00604 if (same_diagonal())
00605 return var_inverse_householder[0];
00606
00607 AssertIndexRange (i, var_inverse_householder.size());
00608 return var_inverse_householder[i];
00609 }
00610
00611
00612 template <typename number>
00613 inline
00614 LAPACKFullMatrix<number>&
00615 PreconditionBlockBase<number>::inverse_svd(unsigned int i)
00616 {
00617 Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());
00618
00619 if (same_diagonal())
00620 return var_inverse_svd[0];
00621
00622 AssertIndexRange (i, var_inverse_svd.size());
00623 return var_inverse_svd[i];
00624 }
00625
00626
00627 template <typename number>
00628 inline
00629 FullMatrix<number>&
00630 PreconditionBlockBase<number>::diagonal(unsigned int i)
00631 {
00632 Assert(store_diagonals(), ExcDiagonalsNotStored());
00633
00634 if (same_diagonal())
00635 return var_diagonal[0];
00636
00637 Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
00638 return var_diagonal[i];
00639 }
00640
00641
00642 template <typename number>
00643 inline bool
00644 PreconditionBlockBase<number>::same_diagonal() const
00645 {
00646 return var_same_diagonal;
00647 }
00648
00649
00650 template <typename number>
00651 inline bool
00652 PreconditionBlockBase<number>::store_diagonals() const
00653 {
00654 return var_store_diagonals;
00655 }
00656
00657
00658 template <typename number>
00659 inline void
00660 PreconditionBlockBase<number>::inverses_computed(bool x)
00661 {
00662 var_inverses_ready = x;
00663 }
00664
00665
00666 template <typename number>
00667 inline bool
00668 PreconditionBlockBase<number>::inverses_ready() const
00669 {
00670 return var_inverses_ready;
00671 }
00672
00673
00674 template <typename number>
00675 inline void
00676 PreconditionBlockBase<number>::log_statistics () const
00677 {
00678 deallog << "PreconditionBlockBase: " << size() << " blocks; ";
00679
00680 if (inversion == svd)
00681 {
00682 unsigned int kermin = 100000000, kermax = 0;
00683 double sigmin = 1.e300, sigmax= -1.e300;
00684 double kappamin = 1.e300, kappamax= -1.e300;
00685
00686 for (unsigned int b=0;b<size();++b)
00687 {
00688 const LAPACKFullMatrix<number>& matrix = inverse_svd(b);
00689 unsigned int k=1;
00690 while (k <= matrix.n_cols() && matrix.singular_value(matrix.n_cols()-k) == 0)
00691 ++k;
00692 const double s0 = matrix.singular_value(0);
00693 const double sm = matrix.singular_value(matrix.n_cols()-k);
00694 const double co = sm/s0;
00695
00696 if (kermin > k) kermin = k-1;
00697 if (kermax < k) kermax = k-1;
00698 if (s0 < sigmin) sigmin = s0;
00699 if (sm > sigmax) sigmax = sm;
00700 if (co < kappamin) kappamin = co;
00701 if (co > kappamax) kappamax = co;
00702 }
00703 deallog << "dim ker [" << kermin << ':' << kermax
00704 << "] sigma [" << sigmin << ':' << sigmax
00705 << "] kappa [" << kappamin << ':' << kappamax << ']' << std::endl;
00706
00707 }
00708 else if (inversion == householder)
00709 {
00710 }
00711 else if (inversion == gauss_jordan)
00712 {
00713 }
00714 else
00715 {
00716 Assert(false, ExcNotImplemented());
00717 }
00718 }
00719
00720
00721 template <typename number>
00722 inline
00723 std::size_t
00724 PreconditionBlockBase<number>::memory_consumption () const
00725 {
00726 std::size_t mem = sizeof(*this);
00727 for (unsigned int i=0; i<var_inverse_full.size(); ++i)
00728 mem += MemoryConsumption::memory_consumption(var_inverse_full[i]);
00729 for (unsigned int i=0; i<var_diagonal.size(); ++i)
00730 mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
00731 return mem;
00732 }
00733
00734
00735 DEAL_II_NAMESPACE_CLOSE
00736
00737 #endif
00738