"Fossies" - the Fresh Open Source Software Archive

Member "dune-istl-2.7.1/dune/istl/test/scalarproductstest.cc" (26 Nov 2020, 4861 Bytes) of package /linux/misc/dune/dune-istl-2.7.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
    2 // vi: set et ts=4 sw=2 sts=2:
    3 #include "config.h"
    4 
    5 #include <dune/common/fvector.hh>
    6 #include <dune/common/parallel/mpihelper.hh>
    7 #include <dune/common/test/testsuite.hh>
    8 #include <dune/common/scalarvectorview.hh>
    9 
   10 #include <dune/istl/bvector.hh>
   11 #include <dune/istl/owneroverlapcopy.hh>
   12 #include <dune/istl/scalarproducts.hh>
   13 
   14 using namespace Dune;
   15 
   16 // scalar ordering doesn't work for complex numbers
   17 template <class ScalarProduct, class BlockVector>
   18 TestSuite scalarProductTest(const ScalarProduct& scalarProduct,
   19                             const size_t numBlocks)
   20 {
   21   TestSuite t;
   22 
   23   static_assert(std::is_same<BlockVector, typename ScalarProduct::domain_type>::value,
   24                 "ScalarProduct does not properly export its vector type!");
   25 
   26   typedef typename ScalarProduct::field_type field_type;
   27   typedef typename ScalarProduct::real_type  real_type;
   28   const real_type myEps((real_type)1e-6);
   29 
   30   static_assert(std::is_same<typename FieldTraits<field_type>::real_type, real_type>::value,
   31                 "real_type does not match field_type");
   32 
   33   typedef typename BlockVector::size_type size_type;
   34 
   35   // empty vectors
   36   BlockVector one(numBlocks);
   37 
   38   const size_type blockSize = Impl::asVector(one[0]).size();
   39 
   40   t.require(numBlocks==one.N());
   41   const size_type length = numBlocks * blockSize; // requires inner block size of VariableBlockVector to be 1!
   42 
   43   std::cout << __func__ << "\t \t ( " << className(one) << " )" << std::endl << std::endl;
   44 
   45   // initialize test vector with data
   46   for(size_type i=0; i < numBlocks; ++i)
   47     one[i] = 1.;
   48 
   49   field_type result = field_type();
   50 
   51   // global operator * tests
   52   result = scalarProduct.dot(one,one);
   53 
   54   // The Euclidean norm of a vector with all entries set to '1' equals its number of entries
   55   t.check(std::abs(result-field_type(length))<= myEps);
   56 
   57   auto sp   = scalarProduct.dot(one,one);
   58   auto norm = scalarProduct.norm(one);
   59 
   60   t.check(std::abs(sp - norm*norm) <=myEps);
   61 
   62   return t;
   63 }
   64 
   65 
   66 int main(int argc, char** argv)
   67 {
   68   MPIHelper::instance(argc, argv);
   69 
   70   TestSuite t;
   71   const size_t BlockSize = 5;
   72   const size_t numBlocks = 10;
   73 
   74   // Test the ScalarProduct class
   75   {
   76     using Vector = BlockVector<FieldVector<double,BlockSize> >;
   77     using ScalarProduct = ScalarProduct<Vector>;
   78     ScalarProduct scalarProduct;
   79     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
   80   }
   81 
   82   {
   83     using Vector = BlockVector<float>;
   84     using ScalarProduct = ScalarProduct<Vector>;
   85     ScalarProduct scalarProduct;
   86     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
   87   }
   88 
   89   {
   90     using Vector = BlockVector<FieldVector<std::complex<float>, 1> >;
   91     using ScalarProduct = ScalarProduct<Vector>;
   92     ScalarProduct scalarProduct;
   93     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
   94   }
   95 
   96   // Test the SeqScalarProduct class
   97   {
   98     using Vector = BlockVector<FieldVector<double,BlockSize> >;
   99     using ScalarProduct = SeqScalarProduct<Vector>;
  100     ScalarProduct scalarProduct;
  101     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
  102   }
  103 
  104   {
  105     using Vector = BlockVector<float>;
  106     using ScalarProduct = SeqScalarProduct<Vector>;
  107     ScalarProduct scalarProduct;
  108     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
  109   }
  110 
  111   {
  112     using Vector = BlockVector<FieldVector<std::complex<float>, 1> >;
  113     using ScalarProduct = SeqScalarProduct<Vector>;
  114     ScalarProduct scalarProduct;
  115     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
  116   }
  117 
  118 #if HAVE_MPI
  119   // Test the ParallelScalarProduct class
  120   {
  121     using Vector = BlockVector<FieldVector<double,BlockSize> >;
  122     using ScalarProduct = ParallelScalarProduct<Vector, OwnerOverlapCopyCommunication<std::size_t,std::size_t> >;
  123     OwnerOverlapCopyCommunication<std::size_t,std::size_t> communicator;
  124     ScalarProduct scalarProduct(communicator,SolverCategory::nonoverlapping);
  125     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
  126   }
  127 
  128   {
  129     using Vector = BlockVector<float>;
  130     using ScalarProduct = ParallelScalarProduct<Vector, OwnerOverlapCopyCommunication<std::size_t,std::size_t> >;
  131     OwnerOverlapCopyCommunication<std::size_t,std::size_t> communicator;
  132     ScalarProduct scalarProduct(communicator,SolverCategory::nonoverlapping);
  133     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
  134   }
  135 
  136   {
  137     using Vector = BlockVector<FieldVector<std::complex<double>, 1> >;
  138     using ScalarProduct = ParallelScalarProduct<Vector, OwnerOverlapCopyCommunication<std::size_t,std::size_t> >;
  139     OwnerOverlapCopyCommunication<std::size_t,std::size_t> communicator;
  140     ScalarProduct scalarProduct(communicator,SolverCategory::nonoverlapping);
  141     scalarProductTest<ScalarProduct, Vector>(scalarProduct,numBlocks);
  142   }
  143 #endif
  144 
  145   return t.exit();
  146 }