* builds the pattern on the fly. Filling the matrix may consume more
* time as with @p{SparseMatrix}, since large memory movements may be
* involved. To help optimizing things, an expected row-length may be
- * provided to the constructor, as well as a mininmum size increment
+ * provided to the constructor, as well as an increment size
* for rows.
*
+ * The storage structure: like with the usual sparse matrix, it is
+ * attempted to store only non-zero elements. these are stored in a
+ * single data array @p{data}. They are ordered by row and inside each
+ * row by column number. Each row is described by its starting point
+ * in the data array and its length. These are stored in the
+ * @p{row_info} array, together with additional useful information.
+ *
+ * Due to the structure, gaps may occur between rows. Whenever a new
+ * entry must be created, an attempt is made to use the gap in its
+ * row. If the gap is full, the row must be extended and all
+ * subsequent rows must be shifted backwards. This is a very expensive
+ * operation and should be avoided as much as possible.
+ *
+ * This is, where the optimization parameters, provided to the
+ * constructor or to the function @p{reinit} come
+ * in. @p{default_row_length} is the amount of entries that will be
+ * allocated for each row on initialization (the actual length of the
+ * rows is still zero). This means, that @p{default_row_length}
+ * entries can be added to this row without shifting other rows. If
+ * less entries are added, the additional memory will be wasted.
+ *
+ * If the space for a row is not sufficient, then it is enlarged by
+ * @p{default_increment} entries. This way, the subsequent rows are
+ * not shifted by single entries very often.
+ *
+ * Finally, the @p{default_reserve} allocates extra space at the end
+ * of the data array. This space is used whenever a row must be
+ * enlarged. Since @p{std::vector} doubles the capacity everytime it
+ * must increase it, this value should allow for all the growth needed.
+ *
+ * Suggested settings: @p{default_row_length} should be the length of
+ * a typical row, for instance the size of the stencil in regular
+ * parts of the grid. Then, @p{default_increment} may be the expected
+ * amount of entries added to the row by having one hanging node. This
+ * way, a good compromise between memory consumption and speed should
+ * be achieved. @p{default_reserve} should then be an estimate for the
+ * number of hanging nodes times @p{default_increment}.
+ *
+ * If the rows are expected to be filled more or less from first to
+ * last, using a @p{default_row_length} may not be such a bad idea.
+ *
* The name of this matrix is in reverence to a publication of the
* Internal Revenue Service of the United States of America. I hope
* some other aliens will appreciate it. By the way, the suffix makes
static const unsigned short
invalid_diagonal = static_cast<unsigned short>(-1);
};
-
-
+
public:
+ /**
+ * Accessor class for iterators
+ */
+ class Accessor
+ {
+ public:
+ /**
+ * Constructor. Since we use
+ * accessors only for read
+ * access, a const matrix
+ * pointer is sufficient.
+ */
+ Accessor (const SparseMatrixEZ<number>*,
+ unsigned int row,
+ unsigned short index);
+
+ /**
+ * Row number of the element
+ * represented by this
+ * object.
+ */
+ unsigned int row() const;
+
+ /**
+ * Index in row of the element
+ * represented by this
+ * object.
+ */
+ unsigned short index() const;
+
+ /**
+ * Column number of the
+ * element represented by
+ * this object.
+ */
+ unsigned int column() const;
+
+ /**
+ * Value of this matrix entry.
+ */
+ number value() const;
+
+ protected:
+ /**
+ * The matrix accessed.
+ */
+ const SparseMatrixEZ<number>* matrix;
+
+ /**
+ * Current row number.
+ */
+ unsigned int a_row;
+
+ /**
+ * Current index in row.
+ */
+ unsigned short a_index;
+ };
+
+ /**
+ * STL conforming iterator.
+ */
+ class const_iterator : private Accessor
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ const_iterator(const SparseMatrixEZ<number>*,
+ unsigned int row,
+ unsigned short index);
+
+ /**
+ * Prefix increment.
+ */
+ const_iterator& operator++ ();
+
+ /**
+ * Postfix increment.
+ */
+ const_iterator& operator++ (int);
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor& operator* () const;
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor* operator-> () const;
+
+ /**
+ * Comparison. True, if
+ * both iterators point to
+ * the same matrix
+ * position.
+ */
+ bool operator == (const const_iterator&) const;
+ /**
+ * Inverse of @p{==}.
+ */
+ bool operator != (const const_iterator&) const;
+
+ /**
+ * Comparison
+ * operator. Result is true
+ * if either the first row
+ * number is smaller or if
+ * the row numbers are
+ * equal and the first
+ * index is smaller.
+ */
+ bool operator < (const const_iterator&) const;
+ };
+
/**
* Type of matrix entries. In analogy to
* the STL container classes.
*/
unsigned int n () const;
- /**
- * Return the number of nonzero
- * elements of this
- * matrix. Actually, it returns
- * the number of entries in the
- * sparsity pattern; if any of
- * the entries should happen to
- * be zero, it is counted anyway.
- */
-// unsigned int n_nonzero_elements () const;
-
- /**
- * Return the number of actually
- * nonzero elements of this
- * matrix.
- *
- * Note, that this function does
- * (in contrary to the
- * @p{n_nonzero_elements}) NOT
- * count all entries of the
- * sparsity pattern but only the
- * ones that are nonzero.
- */
-// unsigned int n_actually_nonzero_elements () const;
-
/**
* Set the element @p{(i,j)} to
* @p{value}. Allocates the entry
number operator () (const unsigned int i,
const unsigned int j) const;
+ /**
+ * Return the value of the entry
+ * (i,j). Returns zero for all
+ * non-existing entries.
+ */
+ number el (const unsigned int i,
+ const unsigned int j) const;
+
/**
* Return the main diagonal element in
* the @p{i}th row. This function throws an
* involve searching for the
* right column number.
*/
- number diag_element (const unsigned int i) const;
+// number diag_element (const unsigned int i) const;
/**
* Same as above, but return a
* writeable reference. You're
* sure you know what you do?
*/
- number & diag_element (const unsigned int i);
+// number & diag_element (const unsigned int i);
/**
* Matrix-vector multiplication:
void precondition_TSOR (Vector<somenumber> &dst,
const Vector<somenumber> &src,
const number om = 1.) const;
+
+ /**
+ * STL-like iterator with the
+ * first entry.
+ */
+ const_iterator begin () const;
+
+ /**
+ * Final iterator.
+ */
+ const_iterator end () const;
+
+ /**
+ * STL-like iterator with the
+ * first entry of row @p{r}.
+ */
+ const_iterator begin (unsigned int r) const;
+
+ /**
+ * Final iterator of row @p{r}.
+ */
+ const_iterator end (unsigned int r) const;
+
+ /**
+ * Return the number of nonzero
+ * elements of this
+ * matrix.
+ */
+// unsigned int n_nonzero_elements () const;
+
+ /**
+ * Return the number of actually
+ * nonzero elements of this
+ * matrix.
+ */
+// unsigned int n_actually_nonzero_elements () const;
/**
* Print the matrix to the given
*/
DeclException0(ExcNoDiagonal);
- private:
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidIndex,
+ int, int,
+ << "The entry with index <" << arg1 << ',' << arg2
+ << "> does not exist.");
+ private:
/**
* Find an entry. Return a
* zero-pointer if the entry does
* Increment when a row grows.
*/
unsigned int increment;
-
- // make all other sparse matrices
- // friends
- template <typename somenumber> friend class SparseMatrix;
+
+ friend class Accessor;
+ friend class const_iterator;
};
{}
-// template <typename number>
-// inline
-// bool
-// SparseMatrixEZ<number>::Entry::operator==(const Entry& e) const
-// {
-// return column == e.column;
-// }
+template <typename number>
+inline
+SparseMatrixEZ<number>::RowInfo::RowInfo(unsigned int start)
+ :
+ start(start), length(0), diagonal(invalid_diagonal)
+{}
+
+
+//----------------------------------------------------------------------//
+template <typename number>
+inline
+SparseMatrixEZ<number>::Accessor::Accessor (
+ const SparseMatrixEZ<number>* matrix,
+ unsigned int r,
+ unsigned short i)
+ :
+ matrix(matrix),
+ a_row(r),
+ a_index(i)
+{}
+
+
+template <typename number>
+inline
+unsigned int
+SparseMatrixEZ<number>::Accessor::row() const
+{
+ return a_row;
+}
+
+
+template <typename number>
+inline
+unsigned int
+SparseMatrixEZ<number>::Accessor::column() const
+{
+ return matrix->data[matrix->row_info[a_row].start+a_index].column;
+}
+
+
+template <typename number>
+inline
+unsigned short
+SparseMatrixEZ<number>::Accessor::index() const
+{
+ return a_index;
+}
-// template <typename number>
-// inline
-// bool
-// SparseMatrixEZ<number>::Entry::operator<(const Entry& e) const
-// {
-// return column < e.column;
-// }
template <typename number>
inline
-SparseMatrixEZ<number>::RowInfo::RowInfo(unsigned int start)
+number
+SparseMatrixEZ<number>::Accessor::value() const
+{
+ return matrix->data[matrix->row_info[a_row].start+a_index].value;
+}
+
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator::const_iterator(
+ const SparseMatrixEZ<number>* matrix,
+ unsigned int r,
+ unsigned short i)
:
- start(start), length(0), diagonal(invalid_diagonal)
+ Accessor(matrix, r, i)
{}
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator&
+SparseMatrixEZ<number>::const_iterator::operator++ ()
+{
+ Assert (a_row < matrix->m(), ExcIteratorPastEnd());
+
+ ++a_index;
+ if (a_index >= matrix->row_info[a_row].length)
+ {
+ a_index = 0;
+ a_row++;
+ }
+ return *this;
+}
+
+
+template <typename number>
+inline
+const SparseMatrixEZ<number>::Accessor&
+SparseMatrixEZ<number>::const_iterator::operator* () const
+{
+ return *this;
+}
+
+
+template <typename number>
+inline
+const SparseMatrixEZ<number>::Accessor*
+SparseMatrixEZ<number>::const_iterator::operator-> () const
+{
+ return this;
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrixEZ<number>::const_iterator::operator == (
+ const const_iterator& other) const
+{
+ return (row() == other->row() && index() == other->index());
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrixEZ<number>::const_iterator::operator != (
+ const const_iterator& other) const
+{
+ return ! (*this == other);
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrixEZ<number>::const_iterator::operator < (
+ const const_iterator& other) const
+{
+ return (row() < other->row() ||
+ (row() == other->row() && index() < other->index()));
+}
+
+
//----------------------------------------------------------------------//
template <typename number>
inline
};
+
+template <typename number>
+inline
+number SparseMatrixEZ<number>::el (const unsigned int i,
+ const unsigned int j) const
+{
+ Entry* entry = locate(i,j);
+ if (entry)
+ return entry->value;
+ return 0.;
+}
+
+
+
+template <typename number>
+inline
+number SparseMatrixEZ<number>::operator() (const unsigned int i,
+ const unsigned int j) const
+{
+ Entry* entry = locate(i,j);
+ if (entry)
+ return entry->value;
+ Assert(false, ExcInvalidEntry(i,j));
+ return 0.;
+}
+
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::begin () const
+{
+ return const_iterator(this, 0, 0);
+}
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::end () const
+{
+ return const_iterator(this, m(), 0);
+}
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::begin (unsigned int r) const
+{
+ Assert (r<m(), ExcIndexRange(r,0,m()));
+ return const_iterator(this, r, 0);
+}
+
+template <typename number>
+inline
+SparseMatrixEZ<number>::const_iterator
+SparseMatrixEZ<number>::end (unsigned int r) const
+{
+ Assert (r<m(), ExcIndexRange(r,0,m()));
+ return const_iterator(this, r+1, 0);
+}
+
#endif
/*---------------------------- sparse_matrix.h ---------------------------*/
#include <fstream>
+#define PREC_CHECK(solver, method, precond) try \
+ { solver.method (A, u, f, precond); } catch (...) {} \
+ residuals.push_back(control.last_value())
+
template<class MATRIX>
void
-check_vmult_quadratic(const MATRIX& A,
+check_vmult_quadratic(std::vector<double>& residuals,
+ const MATRIX& A,
const char* prefix)
{
deallog.push(prefix);
SolverControl control(10, 1.e-13, false);
SolverRichardson<> rich(control, mem, .01);
SolverRichardson<> prich(control, mem, 1.);
- PreconditionIdentity prec;
+ PreconditionIdentity identity;
PreconditionJacobi<MATRIX> jacobi;
jacobi.initialize(A, .5);
PreconditionSOR<MATRIX> sor;
u = 0.;
f = 1.;
- try { rich.solve(A, u, f, prec); } catch (...) {}
- try { prich.solve(A, u, f, jacobi); } catch (...) {}
- try { prich.solve(A, u, f, ssor); } catch (...) {}
- try { prich.solve(A, u, f, sor); } catch (...) {}
-
+ PREC_CHECK(rich, solve, identity);
+ PREC_CHECK(prich, solve, jacobi);
+ PREC_CHECK(prich, solve, ssor);
+ PREC_CHECK(prich, solve, sor);
+
u = 0.;
deallog << "Transpose" << std::endl;
- try { rich.Tsolve(A, u, f, prec); } catch (...) {}
- try { prich.Tsolve(A, u, f, jacobi); } catch (...) {}
- try { prich.Tsolve(A, u, f, ssor); } catch (...) {}
- try { prich.Tsolve(A, u, f, sor); } catch (...) {}
+ PREC_CHECK(rich, Tsolve, identity);
+ PREC_CHECK(prich, Tsolve, jacobi);
+ PREC_CHECK(prich, Tsolve, ssor);
+ PREC_CHECK(prich, Tsolve, sor);
deallog.pop();
}
+
+template <class MATRIX>
+void
+check_iterator (const MATRIX& A)
+{
+ for (typename MATRIX::const_iterator i = A.begin(); i!= A.end(); ++i)
+ deallog << '\t' << i->row()
+ << '\t' << i->column()
+ << '\t' << i->index()
+ << '\t' << i->value()
+ << std::endl;
+ deallog << "Repeat row 2" << std::endl;
+ for (typename MATRIX::const_iterator i = A.begin(2); i!= A.end(2); ++i)
+ deallog << '\t' << i->row()
+ << '\t' << i->column()
+ << '\t' << i->index()
+ << '\t' << i->value()
+ << std::endl;
+}
+
+
int main()
{
std::ofstream logfile("sparse_matrices.output");
FDMatrix testproblem (size, size);
unsigned int dim = (size-1)*(size-1);
+ std::vector<double> A_res;
+ std::vector<double> E_res;
+
deallog << "Structure" << std::endl;
SparsityPattern structure(dim, dim, 5);
testproblem.five_point_structure(structure);
SparseMatrix<double> A(structure);
deallog << "Assemble" << std::endl;
testproblem.five_point(A, true);
- check_vmult_quadratic(A, "5-SparseMatrix<double>");
+ check_vmult_quadratic(A_res, A, "5-SparseMatrix<double>");
SparseMatrixEZ<double> E(dim,dim,row_length,2);
deallog << "Assemble" << std::endl;
testproblem.five_point(E, true);
- check_vmult_quadratic(E, "5-SparseMatrixEZ<double>");
-
+ check_vmult_quadratic(E_res, E, "5-SparseMatrixEZ<double>");
+#ifdef DEBUG
+ check_iterator(E);
+#endif
+
A.clear();
deallog << "Structure" << std::endl;
structure.reinit(dim, dim, 9);
A.reinit(structure);
deallog << "Assemble" << std::endl;
testproblem.nine_point(A);
- check_vmult_quadratic(A, "9-SparseMatrix<double>");
+ check_vmult_quadratic(A_res, A, "9-SparseMatrix<double>");
E.clear();
E.reinit(dim,dim,row_length,2);
deallog << "Assemble" << std::endl;
testproblem.nine_point(E);
- check_vmult_quadratic(E, "9-SparseMatrixEZ<double>");
-
+ check_vmult_quadratic(E_res, E, "9-SparseMatrixEZ<double>");
+
+ for (unsigned int i=0;i<A_res.size();++i)
+ if (std::fabs(A_res[i] - E_res[i]) > 1.e-14)
+ deallog << "SparseMatrix and SparseMatrixEZ differ!!!"
+ << std::endl;
}