]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
replace a restriction by a warning in the documentation
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 Dec 2008 21:26:46 +0000 (21:26 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 Dec 2008 21:26:46 +0000 (21:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@17981 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/block_matrix_array.h
deal.II/lac/source/block_matrix_array.cc

index 0131929cf43365ac8bab26de92cff51cb814b085..b3eedb75e90273252f8a42132caa540bd7dc92ba 100644 (file)
@@ -399,10 +399,6 @@ class BlockMatrixArray : public Subscriptor
 
 /*@}*/
 
-/*! @addtogroup Preconditioners
- *@{
- */
-
 
 /**
  * Inversion of a block-triangular matrix.
@@ -416,9 +412,10 @@ class BlockMatrixArray : public Subscriptor
  * performed by the vmult() member function.
  *
  * @note While block indices may be duplicated (see BlockMatrixArray)
- * to add blocks, this is not allowed for diagonal blocks, since
+ * to add blocks, this has to be used with caution, since
  * summing up the inverse of two blocks does not yield the inverse of
- * the sum.
+ * the sum. While the latter would be desirable, we can only perform
+ * the first.
  *
  * The implementation may be a little clumsy, but it should be
  * sufficient as long as the block sizes are much larger than the
@@ -456,6 +453,8 @@ class BlockMatrixArray : public Subscriptor
  * method for solving.
  * @until Error
  *
+ *
+ * @ingroup Preconditioners
  * @author Guido Kanschat, 2001, 2005
  */
 template <typename number = double>
@@ -614,7 +613,6 @@ class BlockTrianglePrecondition
     bool backward;
 };
 
-/*@}*/
 
 #ifndef DOXYGEN
 //---------------------------------------------------------------------------
index 0ee615bcf9bfc06b3ae288ab48d4123220d8000b..4b779c5052e22177fbe36c97bf767c2128ebb007 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2008 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -307,26 +307,32 @@ BlockTrianglePrecondition<number>::do_row (
     m = this->entries.begin();
   typename std::vector<typename BlockMatrixArray<number>::Entry>::const_iterator
     end = this->entries.end();
-  typename std::vector<typename BlockMatrixArray<number>::Entry>::const_iterator
-    diagonal = end;  
+  std::vector<typename std::vector<typename BlockMatrixArray<number>::Entry>::const_iterator>
+    diagonals;
   
   Vector<number>* p_aux = this->mem->alloc();
   Vector<number>& aux = *p_aux;
   
   aux.reinit(dst.block(row_num), true);
-
+  
+                                  // Loop over all entries, since
+                                  // they are not ordered by rows.
   for (; m != end ; ++m)
     {
       const unsigned int i=m->row;
+                                      // Ignore everything not in
+                                      // this row
       if (i != row_num)
        continue;
       const unsigned int j=m->col;
+                                      // Only use the lower (upper)
+                                      // triangle for forward
+                                      // (backward) substitution
       if (((j > i) && !backward) || ((j < i) && backward))
        continue;
       if (j == i)
        {
-         Assert (diagonal == end, ExcMultipleDiagonal(j));
-         diagonal = m;
+         diagonals.push_back(m);
        } else {
          if (m->transpose)
            m->matrix->Tvmult(aux, dst.block(j));
@@ -335,13 +341,37 @@ BlockTrianglePrecondition<number>::do_row (
          dst.block(i).add (-1 * m->prefix, aux);
        }
     }
-  Assert (diagonal != end, ExcNoDiagonal(row_num));
-  
-  if (diagonal->transpose)
-    diagonal->matrix->Tvmult(aux, dst.block(row_num));
+  Assert (diagonals.size() != 0, ExcNoDiagonal(row_num));
+
+                                  // Inverting the diagonal block is
+                                  // simple, if there is only one
+                                  // matrix
+  if (diagonals.size() == 1)
+    {
+      if (diagonals[0]->transpose)
+       diagonals[0]->matrix->Tvmult(aux, dst.block(row_num));
+      else
+       diagonals[0]->matrix->vmult(aux, dst.block(row_num));
+      dst.block(row_num).equ (diagonals[0]->prefix, aux);
+    }
   else
-    diagonal->matrix->vmult(aux, dst.block(row_num));
-  dst.block(row_num).equ (diagonal->prefix, aux);
+    {
+      aux = 0.;
+      for (unsigned int i=0;i<diagonals.size();++i)
+       {
+         m = diagonals[i];
+                                          // First, divide by the current
+                                          // factor, such that we can
+                                          // multiply by it later.
+         aux.scale(1./m->prefix);
+         if (m->transpose)
+           m->matrix->Tvmult_add(aux, dst.block(row_num));
+         else
+           m->matrix->vmult_add(aux, dst.block(row_num));
+         aux.scale(m->prefix);
+       }
+      dst.block(row_num) = aux;
+    }
   
   this->mem->free(p_aux);
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.