include/deal.II/lac/precondition_block_base.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: precondition_block_base.h 25460 2012-04-27 09:19:06Z bangerth @f$
00003 //
00004 //    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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                                        // set the arrays to the right
00422                                        // size. we could do it like this:
00423                                        // var_inverse = vector<>(nblocks,FullMatrix<>())
00424                                        // but this would involve copying many
00425                                        // FullMatrix objects.
00426                                        //
00427                                        // the following is a neat trick which
00428                                        // avoids copying
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3