#include <cstdlib>
#include <cstdio>
#include <iomanip>
+#include <algorithm>
template <typename number>
FullMatrix<number>::FullMatrix (const FullMatrix &m)
{
init (m.dim_image, m.dim_range);
- number * p = &val[0];
- const number * vp = &m.val[0];
- const number * const e = &val[dim_image*dim_range];
-
- while (p!=e)
- *p++ = *vp++;
+ if (dim_range*dim_image != 0)
+ copy (&m.val[0], &m.val[dim_image*dim_range],
+ &val[0]);
};
template <typename number>
template <typename number2>
-void FullMatrix<number>::Tvmult (Vector<number2>& dst, const Vector<number2>& src, const bool adding) const
+void FullMatrix<number>::Tvmult (Vector<number2>& dst,
+ const Vector<number2>& src,
+ const bool adding) const
{
Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m()));
template <typename number>
template <typename number2, typename number3>
-double FullMatrix<number>::residual (Vector<number2>& dst, const Vector<number2>& src,
- const Vector<number3>& right) const
+double FullMatrix<number>::residual (Vector<number2>& dst,
+ const Vector<number2>& src,
+ const Vector<number3>& right) const
{
Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
FullMatrix<number>::operator = (const FullMatrix<number>& m)
{
reinit(m);
-
- number * p = &val[0];
- const number * vp = &m.val[0];
- const number * const e = &val[dim_image*dim_range];
-
- while (p!=e)
- *p++ = *vp++;
-
+ if (dim_range*dim_image != 0)
+ copy (&m.val[0], &m.val[dim_image*dim_range],
+ &val[0]);
+
return *this;
}
FullMatrix<number>::operator = (const FullMatrix<number2>& m)
{
reinit(m);
-
- number * p = &val[0];
- const number2 * vp = &m.val[0];
- const number * const e = &val[dim_image*dim_range];
-
- while (p!=e)
- *p++ = *vp++;
+ if (dim_range*dim_image != 0)
+ copy (&m.val[0], &m.val[dim_image*dim_range],
+ &val[0]);
return *this;
}
template <typename number>
template <typename number2>
void FullMatrix<number>::fill (const FullMatrix<number2>& src,
- const unsigned int i, const unsigned int j)
+ const unsigned int i,
+ const unsigned int j)
{
Assert (n() >= src.n() + j, ExcInvalidDestination(n(), src.n(), j));
Assert (m() >= src.m() + i, ExcInvalidDestination(m(), src.m(), i));
template <typename number>
void FullMatrix<number>::add_row (const unsigned int i,
- const number s, const unsigned int j)
+ const number s,
+ const unsigned int j)
{
for (unsigned int k=0; k<m(); ++k)
el(i,k) += s*el(j,k);
template <typename number>
-void FullMatrix<number>::add_row (const unsigned int i, const number s,
- const unsigned int j, const number t,
- const unsigned int k)
+void FullMatrix<number>::add_row (const unsigned int i,
+ const number s,
+ const unsigned int j,
+ const number t,
+ const unsigned int k)
{
const unsigned int size_m = m();
for (unsigned l=0; l<size_m; ++l)
template <typename number>
template <typename number2>
-void FullMatrix<number>::mmult (FullMatrix<number2>& dst, const FullMatrix<number2>& src) const
+void FullMatrix<number>::mmult (FullMatrix<number2>& dst,
+ const FullMatrix<number2>& src) const
{
Assert (n() == src.m(), ExcDimensionMismatch(n(), src.m()));
unsigned int i,j,k;
bool q = (dim_range==m.dim_range) && (dim_image==m.dim_image);
if (!q) return false;
- for (unsigned int i=0; i<dim_image; ++i)
- for (unsigned int j=0; j<dim_range; ++j)
- if (el(i,j) != m.el(i,j)) return false;
- return true;
+ return equal (&val[0], &val[dim_range*dim_image],
+ &m.val[0]);
};
template <typename number>
void FullMatrix<number>::clear ()
{
- number *val_ptr = &val[0];
- const number *end_ptr = &val[n()*m()];
- while (val_ptr != end_ptr)
- *val_ptr++ = 0.;
+ fill_n (&val[0], n()*m(), 0);
};
* Make a copy with extra off-diagonals.
*
* This constructs objects intended for
- * the application of the ILU(n)-method.
+ * the application of the ILU(n)-method
+ * or other incomplete decompositions.
* Therefore, additional to the original
- * entry structure, space for #extra_cols#
- * side-diagonals is provided.
+ * entry structure, space for
+ * #extra_off_diagonals#
+ * side-diagonals is provided on both
+ * sides of the main diagonal.
+ *
+ * #max_per_row# is the maximum number of
+ * nonzero elements per row which this
+ * structure is to hold. It is assumed
+ * that this number is sufficiently large
+ * to accomodate both the elements in
+ * #original# as well as the new
+ * off-diagonal elements created by this
+ * constructor. You will usually want to
+ * give the same number as you gave for
+ * #original# plus the number of side
+ * diagonals times two. You may however
+ * give a larger value if you wish to add
+ * further nonzero entries for the
+ * decomposition based on other criteria
+ * than their being on side-diagonals.
+ *
+ * This function requires that #original#
+ * refer to a square matrix structure.
+ * It shall be compressed. The matrix
+ * structure is not compressed
+ * after this function finishes.
*/
SparseMatrixStruct(const SparseMatrixStruct& original,
- unsigned int extra_cols);
+ const unsigned int max_per_row,
+ const unsigned int extra_off_diagonals);
/**
* Destructor.
* Exception
*/
DeclException0 (ExcInvalidConstructorCall);
-
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNotSquare);
+
private:
unsigned int max_dim;
unsigned int rows, cols;
/**
* Constructor. Takes the given matrix
- * sparisty structure to represent the
+ * sparsity structure to represent the
* sparsity pattern of this matrix. You
* can change the sparsity pattern later
* on by calling the #reinit# function.
* To remember: the matrix is of dimension
* $m \times n$.
*/
+
unsigned int n () const;
/**
*/
number diag_element (const unsigned int i) const;
+ /**
+ * Same as above, but return a writeable
+ * reference. You're sure what you do?
+ */
+ number & diag_element (const unsigned int i);
+
/**
* This is kind of an expert mode: get
* access to the #i#th element of this
+template <typename number>
+inline
+number & SparseMatrix<number>::diag_element (const unsigned int i) {
+ Assert (cols != 0, ExcMatrixNotInitialized());
+ Assert (m() == n(), ExcMatrixNotSquare());
+ Assert (i<max_len, ExcInvalidIndex1(i));
+
+ // Use that the first element in each
+ // row of a square matrix is the main
+ // diagonal
+ return val[cols->rowstart[i]];
+};
+
+
+
template <typename number>
inline
number SparseMatrix<number>::global_entry (const unsigned int j) const {
Assert (val != 0, ExcMatrixNotInitialized());
Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
- number *val_ptr = &val[0];
- const somenumber *matrix_ptr = &matrix.val[0];
- const number *const end_ptr = &val[cols->vec_len];
-
- while (val_ptr != end_ptr)
- *val_ptr++ = *matrix_ptr++;
+ copy (&matrix.val[0], &matrix.val[cols->vec_len],
+ &val[0]);
is_ilu = false;
return *this;
{
if (v.dim != dim)
reinit (v.dim, true);
-
if (dim!=0)
copy (v.begin(), v.end(), begin());
+
return *this;
}
#include <iomanip>
#include <algorithm>
#include <cmath>
-
+#include <numeric>
SparseMatrixStruct::SparseMatrixStruct () :
-SparseMatrixStruct::SparseMatrixStruct (const unsigned int m, const unsigned int n,
+SparseMatrixStruct::SparseMatrixStruct (const unsigned int m,
+ const unsigned int n,
const unsigned int max_per_row)
: max_dim(0),
max_vec_len(0),
};
+SparseMatrixStruct::SparseMatrixStruct (const SparseMatrixStruct &original,
+ const unsigned int max_per_row,
+ const unsigned int extra_off_diagonals)
+ : max_dim(0),
+ max_vec_len(0),
+ rowstart(0),
+ colnums(0)
+{
+ Assert (original.rows==original.cols, ExcNotSquare());
+ Assert (original.is_compressed(), ExcNotCompressed());
+
+ reinit (original.rows, original.cols, max_per_row);
+
+ // now copy the entries from
+ // the other object
+ for (unsigned int row=0; row<original.rows; ++row)
+ {
+ // copy the elements of this row
+ // of the other object
+ //
+ // note that the first object actually
+ // is the main-diagonal element,
+ // which we need not copy
+ //
+ // we do the copying in two steps:
+ // first we note that the elements in
+ // #original# are sorted, so we may
+ // first copy all the elements up to
+ // the first side-diagonal one which
+ // is to be filled in. then we insert
+ // the side-diagonals, finally copy
+ // the rest from that element onwards
+ // which is not a side-diagonal any
+ // more.
+ const int * const
+ original_row_start = &original.colnums[original.rowstart[row]] + 1;
+ // thw following requires that
+ // #original# be compressed since
+ // otherwise there might be -1's
+ const int * const
+ original_row_end = &original.colnums[original.rowstart[row+1]];
+
+ const int * const
+ original_last_before_side_diagonals = lower_bound (original_row_start,
+ original_row_end,
+ static_cast<int>
+ (row
+ -extra_off_diagonals));
+ const int * const
+ original_first_after_side_diagonals = upper_bound (original_row_start,
+ original_row_end,
+ static_cast<int>
+ (row
+ +extra_off_diagonals));
+
+ int * next_free_slot = &colnums[rowstart[row]] + 1;
+
+ // copy elements before side-diagonals
+ next_free_slot = copy (original_row_start,
+ original_last_before_side_diagonals,
+ next_free_slot);
+
+ // insert left and right side-diagonals
+ for (unsigned int i=1; i<=min(row,extra_off_diagonals);
+ ++i, ++next_free_slot)
+ *next_free_slot = row-i;
+ for (unsigned int i=1; i<=min(extra_off_diagonals, rows-row-1);
+ ++i, ++next_free_slot)
+ *next_free_slot = row+i;
+
+ // copy rest
+ next_free_slot = copy (original_first_after_side_diagonals,
+ original_row_end,
+ next_free_slot);
+
+ // this error may happen if the
+ // sum of previous elements per row
+ // and those of the new diagonals
+ // exceeds the maximum number of
+ // elements per row given to this
+ // constructor
+ Assert (next_free_slot <= &colnums[rowstart[row+1]],
+ ExcNotEnoughSpace (0,rowstart[row+1]-rowstart[row]));
+ };
+};
+
+
SparseMatrixStruct::~SparseMatrixStruct ()
{
// all other entries are sorted, so
// we can use a binary seach algorithm
- const int* p = lower_bound (&colnums[rowstart[i]+1],
- &colnums[rowstart[i+1]],
- static_cast<signed int>(j));
+ const int * const p = lower_bound (&colnums[rowstart[i]+1],
+ &colnums[rowstart[i+1]],
+ static_cast<signed int>(j));
if ((*p == static_cast<signed int>(j)) &&
(p != &colnums[rowstart[i+1]]))
return (p - &colnums[0]);