]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Use faster copy operations for vectors and matrices. Some other small updates also.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 25 Apr 1999 22:15:29 +0000 (22:15 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 25 Apr 1999 22:15:29 +0000 (22:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@1205 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/include/lac/vector.templates.h
deal.II/lac/source/sparse_matrix.cc

index 1289e6743ef293df62aed0a365f0275d3cda0e18..fced113f2fcc6aa3cf509e7186001277f2f048e9 100644 (file)
@@ -8,6 +8,7 @@
 #include <cstdlib>
 #include <cstdio>
 #include <iomanip>
+#include <algorithm>
 
 
 template <typename number>
@@ -28,12 +29,9 @@ 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]);
 };
 
 
@@ -373,7 +371,9 @@ void FullMatrix<number>::gsmult (Vector<number2>& dst, const Vector<number2>& sr
 
 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()));
@@ -396,8 +396,9 @@ void FullMatrix<number>::Tvmult (Vector<number2>& dst, const Vector<number2>& sr
 
 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()));
@@ -464,14 +465,10 @@ FullMatrix<number>&
 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;
 }
 
@@ -483,13 +480,9 @@ FullMatrix<number>&
 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;
 }
@@ -499,7 +492,8 @@ FullMatrix<number>::operator = (const FullMatrix<number2>& m)
 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));
@@ -513,7 +507,8 @@ void FullMatrix<number>::fill (const FullMatrix<number2>& src,
 
 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);
@@ -522,9 +517,11 @@ void FullMatrix<number>::add_row (const unsigned int i,
 
 
 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)
@@ -590,7 +587,8 @@ void FullMatrix<number>::diagadd (const number src)
 
 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;
@@ -1151,10 +1149,8 @@ FullMatrix<number>::operator == (const FullMatrix<number> &m) const
   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]);
 };
 
 
@@ -1202,10 +1198,7 @@ FullMatrix<number>::norm2 () const
 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);
 };
 
 
index 00a4b7248764f164ce8556ee2b374a748828b282..020600c685d99517645003a5fd28d047582c1a47 100644 (file)
@@ -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<number>::diag_element (const unsigned int i) const {
 
 
 
+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 {
index 4688022c34b4e396d4f4656ac74e064e5bb418ae..54565fb68522bcfaa7becf96d5e8d318fcc32610 100644 (file)
@@ -143,12 +143,8 @@ SparseMatrix<number>::copy_from (const SparseMatrix<somenumber> &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;
index d5d6f360f4bcfd7e4ead1b5972972264dc89c7d1..f579a75c5721f1095597e268f55b31fb24cce2aa 100644 (file)
@@ -530,9 +530,9 @@ Vector<Number>& Vector<Number>::operator = (const Vector<Number>& v)
 {
   if (v.dim != dim)
     reinit (v.dim, true);
-
   if (dim!=0)
     copy (v.begin(), v.end(), begin());
+  
   return *this;
 }
 
index 5e16cd89c9548fc3f43bd7f5e2c2f5b1d74471bf..e52db7de6fd0c3f794917113d72199d442b30830 100644 (file)
@@ -14,7 +14,7 @@
 #include <iomanip>
 #include <algorithm>
 #include <cmath>
-
+#include <numeric>
 
 
 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<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 ()
 {
@@ -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<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]);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.