]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Some minor changes.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Sep 2008 09:40:49 +0000 (09:40 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Sep 2008 09:40:49 +0000 (09:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@16773 0785d39b-7218-0410-832d-ea1e28bc413d

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

index ca2eaf08a837219f891ddd87e17e1b6ab6be88c0..1f037f175b128cab90a4cb0467cc8f8efa3aad91 100755 (executable)
@@ -334,7 +334,9 @@ namespace TrilinosWrappers
       typedef TrilinosScalar value_type;
       
                                        /**
-                                        * Default constructor.
+                                        * Default
+                                        * constructor. Generates an
+                                        * empty (zero-size) matrix.
                                         */
       SparseMatrix ();
 
@@ -491,7 +493,7 @@ namespace TrilinosWrappers
                   const Epetra_Map       &input_col_map,
                   const SparsityPattern  &sparsity_pattern);
 
-                                      /** 
+                                      /**
                                        * This function copies the
                                        * content in
                                        * <tt>sparse_matrix</tt> to
@@ -642,11 +644,12 @@ namespace TrilinosWrappers
                                         * the matrix without having to
                                         * read entries (such as the
                                         * locations of non-zero
-                                        * elements) from it -- without
-                                        * this o peration, removing
-                                        * constraints on parallel
-                                        * matrices is a rather
-                                        * complicated procedure.
+                                        * elements) from it &mdash;
+                                        * without this operation,
+                                        * removing constraints on
+                                        * parallel matrices is a
+                                        * rather complicated
+                                        * procedure.
                                         *
                                         * The second parameter can be
                                         * used to set the diagonal
@@ -683,20 +686,17 @@ namespace TrilinosWrappers
                                         * may be an expensive
                                         * operation and you should
                                         * always take care where to
-                                        * call this function. In
-                                        * contrast to the respective
-                                        * function in the @p
-                                        * SparseMatrix class, we don't
-                                        * throw an exception if the
-                                        * respective entry doesn't
-                                        * exist in the sparsity
-                                        * pattern of this class, since
-                                        * Trilinos does not transmit
-                                        * this information.  On the
-                                        * other hand, an exception
-                                        * will be thrown when the
-                                        * requested element is not
-                                        * saved on the calling
+                                        * call this function. As in
+                                        * the &p SparseMatrix<Number>
+                                        * class, we throw an exception
+                                        * if the respective entry
+                                        * doesn't exist in the
+                                        * sparsity pattern of this
+                                        * class, which is requested
+                                        * from Trilinos. Moreover, an
+                                        * exception will be thrown
+                                        * when the requested element
+                                        * is not saved on the calling
                                         * process.
                                         *
                                         * This function is therefore
@@ -1163,6 +1163,14 @@ namespace TrilinosWrappers
                       << "An error with error number " << arg1
                       << " occured while calling a Trilinos function");
 
+                                      /**
+                                       * Exception
+                                       */
+      DeclException2 (ExcInvalidIndex,
+                     int, int,
+                     << "The entry with index <" << arg1 << ',' << arg2
+                     << "> does not exist.");
+
                                        /**
                                         * Exception
                                         */
@@ -1248,6 +1256,9 @@ namespace TrilinosWrappers
                                         * elements.  The actual type,
                                         * a sparse matrix, is set in
                                         * the constructor.
+                                       *
+                                       * TODO: This object should
+                                       * finally become private.
                                         */
       std::auto_ptr<Epetra_FECrsMatrix> matrix;
 
@@ -1285,6 +1296,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     unsigned int
     const_iterator::Accessor::column() const
@@ -1294,6 +1306,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     unsigned int
     const_iterator::Accessor::index() const
@@ -1303,6 +1316,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     TrilinosScalar
     const_iterator::Accessor::value() const
@@ -1312,6 +1326,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     const_iterator::
     const_iterator(const SparseMatrix *matrix,
@@ -1351,6 +1366,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     const_iterator
     const_iterator::operator++ (int)
@@ -1361,6 +1377,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     const const_iterator::Accessor &
     const_iterator::operator* () const
@@ -1369,6 +1386,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     const const_iterator::Accessor *
     const_iterator::operator-> () const
@@ -1377,6 +1395,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     bool
     const_iterator::
@@ -1387,6 +1406,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     bool
     const_iterator::
@@ -1396,6 +1416,7 @@ namespace TrilinosWrappers
     }
 
 
+
     inline
     bool
     const_iterator::
@@ -1407,17 +1428,8 @@ namespace TrilinosWrappers
     }
     
   }
-  
-  
-  inline
-  TrilinosScalar
-  SparseMatrix::operator() (const unsigned int i,
-                          const unsigned int j) const
-  {
-    return el(i,j);
-  }
 
-  
 
   inline
   SparseMatrix::const_iterator
@@ -1427,6 +1439,7 @@ namespace TrilinosWrappers
   }
 
 
+
   inline
   SparseMatrix::const_iterator
   SparseMatrix::end() const
@@ -1435,6 +1448,7 @@ namespace TrilinosWrappers
   }
 
 
+
   inline
   SparseMatrix::const_iterator
   SparseMatrix::begin(const unsigned int r) const
@@ -1447,6 +1461,7 @@ namespace TrilinosWrappers
   }
 
 
+
   inline
   SparseMatrix::const_iterator
   SparseMatrix::end(const unsigned int r) const
index 0d00c4fb091bea3696f20830d1e521c70f103578..9488901c6b1cf88681d6cd708bcbd5afc02f852a 100755 (executable)
@@ -376,9 +376,9 @@ namespace TrilinosWrappers
                                  // flush buffers
     int ierr;
     if (row_map.SameAs(col_map))
-      ierr = matrix->GlobalAssemble ();
+      ierr = matrix->GlobalAssemble (true);
     else
-      ierr = matrix->GlobalAssemble (col_map, row_map);
+      ierr = matrix->GlobalAssemble (col_map, row_map, true);
     
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
@@ -530,6 +530,82 @@ namespace TrilinosWrappers
 
 
 
+  TrilinosScalar
+  SparseMatrix::operator() (const unsigned int i,
+                           const unsigned int j) const
+  {
+                                     // Extract local indices in
+                                     // the matrix.
+    int trilinos_i = matrix->LRID(i), trilinos_j = matrix->LRID(j);
+    TrilinosScalar value = 0.;
+
+                                     // If the data is not on the
+                                     // present processor, we throw
+                                     // an exception. This is on of
+                                     // the two tiny differences to
+                                     // the el(i,j) call, which does
+                                     // not throw any assertions.
+    if ((trilinos_i == -1 ) || (trilinos_j == -1))
+      {
+       Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
+                                                 local_range().second));
+      }
+    else
+      {
+                                     // Check whether the matrix 
+                                     // already is transformed to
+                                     // local indices.
+       if (matrix->Filled() == false)
+         matrix->FillComplete(true);
+
+                                     // Prepare pointers for extraction
+                                     // of a view of the row.
+       int nnz_present = matrix->NumMyEntries(trilinos_i);
+       int nnz_extracted;
+       int *col_indices;
+       TrilinosScalar *values;
+
+                                     // Generate the view and make
+                                     // sure that we have not generated
+                                     // an error.
+       int ierr = matrix->ExtractMyRowView(trilinos_i, nnz_extracted,
+                                           values, col_indices);
+       Assert (ierr==0, ExcTrilinosError(ierr));
+
+       Assert (nnz_present == nnz_extracted,
+               ExcDimensionMismatch(nnz_present, nnz_extracted));
+
+                                     // Search the index where we
+                                     // look for the value, and then
+                                     // finally get it.
+      
+       int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present,
+                                trilinos_j);
+      
+       int el_index = (int)(el_find - col_indices);
+
+                                       // This is actually the only
+                                       // difference to the el(i,j)
+                                       // function, which means that
+                                       // we throw an exception in
+                                       // this case instead of just
+                                       // returning zero for an
+                                       // element that is not present
+                                       // in the sparsity pattern.
+       if (!el_find)
+         {
+           Assert (false, ExcInvalidIndex (i,j));
+         }
+       else
+         value = values[el_index];
+
+      }
+
+    return value;
+  }
+
+
+
   TrilinosScalar
   SparseMatrix::el (const unsigned int i,
                    const unsigned int j) const
@@ -544,7 +620,9 @@ namespace TrilinosWrappers
                                      // continue. Just print out
                                      // zero.
 
-                                     // TODO: Is this reasonable?
+                                     // TODO: Is this reasonable? Or
+                                     // should we retain the assert
+                                     // call?
     if ((trilinos_i == -1 ) || (trilinos_j == -1))
       {
        return 0.;
@@ -689,7 +767,7 @@ namespace TrilinosWrappers
   TrilinosScalar
   SparseMatrix::l1_norm () const
   {
-    if (!matrix->Filled())
+    if (matrix->Filled() == false)
       matrix->FillComplete();
 
     TrilinosScalar result = matrix->NormOne();
@@ -702,7 +780,7 @@ namespace TrilinosWrappers
   TrilinosScalar
   SparseMatrix::linfty_norm () const
   {
-    if (!matrix->Filled())
+    if (matrix->Filled() == false)
       matrix->FillComplete();
 
     TrilinosScalar result = matrix->NormInf();
@@ -715,7 +793,7 @@ namespace TrilinosWrappers
   TrilinosScalar
   SparseMatrix::frobenius_norm () const
   {
-    if (!matrix->Filled())
+    if (matrix->Filled() == false)
       matrix->FillComplete();
 
     TrilinosScalar result = matrix->NormFrobenius();

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.