// inverse of a matrix by calling an
// iterative solver.
#include <deal.II/lac/iterative_inverse.h>
+#include <lac/schur_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
// rather only comment on
// implementational aspects.
-
- // @sect4{The <code>SchurComplement</code> class template}
-
- // The next class is the Schur
- // complement class. Its rationale
- // has also been discussed in length
- // in the introduction. The only
- // things we would like to note is
- // that the class, too, is derived
- // from the <code>Subscriptor</code> class and
- // that as mentioned above it stores
- // pointers to the entire block
- // matrix and the inverse of the mass
- // matrix block using
- // <code>SmartPointer</code> objects.
- //
- // The <code>vmult</code> function requires
- // two temporary vectors that we do
- // not want to re-allocate and free
- // every time we call this
- // function. Since here, we have full
- // control over the use of these
- // vectors (unlike above, where a
- // class called by the <code>vmult</code>
- // function required these vectors,
- // not the <code>vmult</code> function
- // itself), we allocate them
- // directly, rather than going
- // through the <code>VectorMemory</code>
- // mechanism. However, again, these
- // member variables do not carry any
- // state between successive calls to
- // the member functions of this class
- // (i.e., we never care what values
- // they were set to the last time a
- // member function was called), we
- // mark these vectors as <code>mutable</code>.
- //
- // The rest of the (short)
- // implementation of this class is
- // straightforward if you know the
- // order of matrix-vector
- // multiplications performed by the
- // <code>vmult</code> function:
-class SchurComplement : public Subscriptor
-{
- public:
- SchurComplement (const BlockSparseMatrix<double> &A,
- const IterativeInverse<Vector<double> > &Minv);
-
- void vmult (Vector<double> &dst,
- const Vector<double> &src) const;
-
- private:
- const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
- const SmartPointer<const IterativeInverse<Vector<double> > > m_inverse;
-
- mutable Vector<double> tmp1, tmp2;
-};
-
-
-SchurComplement::SchurComplement (const BlockSparseMatrix<double> &A,
- const IterativeInverse<Vector<double> > &Minv)
- :
- system_matrix (&A),
- m_inverse (&Minv),
- tmp1 (A.block(0,0).m()),
- tmp2 (A.block(0,0).m())
-{}
-
-
-void SchurComplement::vmult (Vector<double> &dst,
- const Vector<double> &src) const
-{
- system_matrix->block(0,1).vmult (tmp1, src);
- m_inverse->vmult (tmp2, tmp1);
- system_matrix->block(1,0).vmult (dst, tmp2);
-}
-
-
// @sect4{The <code>ApproximateSchurComplement</code> class template}
// The third component of our solver
m_inverse.solver.select("cg");
static ReductionControl inner_control(1000, 0., 1.e-13);
m_inverse.solver.set_control(inner_control);
-
+
+ SchurComplement<IterativeInverse<Vector<double> >,
+ SparseMatrix<double>,
+ SparseMatrix<double>,
+ SparseMatrix<double> >
+ schur_complement (m_inverse, system_matrix.block(1,0),
+ system_matrix.block(1,0));
+
Vector<double> tmp (solution.block(0).size());
// Now on to the first
schur_rhs -= system_rhs.block(1);
- SchurComplement
- schur_complement (system_matrix, m_inverse);
ApproximateSchurComplement
approximate_schur_complement (system_matrix);
* call to LAPACKFullMatrix::compute_inverse_svd().
*/
double threshold;
+
+ /**
+ * The order in which blocks
+ * should be traversed. This
+ * vector can initiate
+ * several modes of
+ * execution:
+ *
+ * <ol>
+ * <li>If the length of the
+ * vector is zero, then the
+ * relaxation method will be
+ * exectued from first to
+ * last block.
+ * <li> If the length is one,
+ * then the inner vector must
+ * have the same size as
+ * the number of blocks. The
+ * relaxation method is
+ * applied in the order given
+ * in this vector.
+ * <li> If the outer vector
+ * has length greater one,
+ * then the relaxation method
+ * is applied several times,
+ * each time in the order
+ * given by the inner vector
+ * of the corresponding
+ * index. This mode can for
+ * instance be used for ADI
+ * methods and similar
+ * direction sweeps.
+ * </ol>
+ */
+ std::vector<std::vector<unsigned int> > order;
};
/**
const MATRIX &M=*this->A;
Vector<number2> b_cell, x_cell;
+ const bool permutation_empty = additional_data->order.size() == 0;
+ const unsigned int n_permutations = (permutation_empty)
+ ? 1U : additional_data->order.size();
const unsigned int n_blocks = additional_data->block_list.size();
- for (unsigned int bi=0;bi<n_blocks;++bi)
+
+ for (unsigned int perm=0; perm<n_permutations;++perm)
{
- const unsigned int block = backward ? (n_blocks - bi - 1) : bi;
- const unsigned int bs = additional_data->block_list.block_size(block);
-
- b_cell.reinit(bs);
- x_cell.reinit(bs);
- // Collect off-diagonal parts
- BlockList::const_iterator row = additional_data->block_list.begin(block);
- for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)
+ for (unsigned int bi=0;bi<n_blocks;++bi)
{
- b_cell(row_cell) = src(*row);
- for (typename MATRIX::const_iterator entry = M.begin(*row);
- entry != M.end(*row); ++entry)
- b_cell(row_cell) -= entry->value() * prev(entry->column());
+ unsigned int block = backward ? (n_blocks - bi - 1) : bi;
+ if (!permutation_empty)
+ block = additional_data->order[perm][block];
+
+ const unsigned int bs = additional_data->block_list.block_size(block);
+
+ b_cell.reinit(bs);
+ x_cell.reinit(bs);
+ // Collect off-diagonal parts
+ BlockList::const_iterator row = additional_data->block_list.begin(block);
+ for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)
+ {
+ b_cell(row_cell) = src(*row);
+ for (typename MATRIX::const_iterator entry = M.begin(*row);
+ entry != M.end(*row); ++entry)
+ b_cell(row_cell) -= entry->value() * prev(entry->column());
+ }
+ // Apply inverse diagonal
+ this->inverse_vmult(block, x_cell, b_cell);
+ // Store in result vector
+ row=additional_data->block_list.begin(block);
+ for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)
+ dst(*row) = prev(*row) + additional_data->relaxation * x_cell(row_cell);
}
- // Apply inverse diagonal
- this->inverse_vmult(block, x_cell, b_cell);
- // Store in result vector
- row=additional_data->block_list.begin(block);
- for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)
- dst(*row) = prev(*row) + additional_data->relaxation * x_cell(row_cell);
}
}