]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Provide a correct implemention of TrilinosWrappers::VectorBase::size() even when...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 18 May 2009 07:30:16 +0000 (07:30 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 18 May 2009 07:30:16 +0000 (07:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@18864 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 9447da225dd3b0560f9402ca0d9ab9a1587d7583..6f96a1378e016398e5e8d22d5f8c27f711f209f2 100755 (executable)
@@ -2203,6 +2203,41 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  void
+  SparseMatrix::compress ()
+  {
+                                 // flush buffers
+    int ierr;
+    ierr = matrix->GlobalAssemble (col_map, row_map, true);
+    
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = matrix->OptimizeStorage ();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    last_action = Zero;
+
+    compressed = true;
+  }
+
+
+
+  inline
+  SparseMatrix &
+  SparseMatrix::operator = (const double d)
+  {
+    Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
+    compress ();
+
+    const int ierr = matrix->PutScalar(d);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
                                        // Inline the set() and add()
                                        // functions, since they will be 
                                         // called frequently, and the 
@@ -2631,6 +2666,248 @@ namespace TrilinosWrappers
 
 
 
+                                  // inline "simple" functions that are
+                                  // called frequently and do only involve
+                                  // a call to some Trilinos function.
+  inline
+  unsigned int
+  SparseMatrix::m () const
+  {
+    return matrix -> NumGlobalRows();
+  }
+
+
+
+  inline
+  unsigned int
+  SparseMatrix::n () const
+  {
+    return matrix -> NumGlobalCols();
+  }
+
+
+
+  inline
+  unsigned int
+  SparseMatrix::local_size () const
+  {
+    return matrix -> NumMyRows();
+  }
+
+
+
+  inline
+  std::pair<unsigned int, unsigned int>
+  SparseMatrix::local_range () const
+  {
+    unsigned int begin, end;
+    begin = matrix -> RowMap().MinMyGID();
+    end = matrix -> RowMap().MaxMyGID()+1;
+    
+    return std::make_pair (begin, end);
+  }
+
+
+
+  inline
+  unsigned int
+  SparseMatrix::n_nonzero_elements () const
+  {
+    return matrix->NumGlobalNonzeros();
+  }
+
+
+
+  inline
+  TrilinosScalar
+  SparseMatrix::l1_norm () const
+  {
+    if (matrix->Filled() == false)
+      matrix->GlobalAssemble(col_map, row_map, true);
+
+    return matrix->NormOne();
+  }
+  
+  
+
+  inline
+  TrilinosScalar
+  SparseMatrix::linfty_norm () const
+  {
+    if (matrix->Filled() == false)
+      matrix->GlobalAssemble(col_map, row_map, true);
+
+    return matrix->NormInf();
+  }
+
+
+
+  inline
+  TrilinosScalar
+  SparseMatrix::frobenius_norm () const
+  {
+    if (matrix->Filled() == false)
+      matrix->GlobalAssemble(col_map, row_map, true);
+
+    return matrix->NormFrobenius();
+  }
+
+
+
+  inline
+  SparseMatrix &
+  SparseMatrix::operator *= (const TrilinosScalar a)
+  {
+    const int ierr = matrix->Scale (a);
+    Assert (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  inline
+  SparseMatrix &
+  SparseMatrix::operator /= (const TrilinosScalar a)
+  {
+    Assert (a !=0, ExcDivideByZero());
+
+    const TrilinosScalar factor = 1./a;
+
+    const int ierr = matrix->Scale (factor);
+    Assert (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  inline
+  void
+  SparseMatrix::vmult (VectorBase       &dst,
+                      const VectorBase &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+    if (matrix->Filled() == false)
+      matrix->GlobalAssemble(col_map, row_map, true);
+
+    Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true,
+           ExcMessage ("Column map of matrix does not fit with vector map!"));
+    Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true,
+           ExcMessage ("Row map of matrix does not fit with vector map!"));
+
+    const int ierr = matrix->Multiply (false, src.trilinos_vector(), 
+                                      dst.trilinos_vector());
+    Assert (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  SparseMatrix::Tvmult (VectorBase       &dst,
+                       const VectorBase &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+    if (matrix->Filled() == false)
+      matrix->GlobalAssemble(col_map, row_map, true);
+
+    Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true,
+           ExcMessage ("Column map of matrix does not fit with vector map!"));
+    Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true,
+           ExcMessage ("Row map of matrix does not fit with vector map!"));
+
+    const int ierr = matrix->Multiply (true, src.trilinos_vector(), 
+                                      dst.trilinos_vector());
+    Assert (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  SparseMatrix::vmult_add (VectorBase       &dst,
+                          const VectorBase &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+                                  // Choose to reinit the vector with fast
+                                  // argument set, which does not overwrite
+                                  // the content -- this is what we need
+                                  // since we're going to overwrite that
+                                  // anyway in the vmult operation.
+    temp_vector.reinit(dst, true);
+
+    vmult (temp_vector, src);
+    dst += temp_vector;
+  }
+
+
+
+  inline
+  void
+  SparseMatrix::Tvmult_add (VectorBase       &dst,
+                           const VectorBase &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+    temp_vector.reinit(dst, true);
+
+    vmult (temp_vector, src);
+    dst += temp_vector;
+  }
+
+
+
+  inline
+  TrilinosScalar
+  SparseMatrix::matrix_norm_square (const VectorBase &v) const
+  {
+    Assert (row_map.SameAs(col_map),
+           ExcDimensionMismatch(row_map.NumGlobalElements(),
+                                col_map.NumGlobalElements()));
+
+    temp_vector.reinit(v);
+
+    vmult (temp_vector, v);
+    return temp_vector*v;
+  }
+
+
+
+  inline
+  TrilinosScalar
+  SparseMatrix::matrix_scalar_product (const VectorBase &u,
+                                      const VectorBase &v) const
+  {
+    Assert (row_map.SameAs(col_map),
+           ExcDimensionMismatch(row_map.NumGlobalElements(),
+                                col_map.NumGlobalElements()));
+
+    temp_vector.reinit(v);
+
+    vmult (temp_vector, v);
+    return u*temp_vector;
+  }
+
+
+
+  inline
+  TrilinosScalar
+  SparseMatrix::residual (VectorBase       &dst,
+                         const VectorBase &x,
+                         const VectorBase &b) const
+  {
+    vmult (dst, x);
+    dst -= b;
+    dst *= -1.;
+
+    return dst.l2_norm();
+  }
+
+
   inline
   const Epetra_CrsMatrix &
   SparseMatrix::trilinos_matrix () const
index d86ca5113fdcf3c4c59632b9ce916504f50f3015..923e32657bcd3ac87b62988ed5e70ef0f66b4d11 100644 (file)
@@ -1105,6 +1105,53 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  void
+  VectorBase::reinit (const VectorBase &v,
+                     const bool        fast)
+  {
+    Assert (&*vector != 0,
+           ExcMessage("Vector has not been constructed properly."));
+
+    if (fast == false || local_range() != v.local_range())
+      vector = std::auto_ptr<Epetra_FEVector>(new Epetra_FEVector(*v.vector));
+  }
+
+
+
+  inline
+  void
+  VectorBase::compress ()
+  {
+                                // Now pass over the information about
+                                // what we did last to the vector.
+    const int ierr = vector->GlobalAssemble(last_action);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+    last_action = Zero;
+  
+    compressed = true;
+  }
+
+
+
+  inline
+  VectorBase &
+  VectorBase::operator = (const TrilinosScalar s)
+  {
+
+    Assert (numbers::is_finite(s),
+           ExcMessage("The given value is not finite but either "
+                      "infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->PutScalar(s);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
   inline
   void
   VectorBase::set (const std::vector<unsigned int>    &indices,
@@ -1220,6 +1267,24 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  unsigned int
+  VectorBase::size () const
+  {
+    return (unsigned int) vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
+  }
+
+
+
+  inline
+  unsigned int
+  VectorBase::local_size () const
+  {
+    return (unsigned int) vector->Map().NumMyElements();
+  }
+
+
+
   inline
   std::pair<unsigned int, unsigned int>
   VectorBase::local_range () const
@@ -1232,6 +1297,488 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  TrilinosScalar
+  VectorBase::operator * (const VectorBase &vec) const
+  {
+    Assert (local_range() == vec.local_range(),
+           ExcDimensionMismatch(size(), vec.size()));
+
+    TrilinosScalar result;
+
+    const int ierr = vector->Dot(*(vec.vector), &result);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return result;
+  }
+
+
+
+  inline
+  VectorBase::real_type
+  VectorBase::norm_sqr () const
+  {
+    const TrilinosScalar d = l2_norm();
+    return d*d;
+  }
+
+
+
+  inline
+  TrilinosScalar
+  VectorBase::mean_value () const
+  {
+    TrilinosScalar mean;
+
+    const int ierr = vector->MeanValue (&mean);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return mean;
+  }
+
+
+
+  inline
+  VectorBase::real_type
+  VectorBase::l1_norm () const
+  {
+    TrilinosScalar d;
+
+    const int ierr = vector->Norm1 (&d);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return d;
+  }
+
+
+
+  inline
+  VectorBase::real_type
+  VectorBase::l2_norm () const
+  {
+    TrilinosScalar d;
+
+    const int ierr = vector->Norm2 (&d);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return d;
+  }
+
+
+
+  inline
+  VectorBase::real_type
+  VectorBase::lp_norm (const TrilinosScalar p) const
+  {
+                                        // get a representation of the
+                                        // vector and loop over all
+                                        // the elements
+    TrilinosScalar *start_ptr;
+    int leading_dimension;
+    int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    TrilinosScalar norm = 0;
+    TrilinosScalar sum=0;
+
+    const TrilinosScalar * ptr  = start_ptr;
+
+                                       // add up elements 
+                                      // TODO: This
+                                      // won't work in parallel like
+                                      // this. Find out a better way to
+                                      // this in that case.
+    while (ptr != start_ptr+size())
+      sum += std::pow(std::fabs(*ptr++), p);
+
+    norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
+
+    return norm;
+  }
+
+
+
+  inline
+  VectorBase::real_type
+  VectorBase::linfty_norm () const
+  {
+    TrilinosScalar d;
+
+    const int ierr = vector->NormInf (&d);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return d;
+  }
+
+
+
+                                  // inline also scalar products, vector
+                                  // additions etc. since they are all
+                                  // representable by a single Trilinos
+                                  // call. This reduces the overhead of the
+                                  // wrapper class.
+  inline
+  VectorBase &
+  VectorBase::operator *= (const TrilinosScalar a)
+  {
+
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->Scale(a);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  inline
+  VectorBase &
+  VectorBase::operator /= (const TrilinosScalar a)
+  {
+
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const TrilinosScalar factor = 1./a;
+
+    Assert (numbers::is_finite(factor),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->Scale(factor);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  inline
+  VectorBase &
+  VectorBase::operator += (const VectorBase &v)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+
+    const int ierr = vector->Update (1.0, *(v.vector), 1.0);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  inline
+  VectorBase &
+  VectorBase::operator -= (const VectorBase &v)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+
+    const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  inline
+  void
+  VectorBase::add (const TrilinosScalar s)
+  {
+
+    Assert (numbers::is_finite(s),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    unsigned int n_local = local_size();
+    for (unsigned int i=0; i<n_local; i++)
+      (*vector)[0][i] += s;
+  }
+
+
+
+  inline
+  void
+  VectorBase::add (const TrilinosScalar  a,
+                  const VectorBase     &v)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->Update(a, *(v.vector), 1.);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  VectorBase::add (const TrilinosScalar  a,
+                  const VectorBase     &v,
+                  const TrilinosScalar  b,
+                  const VectorBase     &w)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+    Assert (size() == w.size(),
+           ExcDimensionMismatch(size(), w.size()));
+
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(b),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  VectorBase::sadd (const TrilinosScalar  s,
+                   const VectorBase     &v)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+
+    Assert (numbers::is_finite(s),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->Update(1., *(v.vector), s);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  VectorBase::sadd (const TrilinosScalar  s,
+                   const TrilinosScalar  a,
+                   const VectorBase     &v)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+
+    Assert (numbers::is_finite(s),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->Update(a, *(v.vector), s);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  VectorBase::sadd (const TrilinosScalar  s,
+                   const TrilinosScalar  a,
+                   const VectorBase     &v,
+                   const TrilinosScalar  b,
+                   const VectorBase     &w)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+    Assert (size() == w.size(),
+           ExcDimensionMismatch(size(), w.size()));
+
+
+    Assert (numbers::is_finite(s),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(b),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  VectorBase::sadd (const TrilinosScalar  s,
+                   const TrilinosScalar  a,
+                   const VectorBase     &v,
+                   const TrilinosScalar  b,
+                   const VectorBase     &w,
+                   const TrilinosScalar  c,
+                   const VectorBase     &x)
+  {
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+    Assert (size() == w.size(),
+           ExcDimensionMismatch(size(), w.size()));
+    Assert (size() == x.size(),
+           ExcDimensionMismatch(size(), x.size()));
+
+    Assert (numbers::is_finite(s),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(b),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(c),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+                                        // Update member can only
+                                       // input two other vectors so
+                                       // do it in two steps
+    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    const int jerr = vector->Update(c, *(x.vector), 1.);
+    Assert (jerr == 0, ExcTrilinosError(jerr));
+  }
+
+
+
+  inline
+  void
+  VectorBase::scale (const VectorBase &factors)
+  {
+    Assert (size() == factors.size(),
+           ExcDimensionMismatch(size(), factors.size()));
+
+    const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  inline
+  void
+  VectorBase::equ (const TrilinosScalar  a,
+                  const VectorBase     &v)
+  {
+
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+                                  // If we don't have the same map, copy.
+    if (local_range() != v.local_range())
+      {
+       vector.reset();
+       last_action = Zero;
+       *vector = *v.vector;
+       *this *= a;
+      }
+    else
+      {
+                                  // Otherwise, just update
+       Assert (vector->Map().SameAs(v.vector->Map()) == true,
+               ExcMessage ("The Epetra maps in the assignment operator ="
+                           " do not match, even though the local_range "
+                           " seems to be the same. Check vector setup!"));
+       int ierr;
+       ierr = vector->GlobalAssemble(last_action);
+       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+       ierr = vector->Update(a, *v.vector, 0.0);
+       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+       last_action = Zero;
+      }
+
+  }
+
+
+
+  inline
+  void
+  VectorBase::equ (const TrilinosScalar  a,
+                  const VectorBase     &v,
+                  const TrilinosScalar  b,
+                  const VectorBase     &w)
+  {
+
+    Assert (v.size() == w.size(),
+           ExcDimensionMismatch (v.size(), w.size()));
+
+    Assert (numbers::is_finite(a),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+    Assert (numbers::is_finite(b),
+           ExcMessage("The given value is not finite but "
+                      "either infinite or Not A Number (NaN)"));
+
+                                  // If we don't have the same map, copy.
+    if (local_range() != v.local_range())
+      {
+       vector.reset();
+       last_action = Zero;
+       *vector = *v.vector;
+       sadd (a, b, w);
+      }
+    else
+      {
+                                  // Otherwise, just update
+       Assert (vector->Map().SameAs(v.vector->Map()) == true,
+               ExcMessage ("The Epetra maps in the assignment operator ="
+                           " do not match, even though the local_range "
+                           " seems to be the same. Check vector setup!"));
+       int ierr;
+       ierr = vector->GlobalAssemble(last_action);
+       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+       ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
+       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+       last_action = Zero;
+      }
+  }
+
+
+
+  inline
+  void
+  VectorBase::ratio (const VectorBase &v,
+                    const VectorBase &w)
+  {
+    Assert (v.size() == w.size(),
+           ExcDimensionMismatch (v.size(), w.size()));
+
+    Assert (size() == w.size(),
+           ExcDimensionMismatch (size(), w.size()));
+
+    const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector), 
+                                               *(v.vector), 0.0);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
   inline
   const Epetra_MultiVector &
   VectorBase::trilinos_vector () const
index 4ac73a823d148c90ad7af1b2887461d2832e49b5..e3ce75b68bf1d113a6e9a719ed88befd111be715 100755 (executable)
@@ -705,40 +705,6 @@ namespace TrilinosWrappers
 
 
 
-  void
-  SparseMatrix::compress ()
-  {
-                                 // flush buffers
-    int ierr;
-    ierr = matrix->GlobalAssemble (col_map, row_map, true);
-    
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    ierr = matrix->OptimizeStorage ();
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    last_action = Zero;
-
-    compressed = true;
-  }
-
-
-
-  SparseMatrix &
-  SparseMatrix::operator = (const double d)
-  {
-    Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
-
-    compress ();
-
-    const int ierr = matrix->PutScalar(d);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
   void
   SparseMatrix::clear_row (const unsigned int   row,
                           const TrilinosScalar new_diag_value)
@@ -958,57 +924,6 @@ namespace TrilinosWrappers
 
 
 
-  unsigned int
-  SparseMatrix::m () const
-  {
-    int n_rows = matrix -> NumGlobalRows();
-
-    return n_rows;
-  }
-
-
-
-  unsigned int
-  SparseMatrix::n () const
-  {
-    unsigned int n_cols = matrix -> NumGlobalCols();
-    return n_cols;
-  }
-
-
-
-  unsigned int
-  SparseMatrix::local_size () const
-  {
-    unsigned int n_rows = matrix -> NumMyRows();
-
-    return n_rows;
-  }
-
-
-
-  std::pair<unsigned int, unsigned int>
-  SparseMatrix::local_range () const
-  {
-    unsigned int begin, end;
-    begin = matrix -> RowMap().MinMyGID();
-    end = matrix -> RowMap().MaxMyGID()+1;
-    
-    return std::make_pair (begin, end);
-  }
-
-
-
-  unsigned int
-  SparseMatrix::n_nonzero_elements () const
-  {
-    unsigned int nnz = matrix->NumGlobalNonzeros();
-
-    return nnz;
-  }
-
-
-
   unsigned int
   SparseMatrix::row_length (const unsigned int row) const
   {
@@ -1033,191 +948,6 @@ namespace TrilinosWrappers
 
 
 
-  TrilinosScalar
-  SparseMatrix::l1_norm () const
-  {
-    if (matrix->Filled() == false)
-      matrix->GlobalAssemble(col_map, row_map, true);
-
-    TrilinosScalar result = matrix->NormOne();
-
-    return result;
-  }
-  
-  
-
-  TrilinosScalar
-  SparseMatrix::linfty_norm () const
-  {
-    if (matrix->Filled() == false)
-      matrix->GlobalAssemble(col_map, row_map, true);
-
-    TrilinosScalar result = matrix->NormInf();
-
-    return result;
-  }
-
-
-
-  TrilinosScalar
-  SparseMatrix::frobenius_norm () const
-  {
-    if (matrix->Filled() == false)
-      matrix->GlobalAssemble(col_map, row_map, true);
-
-    TrilinosScalar result = matrix->NormFrobenius();
-
-    return result;
-  }
-
-
-
-  SparseMatrix &
-  SparseMatrix::operator *= (const TrilinosScalar a)
-  {
-    const int ierr = matrix->Scale (a);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
-  SparseMatrix &
-  SparseMatrix::operator /= (const TrilinosScalar a)
-  {
-    Assert (a !=0, ExcDivideByZero());
-
-    const TrilinosScalar factor = 1./a;
-
-    const int ierr = matrix->Scale (factor);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
-  void
-  SparseMatrix::vmult (VectorBase       &dst,
-                      const VectorBase &src) const
-  {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-
-    if (matrix->Filled() == false)
-      matrix->GlobalAssemble(col_map, row_map, true);
-
-    Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true,
-           ExcMessage ("Column map of matrix does not fit with vector map!"));
-    Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true,
-           ExcMessage ("Row map of matrix does not fit with vector map!"));
-
-    const int ierr = matrix->Multiply (false, src.trilinos_vector(), 
-                                      dst.trilinos_vector());
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  SparseMatrix::Tvmult (VectorBase       &dst,
-                       const VectorBase &src) const
-  {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-
-    if (matrix->Filled() == false)
-      matrix->GlobalAssemble(col_map, row_map, true);
-
-    Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true,
-           ExcMessage ("Column map of matrix does not fit with vector map!"));
-    Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true,
-           ExcMessage ("Row map of matrix does not fit with vector map!"));
-
-    const int ierr = matrix->Multiply (true, src.trilinos_vector(), 
-                                      dst.trilinos_vector());
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  SparseMatrix::vmult_add (VectorBase       &dst,
-                          const VectorBase &src) const
-  {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-
-                                  // Choose to reinit the vector with fast
-                                  // argument set, which does not overwrite
-                                  // the content -- this is what we need
-                                  // since we're going to overwrite that
-                                  // anyway in the vmult operation.
-    temp_vector.reinit(dst, true);
-
-    vmult (temp_vector, src);
-    dst += temp_vector;
-  }
-
-
-
-  void
-  SparseMatrix::Tvmult_add (VectorBase       &dst,
-                           const VectorBase &src) const
-  {
-    Assert (&src != &dst, ExcSourceEqualsDestination());
-
-    temp_vector.reinit(dst, true);
-
-    vmult (temp_vector, src);
-    dst += temp_vector;
-  }
-
-
-
-  TrilinosScalar
-  SparseMatrix::matrix_norm_square (const VectorBase &v) const
-  {
-    Assert (row_map.SameAs(col_map),
-           ExcDimensionMismatch(row_map.NumGlobalElements(),
-                                col_map.NumGlobalElements()));
-
-    temp_vector.reinit(v);
-
-    vmult (temp_vector, v);
-    return temp_vector*v;
-  }
-
-
-
-  TrilinosScalar
-  SparseMatrix::matrix_scalar_product (const VectorBase &u,
-                                      const VectorBase &v) const
-  {
-    Assert (row_map.SameAs(col_map),
-           ExcDimensionMismatch(row_map.NumGlobalElements(),
-                                col_map.NumGlobalElements()));
-
-    temp_vector.reinit(v);
-
-    vmult (temp_vector, v);
-    return u*temp_vector;
-  }
-
-
-
-  TrilinosScalar
-  SparseMatrix::residual (VectorBase       &dst,
-                         const VectorBase &x,
-                         const VectorBase &b) const
-  {
-    vmult (dst, x);
-    dst -= b;
-    dst *= -1.;
-
-    return dst.l2_norm();
-  }
-
-
-
   void
   SparseMatrix::mmult (SparseMatrix       &C,
                       const SparseMatrix &B,
@@ -1625,15 +1355,14 @@ namespace TrilinosWrappers
   
 
 
-                                 // TODO: Currently this only flips a
-                                 // flag that tells Trilinos that any
-                                 // application should be done with
-                                 // the transpose. However, the
-                                 // matrix structure is not
-                                 // reset. Can we leave it like this?
   void
   SparseMatrix::transpose () 
   {
+                                 // This only flips a flag that tells
+                                 // Trilinos that any vmult operation
+                                 // should be done with the
+                                 // transpose. However, the matrix
+                                 // structure is not reset.
     int ierr;
 
     if (!matrix->UseTranspose())
index 9039cbb9b2bf442842a658938041e67c0ff5050d..8fea796f7250ce99cb315c44e834ed5bf71ccad0 100644 (file)
@@ -106,52 +106,6 @@ namespace TrilinosWrappers
 
 
 
-  void
-  VectorBase::reinit (const VectorBase &v,
-                     const bool        fast)
-  {
-    Assert (&*vector != 0,
-           ExcMessage("Vector has not been constructed properly."));
-
-    if (fast == false || local_range() != v.local_range())
-      vector = std::auto_ptr<Epetra_FEVector>(new Epetra_FEVector(*v.vector));
-  }
-
-
-
-  void
-  VectorBase::compress ()
-  {
-                                // Now pass over the information about
-                                // what we did last to the vector.
-    const int ierr = vector->GlobalAssemble(last_action);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-    last_action = Zero;
-  
-    compressed = true;
-  }
-
-
-
-  VectorBase &
-  VectorBase::operator = (const TrilinosScalar s)
-  {
-
-    Assert (numbers::is_finite(s),
-           ExcMessage("The given value is not finite but either "
-                      "infinite or Not A Number (NaN)"));
-
-    compress();
-
-    const int ierr = vector->PutScalar(s);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
   VectorBase &
   VectorBase::operator = (const VectorBase &v)
   {
@@ -202,62 +156,11 @@ namespace TrilinosWrappers
                                     // particularly frequent operation
     std::pair<unsigned int, unsigned int>
       local_range = this->local_range ();
-    std::vector<unsigned int>   indices (local_range.second -
-                                        local_range.first);
-    std::vector<TrilinosScalar> values (local_range.second -
-                                       local_range.first);
+    for (unsigned int i=local_range.first; i<local_range.second; ++i)
+      (*vector)[0][i-local_range.first] = v(i);
 
-    for (unsigned int i=0; i<local_range.second-local_range.first; ++i)
-      {
-       indices[i] = i + local_range.first;
-       values[i] = v(i + local_range.first);
-      }
-    
-    set (indices.size(), &indices[0], &values[0]);
-      
     return *this;
   }
-  
-
-
-  bool
-  VectorBase::operator == (const VectorBase &v) const
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-    unsigned int i;
-    for (i=0; i<size(); i++) 
-      if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
-
-    return true;
-  }
-
-
-
-  bool
-  VectorBase::operator != (const VectorBase &v) const
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-
-    return (!(*this==v));
-  }
-
-
-
-  unsigned int
-  VectorBase::size () const
-  {
-    return (unsigned int) vector->Map().NumGlobalElements();
-  }
-
-
-
-  unsigned int
-  VectorBase::local_size () const
-  {
-    return (unsigned int) vector->Map().NumMyElements();
-  }
 
 
 
@@ -315,110 +218,53 @@ namespace TrilinosWrappers
 
 
 
-  TrilinosScalar
-  VectorBase::operator * (const VectorBase &vec) const
-  {
-    Assert (local_range() == vec.local_range(),
-           ExcDimensionMismatch(size(), vec.size()));
-
-    TrilinosScalar result;
-
-    const int ierr = vector->Dot(*(vec.vector), &result);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return result;
-  }
-
-
-
-  VectorBase::real_type
-  VectorBase::norm_sqr () const
-  {
-    const TrilinosScalar d = l2_norm();
-    return d*d;
-  }
-
-
-
-  TrilinosScalar
-  VectorBase::mean_value () const
-  {
-    TrilinosScalar mean;
-
-    const int ierr = vector->MeanValue (&mean);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return mean;
-  }
-
-
-
-  VectorBase::real_type
-  VectorBase::l1_norm () const
+  inline
+  void
+  VectorBase::add (const VectorBase &v,
+                  const bool        allow_different_maps)
   {
-    TrilinosScalar d;
-
-    const int ierr = vector->Norm1 (&d);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return d;
-  }
-
-
+    if (allow_different_maps == false)
+      *this += v;
+    else
+      {
+       AssertThrow (size() == v.size(),
+                    ExcDimensionMismatch (size(), v.size()));
 
-  VectorBase::real_type
-  VectorBase::l2_norm () const
-  {
-    TrilinosScalar d;
+       Epetra_Import data_exchange (vector->Map(), v.vector->Map());
 
-    const int ierr = vector->Norm2 (&d);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+       int ierr = vector->Import(*v.vector, data_exchange, Add);
+       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-    return d;
+       last_action = Insert;
+      }
   }
 
 
 
-  VectorBase::real_type
-  VectorBase::lp_norm (const TrilinosScalar p) const
+  bool
+  VectorBase::operator == (const VectorBase &v) const
   {
-                                        // get a representation of the
-                                        // vector and loop over all
-                                        // the elements
-    TrilinosScalar *start_ptr;
-    int leading_dimension;
-    int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    TrilinosScalar norm = 0;
-    TrilinosScalar sum=0;
-
-    const TrilinosScalar * ptr  = start_ptr;
-
-                                       // add up elements 
-                                      // TODO: This
-                                      // won't work in parallel like
-                                      // this. Find out a better way to
-                                      // this in that case.
-    while (ptr != start_ptr+size())
-      sum += std::pow(std::fabs(*ptr++), p);
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
+    if (local_size() != v.local_size())
+      return false;
 
-    norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
+    unsigned int i;
+    for (i=0; i<local_size(); i++) 
+      if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
 
-    return norm;
+    return true;
   }
 
 
 
-  VectorBase::real_type
-  VectorBase::linfty_norm () const
+  bool
+  VectorBase::operator != (const VectorBase &v) const
   {
-    TrilinosScalar d;
-
-    const int ierr = vector->NormInf (&d);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+    Assert (size() == v.size(),
+           ExcDimensionMismatch(size(), v.size()));
 
-    return d;
+    return (!(*this==v));
   }
 
 
@@ -487,380 +333,6 @@ namespace TrilinosWrappers
 
 
 
-  VectorBase &
-  VectorBase::operator *= (const TrilinosScalar a)
-  {
-
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const int ierr = vector->Scale(a);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
-  VectorBase &
-  VectorBase::operator /= (const TrilinosScalar a)
-  {
-
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const TrilinosScalar factor = 1./a;
-
-    Assert (numbers::is_finite(factor),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const int ierr = vector->Scale(factor);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
-  VectorBase &
-  VectorBase::operator += (const VectorBase &v)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-
-    const int ierr = vector->Update (1.0, *(v.vector), 1.0);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
-  VectorBase &
-  VectorBase::operator -= (const VectorBase &v)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-
-    const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    return *this;
-  }
-
-
-
-  void
-  VectorBase::add (const TrilinosScalar s)
-  {
-
-    Assert (numbers::is_finite(s),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    unsigned int n_local = local_size();
-    int ierr;
-
-    for (unsigned int i=0; i<n_local; i++)
-      {
-       ierr = vector->SumIntoMyValue(i,0,s);
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-      }
-  }
-
-
-
-  void
-  VectorBase::add (const VectorBase &v,
-                  const bool        allow_different_maps)
-  {
-    if (allow_different_maps == false)
-      *this += v;
-    else
-      {
-       AssertThrow (size() == v.size(),
-                    ExcDimensionMismatch (size(), v.size()));
-
-       Epetra_Import data_exchange (vector->Map(), v.vector->Map());
-
-       int ierr = vector->Import(*v.vector, data_exchange, Add);
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-       last_action = Insert;
-      }
-  }
-
-
-
-  void
-  VectorBase::add (const TrilinosScalar  a,
-                  const VectorBase     &v)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const int ierr = vector->Update(a, *(v.vector), 1.);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  VectorBase::add (const TrilinosScalar  a,
-                  const VectorBase     &v,
-                  const TrilinosScalar  b,
-                  const VectorBase     &w)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-    Assert (size() == w.size(),
-           ExcDimensionMismatch(size(), w.size()));
-
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  VectorBase::sadd (const TrilinosScalar  s,
-                   const VectorBase     &v)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-
-    Assert (numbers::is_finite(s),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const int ierr = vector->Update(1., *(v.vector), s);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  VectorBase::sadd (const TrilinosScalar  s,
-                   const TrilinosScalar  a,
-                   const VectorBase     &v)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-
-    Assert (numbers::is_finite(s),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const int ierr = vector->Update(a, *(v.vector), s);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  VectorBase::sadd (const TrilinosScalar  s,
-                   const TrilinosScalar  a,
-                   const VectorBase     &v,
-                   const TrilinosScalar  b,
-                   const VectorBase     &w)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-    Assert (size() == w.size(),
-           ExcDimensionMismatch(size(), w.size()));
-
-
-    Assert (numbers::is_finite(s),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  VectorBase::sadd (const TrilinosScalar  s,
-                   const TrilinosScalar  a,
-                   const VectorBase     &v,
-                   const TrilinosScalar  b,
-                   const VectorBase     &w,
-                   const TrilinosScalar  c,
-                   const VectorBase     &x)
-  {
-    Assert (size() == v.size(),
-           ExcDimensionMismatch(size(), v.size()));
-    Assert (size() == w.size(),
-           ExcDimensionMismatch(size(), w.size()));
-    Assert (size() == x.size(),
-           ExcDimensionMismatch(size(), x.size()));
-
-    Assert (numbers::is_finite(s),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(c),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-                                        // Update member can only
-                                       // input two other vectors so
-                                       // do it in two steps
-    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-    const int jerr = vector->Update(c, *(x.vector), 1.);
-    Assert (jerr == 0, ExcTrilinosError(jerr));
-  }
-
-
-
-  void
-  VectorBase::scale (const VectorBase &factors)
-  {
-    Assert (size() == factors.size(),
-           ExcDimensionMismatch(size(), factors.size()));
-
-    const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-
-  void
-  VectorBase::equ (const TrilinosScalar  a,
-                  const VectorBase     &v)
-  {
-
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-                                  // If we don't have the same map, copy.
-    if (local_range() != v.local_range())
-      {
-       vector.reset();
-       last_action = Zero;
-       *vector = *v.vector;
-       *this *= a;
-      }
-    else
-      {
-                                  // Otherwise, just update
-       Assert (vector->Map().SameAs(v.vector->Map()) == true,
-               ExcMessage ("The Epetra maps in the assignment operator ="
-                           " do not match, even though the local_range "
-                           " seems to be the same. Check vector setup!"));
-       int ierr;
-       ierr = vector->GlobalAssemble(last_action);
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-       ierr = vector->Update(a, *v.vector, 0.0);
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-       last_action = Zero;
-      }
-
-  }
-
-
-
-  void
-  VectorBase::equ (const TrilinosScalar  a,
-                  const VectorBase     &v,
-                  const TrilinosScalar  b,
-                  const VectorBase     &w)
-  {
-
-    Assert (v.size() == w.size(),
-           ExcDimensionMismatch (v.size(), w.size()));
-
-    Assert (numbers::is_finite(a),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b),
-           ExcMessage("The given value is not finite but "
-                      "either infinite or Not A Number (NaN)"));
-
-                                  // If we don't have the same map, copy.
-    if (local_range() != v.local_range())
-      {
-       vector.reset();
-       last_action = Zero;
-       *vector = *v.vector;
-       sadd (a, b, w);
-      }
-    else
-      {
-                                  // Otherwise, just update
-       Assert (vector->Map().SameAs(v.vector->Map()) == true,
-               ExcMessage ("The Epetra maps in the assignment operator ="
-                           " do not match, even though the local_range "
-                           " seems to be the same. Check vector setup!"));
-       int ierr;
-       ierr = vector->GlobalAssemble(last_action);
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-       ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
-       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
-       last_action = Zero;
-      }
-  }
-
-
-
-  void
-  VectorBase::ratio (const VectorBase &v,
-                    const VectorBase &w)
-  {
-    Assert (v.size() == w.size(),
-           ExcDimensionMismatch (v.size(), w.size()));
-
-    Assert (size() == w.size(),
-           ExcDimensionMismatch (size(), w.size()));
-
-    const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector), 
-                                               *(v.vector), 0.0);
-
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-  }
-
-
-  
                                        // TODO: up to now only local
                                         // data printed out! Find a
                                         // way to neatly output

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.