]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
The function TrilinosWrappers::SparseMatrix::add() that sets more than one element...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Nov 2008 16:25:46 +0000 (16:25 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Nov 2008 16:25:46 +0000 (16:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@17697 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/source/trilinos_sparse_matrix.cc

index 2b3fb37901883d9540a892c0fb9af784d9e527d8..92db2a99aa3fc2d76d85b2dd9dda16292fd27b74 100644 (file)
@@ -305,10 +305,20 @@ namespace TrilinosWrappers
                                         * FullMatrix<double> into the sparse
                                         * matrix, according to row_indices
                                         * and col_indices.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
                                        */
       void set (const std::vector<unsigned int>  &row_indices,
                const std::vector<unsigned int>  &col_indices,
-               const FullMatrix<TrilinosScalar> &full_matrix);
+               const FullMatrix<TrilinosScalar> &full_matrix,
+               const bool                        elide_zero_values = false);
 
                                        /**
                                         * Set several elements in the
@@ -316,25 +326,42 @@ namespace TrilinosWrappers
                                         * column indices as given by
                                         * <tt>col_indices</tt> to the
                                         * respective value.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
                                        */
-      void set (const unsigned int                row,
+      void set (const unsigned int                 row,
                const std::vector<unsigned int>   &col_indices,
-               const std::vector<TrilinosScalar> &values);
+               const std::vector<TrilinosScalar> &values,
+               const bool                         elide_zero_values = false);
 
                                        /**
                                         * Set several elements to values
-                                        * given by <tt>values</tt> at
-                                        * locations specified by row_indices
-                                        * and col_indices in the sparse
-                                        * matrix. The array <tt>values</tt>
-                                        * is assumed to store data in C
-                                        * style, i.e., row-wise.
+                                        * given by <tt>values</tt> in the
+                                        * given row at column locations
+                                        * specified by col_indices in the
+                                        * sparse matrix.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
                                        */
-      void set (const unsigned int    n_rows,
-               const unsigned int   *row_indices,
+      void set (const unsigned int    row,
                const unsigned int    n_cols,
                const unsigned int   *col_indices,
-               const TrilinosScalar *values);
+               const TrilinosScalar *values,
+               const bool            elide_zero_values = false);
 
                                        /**
                                         * Add @p value to the element
@@ -350,10 +377,21 @@ namespace TrilinosWrappers
                                         * into the sparse matrix, at
                                         * locations according to row_indices
                                         * and col_indices.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be added anyway or
+                                       * these should be filtered away and
+                                       * only non-zero data is added. The
+                                       * default value is <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
                                        */
       void add (const std::vector<unsigned int>  &row_indices,
                const std::vector<unsigned int>  &col_indices,
-               const FullMatrix<TrilinosScalar> &full_matrix);
+               const FullMatrix<TrilinosScalar> &full_matrix,
+               const bool                        elide_zero_values = true);
 
                                        /**
                                         * Set several elements in the
@@ -361,25 +399,44 @@ namespace TrilinosWrappers
                                         * column indices as given by
                                         * <tt>col_indices</tt> to the
                                         * respective value.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be added anyway or
+                                       * these should be filtered away and
+                                       * only non-zero data is added. The
+                                       * default value is <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
                                        */
-      void add (const unsigned int                row,
+      void add (const unsigned int                 row,
                const std::vector<unsigned int>   &col_indices,
-               const std::vector<TrilinosScalar> &values);
+               const std::vector<TrilinosScalar> &values,
+               const bool                         elide_zero_values = true);
 
                                        /**
                                         * Set several elements to values
-                                        * given by <tt>values</tt> at
-                                        * locations specified by row_indices
-                                        * and col_indices in the sparse
-                                        * matrix. The array <tt>values</tt>
-                                        * is assumed to store data in C
-                                        * style, i.e., row-wise.
+                                        * given by <tt>values</tt> in the
+                                        * given row at column locations
+                                        * specified by col_indices in the
+                                        * sparse matrix.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be added anyway or
+                                       * these should be filtered away and
+                                       * only non-zero data is added. The
+                                       * default value is <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
                                        */
-      void add (const unsigned int    n_rows,
-               const unsigned int   *row_indices,
+      void add (const unsigned int    row,
                const unsigned int    n_cols,
                const unsigned int   *col_indices,
-               const TrilinosScalar *values);
+               const TrilinosScalar *values,
+               const bool            elide_zero_values = true);
 
                                        /**
                                         * Matrix-vector multiplication:
@@ -720,7 +777,7 @@ namespace TrilinosWrappers
                                        * writing data from local to global
                                        * matrices.
                                        */
-      std::vector<unsigned int> block_col_indices, local_row_length, local_col_indices;
+      std::vector<unsigned int> block_col_indices, local_col_indices, local_row_length;
   };
 
 
@@ -762,15 +819,17 @@ namespace TrilinosWrappers
   void
   BlockSparseMatrix::set (const std::vector<unsigned int>  &row_indices,
                          const std::vector<unsigned int>  &col_indices,
-                         const FullMatrix<TrilinosScalar> &values)
+                         const FullMatrix<TrilinosScalar> &values,
+                         const bool                        elide_zero_values)
   {
     Assert (row_indices.size() == values.m(),
            ExcDimensionMismatch(row_indices.size(), values.m()));
     Assert (col_indices.size() == values.n(),
            ExcDimensionMismatch(col_indices.size(), values.n()));
 
-    set (row_indices.size(), &row_indices[0], 
-        col_indices.size(), &col_indices[0], &values(0,0));
+    for (unsigned int i=0; i<row_indices.size(); ++i)
+      set (row_indices[i], col_indices.size(), &col_indices[0], &values(0,0),
+          elide_zero_values);
   }
 
 
@@ -779,27 +838,32 @@ namespace TrilinosWrappers
   void
   BlockSparseMatrix::set (const unsigned int                 row,
                          const std::vector<unsigned int>   &col_indices,
-                         const std::vector<TrilinosScalar> &values)
+                         const std::vector<TrilinosScalar> &values,
+                         const bool                         elide_zero_values)
   {
     Assert (col_indices.size() == values.size(),
            ExcDimensionMismatch(col_indices.size(), values.size()));
 
-    set (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+    set (row, col_indices.size(), &col_indices[0], &values[0],
+        elide_zero_values);
   }
 
 
 
+                                  // This is a very messy function, since
+                                  // we need to calculate to each position 
+                                  // the location in the global array.
   inline
   void
-  BlockSparseMatrix::set (const unsigned int    n_rows,
-                         const unsigned int   *row_indices,
+  BlockSparseMatrix::set (const unsigned int    row,
                          const unsigned int    n_cols,
                          const unsigned int   *col_indices,
-                         const TrilinosScalar *values)
+                         const TrilinosScalar *values,
+                         const bool            elide_zero_values)
   {
                                   // Resize scratch arrays
+    local_row_length.resize (this->n_block_cols());
     block_col_indices.resize (this->n_block_cols());
-    local_row_length.resize (this->n_block_cols());    
     local_col_indices.resize (n_cols);
 
                                   // Clear the content in local_row_length
@@ -809,10 +873,7 @@ namespace TrilinosWrappers
                                   // Go through the column indices to find
                                   // out which portions of the values
                                   // should be written into which block
-                                  // matrix. This can be done before
-                                  // starting the loop over all the rows,
-                                  // since we assume a rectangular set of
-                                  // matrix data.
+                                  // matrix.
     {
       unsigned int current_block = 0, row_length = 0;
       block_col_indices[0] = 0;
@@ -853,21 +914,19 @@ namespace TrilinosWrappers
                                   // where we should start reading out
                                   // data. Now let's write the data into
                                   // the individual blocks!
-    for (unsigned int i=0; i<n_rows; ++i)
+    const std::pair<unsigned int,unsigned int> 
+      row_index = this->row_block_indices.global_to_local (row);
+    for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
       {
-       const std::pair<unsigned int,unsigned int> 
-         row_index = this->row_block_indices.global_to_local (row_indices[i]);
-       for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
-         {
-           if (local_row_length[block_col] == 0)
-             continue;
+       if (local_row_length[block_col] == 0)
+         continue;
 
-           block(row_index.first, block_col).set 
-             (1, &row_index.second, 
+       block(row_index.first, block_col).set 
+             (row_index.second, 
               local_row_length[block_col], 
               &local_col_indices[block_col_indices[block_col]],
-              &values[n_cols*i + block_col_indices[block_col]]);
-         }
+              &values[block_col_indices[block_col]],
+              elide_zero_values);
       }
   }
 
@@ -892,15 +951,17 @@ namespace TrilinosWrappers
   void
   BlockSparseMatrix::add (const std::vector<unsigned int>  &row_indices,
                          const std::vector<unsigned int>  &col_indices,
-                         const FullMatrix<TrilinosScalar> &values)
+                         const FullMatrix<TrilinosScalar> &values,
+                         const bool                        elide_zero_values)
   {
     Assert (row_indices.size() == values.m(),
            ExcDimensionMismatch(row_indices.size(), values.m()));
     Assert (col_indices.size() == values.n(),
            ExcDimensionMismatch(col_indices.size(), values.n()));
 
-    add (row_indices.size(), &row_indices[0], 
-        col_indices.size(), &col_indices[0], &values(0,0));
+    for (unsigned int i=0; i<row_indices.size(); ++i)
+      add(row_indices[i], col_indices.size(), &col_indices[0], 
+         &values(0,0), elide_zero_values);
   }
 
 
@@ -909,23 +970,25 @@ namespace TrilinosWrappers
   void
   BlockSparseMatrix::add (const unsigned int                 row,
                          const std::vector<unsigned int>   &col_indices,
-                         const std::vector<TrilinosScalar> &values)
+                         const std::vector<TrilinosScalar> &values,
+                         const bool                         elide_zero_values)
   {
     Assert (col_indices.size() == values.size(),
            ExcDimensionMismatch(col_indices.size(), values.size()));
 
-    add (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+    add (row, col_indices.size(), &col_indices[0], &values[0], 
+        elide_zero_values);
   }
 
 
 
   inline
   void
-  BlockSparseMatrix::add (const unsigned int    n_rows,
-                         const unsigned int   *row_indices,
+  BlockSparseMatrix::add (const unsigned int    row,
                          const unsigned int    n_cols,
                          const unsigned int   *col_indices,
-                         const TrilinosScalar *values)
+                         const TrilinosScalar *values,
+                         const bool            elide_zero_values)
   {
                                   // TODO: Look over this to find out
                                   // whether we can do that more
@@ -943,10 +1006,7 @@ namespace TrilinosWrappers
                                   // Go through the column indices to find
                                   // out which portions of the values
                                   // should be written into which block
-                                  // matrix. This can be done before
-                                  // starting the loop over all the rows,
-                                  // since we assume a rectangular set of
-                                  // matrix data.
+                                  // matrix.
     {
       unsigned int current_block = 0, row_length = 0;
       block_col_indices[0] = 0;
@@ -987,21 +1047,19 @@ namespace TrilinosWrappers
                                   // where we should start reading out
                                   // data. Now let's write the data into
                                   // the individual blocks!
-    for (unsigned int i=0; i<n_rows; ++i)
+    const std::pair<unsigned int,unsigned int> 
+      row_index = this->row_block_indices.global_to_local (row);
+    for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
       {
-       const std::pair<unsigned int,unsigned int> 
-         row_index = this->row_block_indices.global_to_local (row_indices[i]);
-       for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
-         {
-           if (local_row_length[block_col] == 0)
-             continue;
+       if (local_row_length[block_col] == 0)
+         continue;
 
-           block(row_index.first, block_col).add 
-             (1, &row_index.second, 
+       block(row_index.first, block_col).add 
+             (row_index.second, 
               local_row_length[block_col], 
               &local_col_indices[block_col_indices[block_col]],
-              &values[n_cols*i + block_col_indices[block_col]]);
-         }
+              &values[block_col_indices[block_col]],
+              elide_zero_values);
       }
   }
 
index d4a7b18ce6f1779c826d2e0bbfef47fb18149ac9..fb226140018c2ea5fcb93ea4718df85b03ee02d1 100755 (executable)
@@ -883,10 +883,20 @@ namespace TrilinosWrappers
                                        * insertion of elements at positions
                                        * which have not been initialized
                                        * will throw an exception.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
                                        */
       void set (const std::vector<unsigned int>  &row_indices,
                const std::vector<unsigned int>  &col_indices,
-               const FullMatrix<TrilinosScalar> &full_matrix);
+               const FullMatrix<TrilinosScalar> &full_matrix,
+               const bool                        elide_zero_values = false);
 
                                        /**
                                         * Set several elements in the
@@ -905,19 +915,27 @@ namespace TrilinosWrappers
                                        * insertion of elements at positions
                                        * which have not been initialized
                                        * will throw an exception.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
                                        */
       void set (const unsigned int                row,
                const std::vector<unsigned int>   &col_indices,
-               const std::vector<TrilinosScalar> &values);
+               const std::vector<TrilinosScalar> &values,
+               const bool                         elide_zero_values = false);
 
                                        /**
                                         * Set several elements to values
-                                        * given by <tt>values</tt> at
-                                        * locations specified by row_indices
-                                        * and col_indices in the sparse
-                                        * matrix. The array <tt>values</tt>
-                                        * is assumed to store data in C
-                                        * style, i.e., row-wise.
+                                        * given by <tt>values</tt> in a
+                                        * given row in columns given by
+                                        * col_indices into the sparse
+                                        * matrix.
                                        *
                                        * This function is able to insert
                                        * new elements into the matrix as
@@ -929,12 +947,21 @@ namespace TrilinosWrappers
                                        * insertion of elements at positions
                                        * which have not been initialized
                                        * will throw an exception.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be inserted anyway
+                                       * or they should be filtered
+                                       * away. The default value is
+                                       * <tt>false</tt>, i.e., even zero
+                                       * values are inserted/replaced.
                                        */
-      void set (const unsigned int    n_rows,
-               const unsigned int   *row_indices,
+      void set (const unsigned int    row,
                const unsigned int    n_cols,
                const unsigned int   *col_indices,
-               const TrilinosScalar *values);
+               const TrilinosScalar *values,
+               const bool            elide_zero_values = false);
 
                                        /**
                                         * Add @p value to the element
@@ -970,10 +997,21 @@ namespace TrilinosWrappers
                                        * throws an exception if an
                                        * entry does not exist in the
                                        * sparsity pattern.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be added anyway or
+                                       * these should be filtered away and
+                                       * only non-zero data is added. The
+                                       * default value is <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
                                        */
       void add (const std::vector<unsigned int>  &row_indices,
                const std::vector<unsigned int>  &col_indices,
-               const FullMatrix<TrilinosScalar> &full_matrix);
+               const FullMatrix<TrilinosScalar> &full_matrix,
+               const bool                        elide_zero_values = true);
 
                                        /**
                                         * Set several elements in the
@@ -990,34 +1028,52 @@ namespace TrilinosWrappers
                                        * throws an exception if an
                                        * entry does not exist in the
                                        * sparsity pattern.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be added anyway or
+                                       * these should be filtered away and
+                                       * only non-zero data is added. The
+                                       * default value is <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
                                        */
       void add (const unsigned int                row,
                const std::vector<unsigned int>   &col_indices,
-               const std::vector<TrilinosScalar> &values);
+               const std::vector<TrilinosScalar> &values,
+               const bool                         elide_zero_values = true);
 
                                        /**
-                                        * Set several elements to values
-                                        * given by <tt>values</tt> at
-                                        * locations specified by row_indices
-                                        * and col_indices in the sparse
-                                        * matrix. The array <tt>values</tt>
-                                        * is assumed to store data in C
-                                        * style, i.e., row-wise.
+                                        * Add an array of values given by
+                                        * <tt>values</tt> in the given
+                                        * global matrix row at columns
+                                        * specified by col_indices in the
+                                        * sparse matrix.
                                        *
                                        * Just as the respective call in
-                                       * deal.II SparseMatrix<Number>
-                                       * class (but in contrast to the
-                                       * situation for PETSc based
-                                       * matrices), this function
-                                       * throws an exception if an
+                                       * deal.II SparseMatrix<Number> class
+                                       * (but in contrast to the situation
+                                       * for PETSc based matrices), this
+                                       * function throws an exception if an
                                        * entry does not exist in the
                                        * sparsity pattern.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be added anyway or
+                                       * these should be filtered away and
+                                       * only non-zero data is added. The
+                                       * default value is <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
                                        */
-      void add (const unsigned int    n_rows,
-               const unsigned int   *row_indices,
+      void add (const unsigned int    row,
                const unsigned int    n_cols,
                const unsigned int   *col_indices,
-               const TrilinosScalar *values);
+               const TrilinosScalar *values,
+               const bool            elide_zero_values = false);
 
                                        /**
                                         * Multiply the entire matrix
@@ -1683,6 +1739,24 @@ namespace TrilinosWrappers
                                        */
       mutable VectorBase temp_vector;
 
+                                      /**
+                                       * An internal array of integer
+                                       * values that is used to store the
+                                       * column indices when
+                                       * adding/inserting local data into
+                                       * the (large) sparse matrix.
+                                       */
+      std::vector<unsigned int> column_indices;
+
+                                      /**
+                                       * An internal array of double values
+                                       * that is used to store the column
+                                       * indices when adding/inserting
+                                       * local data into the (large) sparse
+                                       * matrix.
+                                       */
+      std::vector<TrilinosScalar> column_values;
+
     public:
                                        /**
                                         * A sparse matrix object in
@@ -1954,7 +2028,7 @@ namespace TrilinosWrappers
            ExcMessage("The given value is not finite but either "
                       "infinite or Not A Number (NaN)"));
 
-    set (1, &i, 1, &j, &value);
+    set (i, 1, &j, &value, false);
   }
 
 
@@ -1963,15 +2037,17 @@ namespace TrilinosWrappers
   void
   SparseMatrix::set (const std::vector<unsigned int>  &row_indices,
                     const std::vector<unsigned int>  &col_indices,
-                    const FullMatrix<TrilinosScalar> &values)
+                    const FullMatrix<TrilinosScalar> &values,
+                    const bool                        elide_zero_values)
   {
     Assert (row_indices.size() == values.m(),
            ExcDimensionMismatch(row_indices.size(), values.m()));
     Assert (col_indices.size() == values.n(),
            ExcDimensionMismatch(col_indices.size(), values.n()));
 
-    set (row_indices.size(), &row_indices[0], 
-        col_indices.size(), &col_indices[0], &values(0,0));
+    for (unsigned int i=0; i<row_indices.size(); ++i)
+      set (row_indices[i], col_indices.size(), &col_indices[0], &values(0,0),
+          elide_zero_values);
   }
 
 
@@ -1980,23 +2056,25 @@ namespace TrilinosWrappers
   void
   SparseMatrix::set (const unsigned int                 row,
                     const std::vector<unsigned int>   &col_indices,
-                    const std::vector<TrilinosScalar> &values)
+                    const std::vector<TrilinosScalar> &values,
+                    const bool                         elide_zero_values)
   {
     Assert (col_indices.size() == values.size(),
            ExcDimensionMismatch(col_indices.size(), values.size()));
 
-    set (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+    set (row, col_indices.size(), &col_indices[0], &values[0],
+        elide_zero_values);
   }
 
 
 
   inline
   void
-  SparseMatrix::set (const unsigned int    n_rows,
-                    const unsigned int   *row_indices,
+  SparseMatrix::set (const unsigned int    row,
                     const unsigned int    n_cols,
                     const unsigned int   *col_indices,
-                    const TrilinosScalar *values)
+                    const TrilinosScalar *values,
+                    const bool            elide_zero_values)
   {
     int ierr;
     if (last_action == Add)
@@ -2009,11 +2087,38 @@ namespace TrilinosWrappers
     if (last_action != Insert)
       last_action = Insert;
 
-                                  // Now go through all rows that are
-                                  // present in the input data.
-    for (unsigned int i=0; i<n_rows; ++i)
+    int * col_index_ptr;
+    TrilinosScalar const* col_value_ptr;
+
+                                  // If we don't elide zeros, the pointers
+                                  // are already available...
+    if (elide_zero_values == false)
+      {
+       col_index_ptr = (int*)col_indices;
+       col_value_ptr = values;
+      }
+    else
       {
-       const int row = row_indices[i];
+                                  // Otherwise, extract nonzero values in
+                                  // each row and pass on to the other 
+                                  // function.
+       column_indices.resize(n_cols);
+       column_values.resize(n_cols);
+
+       unsigned int n_columns = 0;
+       for (unsigned int j=0; j<n_cols; ++j)
+         if (values[j] != 0)
+           {
+             column_indices[n_columns] = col_indices[j];
+             column_values[n_columns] = values[j];
+             n_columns++;
+           }
+       Assert(n_columns <= n_cols, ExcInternalError());
+
+       col_index_ptr = (int*)&column_indices[0];
+       col_value_ptr = &column_values[0];
+      }
+
 
                                   // If the calling matrix owns the row to
                                   // which we want to insert values, we
@@ -2029,31 +2134,29 @@ namespace TrilinosWrappers
                                   // add the possibility to insert new values,
                                   // and in the second we just replace
                                   // data.
-       if (row_map.MyGID(row) == true)
+    if (row_map.MyGID(row) == true)
+      {
+       if (matrix->Filled() == false)
          {
-           if (matrix->Filled() == false)
-             {
-               ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
-                                  row, (int)n_cols,
-                                  const_cast<TrilinosScalar*>(&values[i*n_cols]),
-                                  (int*)&col_indices[0]);
+           ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_cols, 
+                                           const_cast<double*>(col_value_ptr),
+                                                               col_index_ptr);
 
                                        // When inserting elements, we do
                                        // not want to create exceptions in
                                        // the case when inserting non-local
                                        // data (since that's what we want 
                                        // to do right now).
-               if (ierr > 0)
-                 ierr = 0;
-             }
-           else
-             ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(
-                                  row, (int)n_cols,
-                                  const_cast<TrilinosScalar*>(&values[i*n_cols]),
-                                  (int*)&col_indices[0]);
+           if (ierr > 0)
+             ierr = 0;
          }
        else
-         {
+         ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row, n_cols,     
+                                          const_cast<double*>(col_value_ptr),
+                                                              col_index_ptr);
+      }
+    else
+      {
                                   // When we're at off-processor data, we
                                   // have to stick with the standard
                                   // Insert/ReplaceGlobalValues
@@ -2064,30 +2167,26 @@ namespace TrilinosWrappers
                                   // call the function we already use,
                                   // which is very unefficient if writing
                                   // one element at a time).
-           compressed = false;
+       compressed = false;
       
-           const TrilinosScalar* value_ptr = &values[i*n_cols];
-
-           if (matrix->Filled() == false)
-             {
-               ierr = matrix->InsertGlobalValues (1, &row, 
-                                           (int)n_cols, (int*)&col_indices[0],
-                                           &value_ptr, 
-                                           Epetra_FECrsMatrix::ROW_MAJOR);
-               if (ierr > 0)
-                 ierr = 0;
-             }
-           else
-             ierr = matrix->ReplaceGlobalValues (1, &row, 
-                                           (int)n_cols, (int*)&col_indices[0],
-                                           &value_ptr, 
-                                           Epetra_FECrsMatrix::ROW_MAJOR);
+       if (matrix->Filled() == false)
+         {
+           ierr = matrix->InsertGlobalValues (1, (int*)&row, 
+                                              (int)n_cols, col_index_ptr,
+                                              &col_value_ptr, 
+                                              Epetra_FECrsMatrix::ROW_MAJOR);
+           if (ierr > 0)
+             ierr = 0;
          }
-
-       Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i],
-                                                       col_indices[0]));
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+       else
+         ierr = matrix->ReplaceGlobalValues (1, (int*)&row, 
+                                             (int)n_cols, col_index_ptr,
+                                             &col_value_ptr, 
+                                             Epetra_FECrsMatrix::ROW_MAJOR);
       }
+
+    Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
 
 
@@ -2105,7 +2204,8 @@ namespace TrilinosWrappers
 
     if (value == 0)
       {
-                                 // we have to do above actions in any case
+                                 // we have to do checkings on Insert/Add
+                                 // in any case
                                  // to be consistent with the MPI
                                  // communication model (see the comments
                                  // in the documentation of
@@ -2128,7 +2228,7 @@ namespace TrilinosWrappers
        return;
       }
     else
-      add (1, &i, 1, &j, &value);
+      add (i, 1, &j, &value, false);
   }
 
 
@@ -2137,15 +2237,17 @@ namespace TrilinosWrappers
   void
   SparseMatrix::add (const std::vector<unsigned int>  &row_indices,
                     const std::vector<unsigned int>  &col_indices,
-                    const FullMatrix<TrilinosScalar> &values)
+                    const FullMatrix<TrilinosScalar> &values,
+                    const bool                        elide_zero_values)
   {
     Assert (row_indices.size() == values.m(),
            ExcDimensionMismatch(row_indices.size(), values.m()));
     Assert (col_indices.size() == values.n(),
            ExcDimensionMismatch(col_indices.size(), values.n()));
 
-    add (row_indices.size(), &row_indices[0], 
-        col_indices.size(), &col_indices[0], &values(0,0));
+    for (unsigned int i=0; i<row_indices.size(); ++i)
+      add (row_indices[i], col_indices.size(), &col_indices[0], &values(0,0),
+          elide_zero_values);
   }
 
 
@@ -2154,23 +2256,25 @@ namespace TrilinosWrappers
   void
   SparseMatrix::add (const unsigned int                 row,
                     const std::vector<unsigned int>   &col_indices,
-                    const std::vector<TrilinosScalar> &values)
+                    const std::vector<TrilinosScalar> &values,
+                    const bool                         elide_zero_values)
   {
     Assert (col_indices.size() == values.size(),
            ExcDimensionMismatch(col_indices.size(), values.size()));
 
-    add (1, &row, col_indices.size(), &col_indices[0], &values[0]);
+    add (row, col_indices.size(), &col_indices[0], &values[0],
+        elide_zero_values);
   }
 
 
 
   inline
   void
-  SparseMatrix::add (const unsigned int    n_rows,
-                    const unsigned int   *row_indices,
+  SparseMatrix::add (const unsigned int    row,
                     const unsigned int    n_cols,
                     const unsigned int   *col_indices,
-                    const TrilinosScalar *values)
+                    const TrilinosScalar *values,
+                    const bool            elide_zero_values)
   {
     int ierr;
     if (last_action == Insert)
@@ -2183,25 +2287,51 @@ namespace TrilinosWrappers
     if (last_action != Add)
       last_action = Add;
 
-                                  // Now go through all rows that are
-                                  // present in the input data.
-    for (unsigned int i=0; i<n_rows; ++i)
+    int * col_index_ptr;
+    TrilinosScalar const* col_value_ptr;
+
+                                  // If we don't elide zeros, the pointers
+                                  // are already available...
+    if (elide_zero_values == false)
       {
-       const int row = row_indices[i];
+       col_index_ptr = (int*)&col_indices[0];
+       col_value_ptr = values;
+      }
+    else
+      {
+                                  // Otherwise, extract nonzero values in
+                                  // each row and pass on to the other 
+                                  // function.
+       column_indices.resize(n_cols);
+       column_values.resize(n_cols);
+
+       unsigned int n_columns = 0;
+       for (unsigned int j=0; j<n_cols; ++j)
+         if (values[j] != 0)
+           {
+             column_indices[n_columns] = col_indices[j];
+             column_values[n_columns] = values[j];
+             n_columns++;
+           }
+       Assert(n_columns <= n_cols, ExcInternalError());
+
+       col_index_ptr = (int*)&column_indices[0];
+       col_value_ptr = &column_values[0];
+      }
+
                                   // If the calling matrix owns the row to
                                   // which we want to add values, we
                                   // can directly call the Epetra_CrsMatrix
                                   // input function, which is much faster
                                   // than the Epetra_FECrsMatrix function.
-       if (row_map.MyGID(row_indices[i]) == true)
-         {
-           ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(
-                                  row, (int)n_cols,
-                                  const_cast<TrilinosScalar*>(&values[i*n_cols]),
-                                  (int*)&col_indices[0]);
-         }
-       else
-         {
+    if (row_map.MyGID(row) == true)
+      {
+       ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_cols,
+                                        const_cast<double*>(col_value_ptr),
+                                                            col_index_ptr);
+      }
+    else
+      {
                                   // When we're at off-processor data, we
                                   // have to stick with the standard
                                   // SumIntoGlobalValues
@@ -2212,20 +2342,16 @@ namespace TrilinosWrappers
                                   // call the function we already use,
                                   // which is very unefficient if writing
                                   // one element at a time).
-           compressed = false;
-      
-           const TrilinosScalar* value_ptr = &values[i*n_cols];
-
-           ierr = matrix->SumIntoGlobalValues (1, &row, 
-                                               (int)n_cols, (int*)&col_indices[0],
-                                               &value_ptr, 
-                                               Epetra_FECrsMatrix::ROW_MAJOR);
-         }
+       compressed = false;
 
-       Assert (ierr <= 0, ExcAccessToNonPresentElement(row_indices[i],
-                                                       col_indices[0]));
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+       ierr = matrix->SumIntoGlobalValues (1, (int*)&row, (int)n_cols, 
+                                           col_index_ptr,
+                                           &col_value_ptr, 
+                                           Epetra_FECrsMatrix::ROW_MAJOR);
       }
+
+    Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
 
 
index 18570f3f9db739bd6446c3bc11f7f511f0834652..d45ebca374019c2b4a9b561369006b498aeeff56 100755 (executable)
@@ -316,7 +316,7 @@ namespace TrilinosWrappers
 //                                   sparsity_pattern.n_cols()));
 
     std::vector<TrilinosScalar> values;
-    std::vector<int>            row_indices;
+    std::vector<unsigned int>   row_indices;
     
     for (unsigned int row=0; row<n_rows; ++row)
       if (row_map.MyGID(row))
@@ -327,10 +327,8 @@ namespace TrilinosWrappers
 
          for (int col=0; col < row_length; ++col)
            row_indices[col] = sparsity_pattern.column_number (row, col);
-
-         matrix->Epetra_CrsMatrix::InsertGlobalValues ((int)row, row_length,
-                                                       &values[0],
-                                                       &row_indices[0]);
+         
+         set (row, row_length, &row_indices[0], &values[0], false);
        }
     
     last_action = Zero;
@@ -416,7 +414,7 @@ namespace TrilinosWrappers
 //                                   sparsity_pattern.n_cols()));
 
     std::vector<TrilinosScalar> values;
-    std::vector<int>    row_indices;
+    std::vector<unsigned int>   row_indices;
     
     for (unsigned int row=0; row<n_rows; ++row)
       if (row_map.MyGID(row))
@@ -433,9 +431,7 @@ namespace TrilinosWrappers
               ++col_num, ++col)
            row_indices[col] = *col_num;
 
-         matrix->Epetra_CrsMatrix::InsertGlobalValues ((int)row, row_length,
-                                                       &values[0],
-                                                       &row_indices[0]);
+         set (row, row_length, &row_indices[0], &values[0], false);
        }
     
     last_action = Zero;
@@ -536,7 +532,7 @@ namespace TrilinosWrappers
                                        &n_entries_per_row[0], false));
 
     std::vector<TrilinosScalar> values;
-    std::vector<int> row_indices;
+    std::vector<unsigned int>   row_indices;
 
     for (unsigned int row=0; row<n_rows; ++row)
       {
@@ -554,11 +550,7 @@ namespace TrilinosWrappers
              ++index;
            }
 
-       const int row_length = index;
-
-       matrix->Epetra_CrsMatrix::InsertGlobalValues ((int)row, row_length,
-                                                     &values[0],
-                                                     &row_indices[0]);
+       set (row, index, &row_indices[0], &values[0], false);
       }
 
     compress();

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.