From: wolf Date: Sun, 25 Apr 1999 22:15:29 +0000 (+0000) Subject: Use faster copy operations for vectors and matrices. Some other small updates also. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4b3f7c50f750333efdeeb66a6a6fae16ec2effcf;p=dealii-svn.git Use faster copy operations for vectors and matrices. Some other small updates also. git-svn-id: https://svn.dealii.org/trunk@1205 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index 1289e6743e..fced113f2f 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -8,6 +8,7 @@ #include #include #include +#include template @@ -28,12 +29,9 @@ template FullMatrix::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]); }; @@ -373,7 +371,9 @@ void FullMatrix::gsmult (Vector& dst, const Vector& sr template template -void FullMatrix::Tvmult (Vector& dst, const Vector& src, const bool adding) const +void FullMatrix::Tvmult (Vector& dst, + const Vector& src, + const bool adding) const { Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m())); @@ -396,8 +396,9 @@ void FullMatrix::Tvmult (Vector& dst, const Vector& sr template template -double FullMatrix::residual (Vector& dst, const Vector& src, - const Vector& right) const +double FullMatrix::residual (Vector& dst, + const Vector& src, + const Vector& right) const { Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); @@ -464,14 +465,10 @@ FullMatrix& FullMatrix::operator = (const FullMatrix& 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; } @@ -483,13 +480,9 @@ FullMatrix& FullMatrix::operator = (const FullMatrix& 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; } @@ -499,7 +492,8 @@ FullMatrix::operator = (const FullMatrix& m) template template void FullMatrix::fill (const FullMatrix& 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)); @@ -513,7 +507,8 @@ void FullMatrix::fill (const FullMatrix& src, template void FullMatrix::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::add_row (const unsigned int i, template -void FullMatrix::add_row (const unsigned int i, const number s, - const unsigned int j, const number t, - const unsigned int k) +void FullMatrix::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::diagadd (const number src) template template -void FullMatrix::mmult (FullMatrix& dst, const FullMatrix& src) const +void FullMatrix::mmult (FullMatrix& dst, + const FullMatrix& src) const { Assert (n() == src.m(), ExcDimensionMismatch(n(), src.m())); unsigned int i,j,k; @@ -1151,10 +1149,8 @@ FullMatrix::operator == (const FullMatrix &m) const bool q = (dim_range==m.dim_range) && (dim_image==m.dim_image); if (!q) return false; - for (unsigned int i=0; i::norm2 () const template void FullMatrix::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); }; diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 00a4b72487..020600c685 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -89,13 +89,39 @@ class SparseMatrixStruct * 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. @@ -364,7 +390,11 @@ class SparseMatrixStruct * Exception */ DeclException0 (ExcInvalidConstructorCall); - + /** + * Exception + */ + DeclException0 (ExcNotSquare); + private: unsigned int max_dim; unsigned int rows, cols; @@ -430,7 +460,7 @@ class SparseMatrix /** * 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. @@ -514,6 +544,7 @@ class SparseMatrix * To remember: the matrix is of dimension * $m \times n$. */ + unsigned int n () const; /** @@ -641,6 +672,12 @@ class SparseMatrix */ 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 @@ -1003,6 +1040,21 @@ number SparseMatrix::diag_element (const unsigned int i) const { +template +inline +number & SparseMatrix::diag_element (const unsigned int i) { + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (irowstart[i]]; +}; + + + template inline number SparseMatrix::global_entry (const unsigned int j) const { diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 4688022c34..54565fb685 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -143,12 +143,8 @@ SparseMatrix::copy_from (const SparseMatrix &matrix) 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; diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index d5d6f360f4..f579a75c57 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -530,9 +530,9 @@ Vector& Vector::operator = (const Vector& v) { if (v.dim != dim) reinit (v.dim, true); - if (dim!=0) copy (v.begin(), v.end(), begin()); + return *this; } diff --git a/deal.II/lac/source/sparse_matrix.cc b/deal.II/lac/source/sparse_matrix.cc index 5e16cd89c9..e52db7de6f 100644 --- a/deal.II/lac/source/sparse_matrix.cc +++ b/deal.II/lac/source/sparse_matrix.cc @@ -14,7 +14,7 @@ #include #include #include - +#include SparseMatrixStruct::SparseMatrixStruct () : @@ -44,7 +44,8 @@ SparseMatrixStruct::SparseMatrixStruct (const SparseMatrixStruct &s) : -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), @@ -67,6 +68,93 @@ SparseMatrixStruct::SparseMatrixStruct (const unsigned int n, }; +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 + (row + -extra_off_diagonals)); + const int * const + original_first_after_side_diagonals = upper_bound (original_row_start, + original_row_end, + static_cast + (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 () { @@ -256,9 +344,9 @@ SparseMatrixStruct::operator () (const unsigned int i, const unsigned int j) con // 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(j)); + const int * const p = lower_bound (&colnums[rowstart[i]+1], + &colnums[rowstart[i+1]], + static_cast(j)); if ((*p == static_cast(j)) && (p != &colnums[rowstart[i+1]])) return (p - &colnums[0]);