* whether this might pose some problems in the inversion of the local matrices.
* Maybe someone would like to check this.
*
- * @author Guido Kanschat, documentation mostly by Wolfgang Bangerth; 1999
+ * @author Guido Kanschat, documentation and extensions by Wolfgang Bangerth; 1999, 2000
*/
template<typename number>
class SparseVanka
* longer than the Vanka
* object. The same is true for
* the matrix.
+ *
+ * The matrix #M# which is passed
+ * here may or may not be the
+ * same matrix for which this
+ * object shall act as
+ * preconditioner. In particular,
+ * it is conceivable that the
+ * preconditioner is build up for
+ * one matrix once, but is used
+ * for subsequent steps in a
+ * nonlinear process as well,
+ * where the matrix changes in
+ * each step slightly.
+ *
+ * If #conserve_mem# is #false#,
+ * then the inverses of the local
+ * systems are computed now, if
+ * the flag is #true#, then they
+ * are computed every time the
+ * preconditioner is
+ * applied. This saves some
+ * memory, but makes
+ * preconditioning very
+ * slow. Note also, that if the
+ * flag is #false#, the the
+ * contents of the matrix #M# at
+ * the time of calling this
+ * constructor are used, while if
+ * the flag is #true#, then the
+ * values in #M# at the time of
+ * preconditioning are used. This
+ * may lead to different results,
+ * obviously, of #M# changes.
*/
SparseVanka(const SparseMatrix<number> &M,
- const vector<bool> &selected);
+ const vector<bool> &selected,
+ const bool conserve_memory = false);
/**
* Destructor.
* Delete all allocated matrices.
*/
~SparseVanka();
-
+
/**
* Do the preconditioning.
* This function takes the residual
void operator() (Vector<number2> &dst,
const Vector<number2> &src) const;
- /**
- * Minimize memory consumption.
- * Activating this option reduces
- * memory needs of the Vanka object
- * to nearly zero. You pay for this
- * by a high increase of computing
- * time, since all local matrices
- * are built up and inverted every
- * time the Vanka operator is applied.
- */
- void conserve_memory();
-
/**
* Exception
*/
* that are tagged in #selected#.
*/
mutable vector<SmartPointer<FullMatrix<float> > > inverses;
+
+ /**
+ * Compute the inverses of all
+ * selected diagonal elements.
+ */
+ void compute_inverses ();
};
template<typename number>
SparseVanka<number>::SparseVanka(const SparseMatrix<number> &M,
- const vector<bool> &selected)
+ const vector<bool> &selected,
+ const bool conserve_mem)
:
matrix(&M),
selected(selected),
- conserve_mem(false),
+ conserve_mem(conserve_mem),
inverses(M.m(),0)
{
Assert (M.m() == M.n(),
ExcMatrixNotSquare ());
+
+ if (conserve_mem == false)
+ compute_inverses ();
}
SparseVanka<number>::~SparseVanka()
{
vector<SmartPointer<FullMatrix<float> > >::iterator i;
- for(i=inverses.begin();i!=inverses.end();++i)
+ for(i=inverses.begin(); i!=inverses.end(); ++i)
{
- FullMatrix<float>* p = *i;
+ FullMatrix<float> *p = *i;
*i = 0;
if (p != 0) delete p;
}
+template <typename number>
+void
+SparseVanka<number>::compute_inverses ()
+{
+ // first define an alias to the sparsity
+ // pattern of the matrix, since this
+ // will be used quite often
+ const SparseMatrixStruct &structure
+ = matrix->get_sparsity_pattern();
+
+ // traverse all rows of the matrix
+ // which are selected
+ for (unsigned int row=0; row< matrix->m() ; ++row)
+ if (selected[row] == true)
+ {
+ const unsigned int row_length = structure.row_length(row);
+ inverses[row] = new FullMatrix<float> (row_length,
+ row_length);
+ // mapping between:
+ // 1 column number of all
+ // entries in this row, and
+ // 2 the position within this
+ // row (as stored in the
+ // sparsematrixstruct object
+ //
+ // since we do not explicitely
+ // consider nonsysmmetric sparsity
+ // patterns, the first element
+ // of each entry simply denotes
+ // all degrees of freedom that
+ // couple with #row#.
+ map<unsigned int, unsigned int> local_index;
+ for (unsigned int i=0; i<row_length; ++i)
+ local_index.insert(pair<unsigned int, unsigned int>
+ (structure.column_number(row, i), i));
+
+ // Build local matrix and rhs
+ for (map<unsigned int, unsigned int>::const_iterator is=local_index.begin();
+ is!=local_index.end(); ++is)
+ {
+ // irow loops over all DoFs that
+ // couple with the present DoF
+ const unsigned int irow = is->first;
+ // index of DoF irow in the matrix
+ // row corresponding to DoF #row#.
+ // runs between 0 and row_length
+ const unsigned int i = is->second;
+ // number of DoFs coupling to
+ // irow (including irow itself)
+ const unsigned int irow_length = structure.row_length(irow);
+
+ // copy rhs
+// b(i) = src(irow);
+
+ // for all the DoFs that irow
+ // couples with
+ for (unsigned int j=0; j<irow_length; ++j)
+ {
+ // col is the number of
+ // this dof
+ const unsigned int col = structure.column_number(irow, j);
+ // find out whether this DoF
+ // (that couples with #irow#,
+ // which itself couples with
+ // #row#) also couples with
+ // #row#.
+ const map<unsigned int, unsigned int>::const_iterator js
+ = local_index.find(col);
+ // if not, then still use
+ // this dof to modify the rhs
+ //
+ // note that if so, we already
+ // have copied the entry above
+ if (js == local_index.end())
+;//b(i) -= matrix->raw_entry(irow,j) * dst(col);
+ else
+ // if so, then build the
+ // matrix out of it
+ (*inverses[row])(i,js->second) = matrix->raw_entry(irow,j);
+ };
+ };
+
+ // Compute new values
+ inverses[row]->gauss_jordan();
+ };
+};
+
+
+
template<typename number>
template<typename number2>
void
for (unsigned int row=0; row< matrix->m() ; ++row)
if (selected[row] == true)
{
- // shall we reconstruct the system
- // of equations for this DoF?
- bool build_matrix = true;
const unsigned int row_length = structure.row_length(row);
- if (!conserve_mem)
- {
- if (inverses[row] == 0)
- // inverse not yet built
- {
- FullMatrix<float> *p = new FullMatrix<float>;
- inverses[row] = p;
- } else {
- // inverse already built
- build_matrix = false;
- }
- }
-
// if we don't store the
// inverse matrices, then alias
// the entry in the global
{
inverses[row] = &local_matrix;
inverses[row]->reinit (row_length, row_length);
- }
- else
- if (build_matrix)
- inverses[row]->reinit (row_length, row_length);
+ };
b.reinit (row_length);
x.reinit (row_length);
else
// if so, then build the
// matrix out of it
- if (build_matrix)
+ if (conserve_mem == true)
(*inverses[row])(i,js->second) = matrix->raw_entry(irow,j);
};
};
// Compute new values
- if (build_matrix)
+ if (conserve_mem == true)
inverses[row]->gauss_jordan();
+
+ // apply preconditioner
inverses[row]->vmult(x,b);
// Distribute new values