diff --git a/CMakeLists.txt b/CMakeLists.txt
index a69df81fe7b87b04ce379ef58791f3ae95fcbebd..50126789283f0173207554466058f6e6752db9ad 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,7 +33,7 @@ if( WITH_TEMPLATE_EXPLICIT_INSTANTIATION STREQUAL "yes" )
endif()
if( WITH_CUDA STREQUAL "yes" )
- AddCompilerFlag( "-DHAVE_NOT_CXX11" )
+ AddCompilerFlag( "-DHAVE_NOT_CXX11 -U_GLIBCXX_ATOMIC_BUILTINS -U_GLIBCXX_USE_INT128" )
else()
AddCompilerFlag( "-std=gnu++0x" )
endif()
diff --git a/install b/install
index 2b9489d97e9e589b113b65761c9ff77ccfe4ff01..029a1902e7722ed7efbba631ff47c2e9c92c7784 100755
--- a/install
+++ b/install
@@ -2,7 +2,7 @@
TARGET=TNL
INSTALL_PREFIX=${HOME}/local
-WITH_CUDA=yes
+WITH_CUDA=no
WITH_CUSPARSE=no
CUDA_ARCHITECTURE=2.0
TEMPLATE_EXPLICIT_INSTANTIATION=yes
diff --git a/src/core/vectors/CMakeLists.txt b/src/core/vectors/CMakeLists.txt
index d41998056a62784e2bed082fc24cb69d8ff7ca8a..f90d399ba6d40d17a45b54c06dff171eca32f1c2 100755
--- a/src/core/vectors/CMakeLists.txt
+++ b/src/core/vectors/CMakeLists.txt
@@ -1,5 +1,6 @@
set( headers tnlVector.h
tnlMultiVector.h
- tnlSharedVector.h )
+ tnlSharedVector.h
+ tnlVectorOperations.h )
INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/core/vectors )
\ No newline at end of file
diff --git a/src/core/vectors/tnlVectorOperations.h b/src/core/vectors/tnlVectorOperations.h
new file mode 100644
index 0000000000000000000000000000000000000000..5d6031156edff2d67065c4a3cee0d1bcd2f1b4a2
--- /dev/null
+++ b/src/core/vectors/tnlVectorOperations.h
@@ -0,0 +1,207 @@
+/***************************************************************************
+ tnlVectorOperations.h - description
+ -------------------
+ begin : Nov 8, 2012
+ copyright : (C) 2012 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLVECTOROPERATIONS_H_
+#define TNLVECTOROPERATIONS_H_
+
+#include <core/cuda/device-check.h>
+#include <core/cuda/cuda-reduction.h>
+#include <core/cuda/reduction-operations.h>
+#include <core/tnlHost.h>
+#include <core/tnlCuda.h>
+
+template< typename Device >
+class tnlVectorOperations{};
+
+template<>
+class tnlVectorOperations< tnlHost >
+{
+ public:
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorMax( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorMin( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorAbsMax( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorAbsMin( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorLpNorm( const Vector& v,
+ const typename Vector :: RealType& p );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorSum( const Vector& v );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceMax( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceMin( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceAbsMax( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceAbsMin( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceLpNorm( const Vector1& v1,
+ const Vector2& v2,
+ const typename Vector1 :: RealType& p );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceSum( const Vector1& v1,
+ const Vector2& v2 );
+ template< typename Vector >
+ static void vectorScalarMultiplication( Vector& v,
+ const typename Vector :: RealType& alpha );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorSdot( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha );
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxmy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha );
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpsby( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const typename Vector1 :: RealType& beta );
+
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpsbz( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta );
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpsbzpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta );
+};
+
+template<>
+class tnlVectorOperations< tnlCuda >
+{
+ public:
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorMax( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorMin( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorAbsMax( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorAbsMin( const Vector& v );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorLpNorm( const Vector& v,
+ const typename Vector :: RealType& p );
+
+ template< typename Vector >
+ static typename Vector :: RealType getVectorSum( const Vector& v );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceMax( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceMin( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceAbsMax( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceAbsMin( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceLpNorm( const Vector1& v1,
+ const Vector2& v2,
+ const typename Vector1 :: RealType& p );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorDifferenceSum( const Vector1& v1,
+ const Vector2& v2 );
+ template< typename Vector >
+ static void vectorScalarMultiplication( Vector& v,
+ const typename Vector :: RealType& alpha );
+
+ template< typename Vector1, typename Vector2 >
+ static typename Vector1 :: RealType getVectorSdot( const Vector1& v1,
+ const Vector2& v2 );
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha );
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxmy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha );
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpsby( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const typename Vector1 :: RealType& beta );
+
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpsbz( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta );
+
+ template< typename Vector1, typename Vector2 >
+ static void vectorSaxpsbzpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta );
+};
+
+#include <implementation/core/vectors/tnlVectorOperationsHost_impl.h>
+#include <implementation/core/vectors/tnlVectorOperationsCuda_impl.h>
+
+#endif /* TNLVECTOROPERATIONS_H_ */
diff --git a/src/implementation/core/CMakeLists.txt b/src/implementation/core/CMakeLists.txt
index 2dbef233e035191e848e3acea417562883f0c95f..844e450e66966b536cd67f6cf2c08e8987fbcfc4 100755
--- a/src/implementation/core/CMakeLists.txt
+++ b/src/implementation/core/CMakeLists.txt
@@ -19,7 +19,8 @@ set( common_SOURCES
${CURRENT_DIR}/tnlTimerCPU.cpp
${CURRENT_DIR}/mfilename.cpp
${CURRENT_DIR}/mpi-supp.cpp
- ${CURRENT_DIR}/tnlHost_impl.cpp )
+ ${CURRENT_DIR}/tnlHost_impl.cpp
+ ${CURRENT_DIR}/tnlCuda.cpp )
IF( BUILD_CUDA )
set( tnl_implementation_core_CUDA__SOURCES
@@ -27,8 +28,7 @@ IF( BUILD_CUDA )
${tnl_implementation_core_cuda_CUDA__SOURCES}
${tnl_implementation_core_vectors_CUDA__SOURCES}
${common_SOURCES}
- ${CURRENT_DIR}/memory-operations_impl.cu
- ${CURRENT_DIR}/tnlCuda.cu
+ ${CURRENT_DIR}/memory-operations_impl.cu
PARENT_SCOPE )
ENDIF()
diff --git a/src/implementation/core/memory-operations_impl.cpp b/src/implementation/core/memory-operations_impl.cpp
index e0bd11517cc14ae3ccae5d0bf0d942e817c90f2f..2038056cca69454582348f2d1857f0350716319d 100644
--- a/src/implementation/core/memory-operations_impl.cpp
+++ b/src/implementation/core/memory-operations_impl.cpp
@@ -75,19 +75,19 @@ template bool setMemoryHost( float* data, const float& value, const long int siz
template bool setMemoryHost( double* data, const double& value, const long int size );
template bool setMemoryHost( long double* data, const long double& value, const long int size );
-template bool setMemoryCuda( char* data, const char& value, const int size );
-template bool setMemoryCuda( int* data, const int& value, const int size );
-template bool setMemoryCuda( long int* data, const long int& value, const int size );
-template bool setMemoryCuda( float* data, const float& value, const int size );
-template bool setMemoryCuda( double* data, const double& value, const int size );
-template bool setMemoryCuda( long double* data, const long double& value, const int size );
-
-template bool setMemoryCuda( char* data, const char& value, const long int size );
-template bool setMemoryCuda( int* data, const int& value, const long int size );
-template bool setMemoryCuda( long int* data, const long int& value, const long int size );
-template bool setMemoryCuda( float* data, const float& value, const long int size );
-template bool setMemoryCuda( double* data, const double& value, const long int size );
-template bool setMemoryCuda( long double* data, const long double& value, const long int size );
+template bool setMemoryCuda( char* data, const char& value, const int size, const int maxGridSize );
+template bool setMemoryCuda( int* data, const int& value, const int size, const int maxGridSize );
+template bool setMemoryCuda( long int* data, const long int& value, const int size, const int maxGridSize );
+template bool setMemoryCuda( float* data, const float& value, const int size, const int maxGridSize );
+template bool setMemoryCuda( double* data, const double& value, const int size, const int maxGridSize );
+template bool setMemoryCuda( long double* data, const long double& value, const int size, const int maxGridSize );
+
+template bool setMemoryCuda( char* data, const char& value, const long int size, const int maxGridSize );
+template bool setMemoryCuda( int* data, const int& value, const long int size, const int maxGridSize );
+template bool setMemoryCuda( long int* data, const long int& value, const long int size, const int maxGridSize );
+template bool setMemoryCuda( float* data, const float& value, const long int size, const int maxGridSize );
+template bool setMemoryCuda( double* data, const double& value, const long int size, const int maxGridSize );
+template bool setMemoryCuda( long double* data, const long double& value, const long int size, const int maxGridSize );
template bool copyMemoryHostToHost( char* destination, const char* source, const int size );
template bool copyMemoryHostToHost( int* destination, const int* source, const int size );
diff --git a/src/implementation/core/tnlCuda.cu b/src/implementation/core/tnlCuda.cpp
similarity index 95%
rename from src/implementation/core/tnlCuda.cu
rename to src/implementation/core/tnlCuda.cpp
index 842a583b3ff730a39268a85b45d742893e0e91e4..b5b7a195fbc130830294b1b82e4d6545b59b6cea 100644
--- a/src/implementation/core/tnlCuda.cu
+++ b/src/implementation/core/tnlCuda.cpp
@@ -1,5 +1,5 @@
/***************************************************************************
- tnlCuda.cu - description
+ tnlCuda.cpp - description
-------------------
begin : Jul 11, 2013
copyright : (C) 2013 by Tomas Oberhuber
diff --git a/src/implementation/core/vectors/CMakeLists.txt b/src/implementation/core/vectors/CMakeLists.txt
index 266d6624a3b53991e240aab2bf0a9c18a7a698fb..b63100ea0bedf536931d2822886456c15e9fde27 100755
--- a/src/implementation/core/vectors/CMakeLists.txt
+++ b/src/implementation/core/vectors/CMakeLists.txt
@@ -4,7 +4,8 @@ SET( headers tnlMultiVector1D_impl.h
tnlMultiVector4D_impl.h
tnlSharedVector_impl.h
tnlVector_impl.h
- vector-operations.h )
+ tnlVectorOperationsHost_impl.h
+ tnlVectorOperationsCuda_impl.h )
SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/core/vectors )
set( common_SOURCES
diff --git a/src/implementation/core/vectors/tnlSharedVector_impl.h b/src/implementation/core/vectors/tnlSharedVector_impl.h
index 0881b3d906a55780853f4558b6fc545733acba50..eebae6b9334ba300b405f951dbb9f301df248297 100644
--- a/src/implementation/core/vectors/tnlSharedVector_impl.h
+++ b/src/implementation/core/vectors/tnlSharedVector_impl.h
@@ -18,7 +18,7 @@
#ifndef TNLSHAREDVECTOR_H_IMPLEMENTATION
#define TNLSHAREDVECTOR_H_IMPLEMENTATION
-#include <implementation/core/vectors/vector-operations.h>
+#include <core/vectors/tnlVectorOperations.h>
template< typename Real,
typename Device,
@@ -112,7 +112,7 @@ template< typename Real,
typename Index >
Real tnlSharedVector< Real, Device, Index > :: max() const
{
- return getVectorMax( *this );
+ return tnlVectorOperations< Device > :: getVectorMax( *this );
}
template< typename Real,
@@ -120,7 +120,7 @@ template< typename Real,
typename Index >
Real tnlSharedVector< Real, Device, Index > :: min() const
{
- return getVectorMin( *this );
+ return tnlVectorOperations< Device > :: getVectorMin( *this );
}
@@ -129,7 +129,7 @@ template< typename Real,
typename Index >
Real tnlSharedVector< Real, Device, Index > :: absMax() const
{
- return getVectorAbsMax( *this );
+ return tnlVectorOperations< Device > :: getVectorAbsMax( *this );
}
template< typename Real,
@@ -137,7 +137,7 @@ template< typename Real,
typename Index >
Real tnlSharedVector< Real, Device, Index > :: absMin() const
{
- return getVectorAbsMin( *this );
+ return tnlVectorOperations< Device > :: getVectorAbsMin( *this );
}
template< typename Real,
@@ -145,7 +145,7 @@ template< typename Real,
typename Index >
Real tnlSharedVector< Real, Device, Index > :: lpNorm( const Real& p ) const
{
- return getVectorLpNorm( *this, p );
+ return tnlVectorOperations< Device > :: getVectorLpNorm( *this, p );
}
@@ -154,7 +154,7 @@ template< typename Real,
typename Index >
Real tnlSharedVector< Real, Device, Index > :: sum() const
{
- return getVectorSum( *this );
+ return tnlVectorOperations< Device > :: getVectorSum( *this );
}
@@ -164,7 +164,7 @@ template< typename Real,
template< typename Vector >
Real tnlSharedVector< Real, Device, Index > :: differenceMax( const Vector& v ) const
{
- return getVectorDifferenceMax( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceMax( *this, v );
}
@@ -174,7 +174,7 @@ template< typename Real,
template< typename Vector >
Real tnlSharedVector< Real, Device, Index > :: differenceMin( const Vector& v ) const
{
- return getVectorDifferenceMin( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceMin( *this, v );
}
@@ -184,7 +184,7 @@ template< typename Real,
template< typename Vector >
Real tnlSharedVector< Real, Device, Index > :: differenceAbsMax( const Vector& v ) const
{
- return getVectorDifferenceAbsMax( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceAbsMax( *this, v );
}
template< typename Real,
@@ -193,7 +193,7 @@ template< typename Real,
template< typename Vector >
Real tnlSharedVector< Real, Device, Index > :: differenceAbsMin( const Vector& v ) const
{
- return getVectorDifferenceAbsMin( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceAbsMin( *this, v );
}
template< typename Real,
@@ -202,7 +202,7 @@ template< typename Real,
template< typename Vector >
Real tnlSharedVector< Real, Device, Index > :: differenceLpNorm( const Vector& v, const Real& p ) const
{
- return getVectorDifferenceLpNorm( *this, v, p );
+ return tnlVectorOperations< Device > :: getVectorDifferenceLpNorm( *this, v, p );
}
@@ -212,7 +212,7 @@ template< typename Real,
template< typename Vector >
Real tnlSharedVector< Real, Device, Index > :: differenceSum( const Vector& v ) const
{
- return getVectorDifferenceSum( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceSum( *this, v );
}
@@ -221,7 +221,7 @@ template< typename Real,
typename Index >
void tnlSharedVector< Real, Device, Index > :: scalarMultiplication( const Real& alpha )
{
- vectorScalarMultiplication( *this, alpha );
+ tnlVectorOperations< Device > :: vectorScalarMultiplication( *this, alpha );
}
@@ -231,7 +231,7 @@ template< typename Real,
template< typename Vector >
Real tnlSharedVector< Real, Device, Index > :: sdot( const Vector& v )
{
- return getVectorSdot( *this, v );
+ return tnlVectorOperations< Device > :: getVectorSdot( *this, v );
}
@@ -242,7 +242,7 @@ template< typename Vector >
void tnlSharedVector< Real, Device, Index > :: saxpy( const Real& alpha,
const Vector& x )
{
- vectorSaxpy( *this, x, alpha );
+ tnlVectorOperations< Device > :: vectorSaxpy( *this, x, alpha );
}
template< typename Real,
@@ -252,7 +252,7 @@ template< typename Vector >
void tnlSharedVector< Real, Device, Index > :: saxmy( const Real& alpha,
const Vector& x )
{
- vectorSaxmy( *this, x, alpha );
+ tnlVectorOperations< Device > :: vectorSaxmy( *this, x, alpha );
}
@@ -264,7 +264,7 @@ void tnlSharedVector< Real, Device, Index > :: saxpsby( const Real& alpha,
const Vector& x,
const Real& beta )
{
- vectorSaxpsby( *this, x, alpha, beta );
+ tnlVectorOperations< Device > :: vectorSaxpsby( *this, x, alpha, beta );
}
template< typename Real,
@@ -276,7 +276,7 @@ void tnlSharedVector< Real, Device, Index > :: saxpsbz( const Real& alpha,
const Real& beta,
const Vector& z )
{
- vectorSaxpsbz( *this, x, alpha, z, beta );
+ tnlVectorOperations< Device > :: vectorSaxpsbz( *this, x, alpha, z, beta );
}
template< typename Real,
@@ -288,7 +288,7 @@ void tnlSharedVector< Real, Device, Index > :: saxpsbzpy( const Real& alpha,
const Real& beta,
const Vector& z )
{
- vectorSaxpsbz( *this, x, alpha, z, beta );
+ tnlVectorOperations< Device > :: vectorSaxpsbz( *this, x, alpha, z, beta );
}
diff --git a/src/implementation/core/vectors/tnlVectorOperationsCuda_impl.h b/src/implementation/core/vectors/tnlVectorOperationsCuda_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..e03eb31b1c02aafaf8359f4d1e4591219d4edbf7
--- /dev/null
+++ b/src/implementation/core/vectors/tnlVectorOperationsCuda_impl.h
@@ -0,0 +1,475 @@
+/***************************************************************************
+ tnlVectorOperationsCuda_impl.h - description
+ -------------------
+ begin : Nov 8, 2012
+ copyright : (C) 2012 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLVECTOROPERATIONSCUDA_IMPL_H_
+#define TNLVECTOROPERATIONSCUDA_IMPL_H_
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorMax( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+ Real result( 0 );
+ tnlParallelReductionMax< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v. getSize(),
+ v. getData(),
+ ( Real* ) 0,
+ result );
+ return result;
+}
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorMin( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionMin< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v. getSize(),
+ v. getData(),
+ ( Real* ) 0,
+ result );
+ return result;
+}
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorAbsMax( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionAbsMax< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v. getSize(),
+ v. getData(),
+ ( Real* ) 0,
+ result );
+ return result;
+}
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorAbsMin( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionAbsMin< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v. getSize(),
+ v. getData(),
+ ( Real* ) 0,
+ result );
+ return result;
+}
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorLpNorm( const Vector& v,
+ const typename Vector :: RealType& p )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+ tnlAssert( p > 0.0,
+ cerr << " p = " << p );
+
+ Real result( 0 );
+ tnlParallelReductionLpNorm< Real, Index > operation;
+ operation. setPower( p );
+ reductionOnCudaDevice( operation,
+ v. getSize(),
+ v. getData(),
+ ( Real* ) 0,
+ result );
+ return pow( result, 1.0 / p );
+}
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorSum( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionSum< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v. getSize(),
+ v. getData(),
+ ( Real* ) 0,
+ result );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorDifferenceMax( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionDiffMax< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v1. getSize(),
+ v1. getData(),
+ v2. getData(),
+ result );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorDifferenceMin( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionDiffMin< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v1. getSize(),
+ v1. getData(),
+ v2. getData(),
+ result );
+ return result;
+}
+
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorDifferenceAbsMax( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionDiffAbsMax< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v1. getSize(),
+ v1. getData(),
+ v2. getData(),
+ result );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorDifferenceAbsMin( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionDiffAbsMin< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v1. getSize(),
+ v1. getData(),
+ v2. getData(),
+ result );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorDifferenceLpNorm( const Vector1& v1,
+ const Vector2& v2,
+ const typename Vector1 :: RealType& p )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( p > 0.0,
+ cerr << " p = " << p );
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionDiffLpNorm< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v1. getSize(),
+ v1. getData(),
+ v2. getData(),
+ result );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorDifferenceSum( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result( 0 );
+ tnlParallelReductionDiffSum< Real, Index > operation;
+ reductionOnCudaDevice( operation,
+ v1. getSize(),
+ v1. getData(),
+ v2. getData(),
+ result );
+ return result;
+}
+
+
+template< typename Vector >
+void tnlVectorOperations< tnlCuda > :: vectorScalarMultiplication( Vector& v,
+ const typename Vector :: RealType& alpha )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+
+#ifdef HAVE_CUDA
+ dim3 blockSize, gridSize;
+ blockSize. x = 512;
+ gridSize. x = v. getSize() / 512 + 1;
+ // TODO: Fix this - the grid size might be limiting for large vectors.
+
+ /*
+ tnlVectorCUDAScalaMultiplicationKernel<<< gridSize, blockSize >>>( v. getSize(),
+ alpha,
+ v. getData() );*/
+#else
+ cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
+#endif
+
+}
+
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorSdot( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result( 0 );
+ /*reductionOnCudaDevice< Real,
+ Real,
+ Index,
+ tnlParallelReductionSdot >
+ ( v1. getSize(),
+ v1. getData(),
+ v2. getData(),
+ result,
+ ( Real ) 0 );*/
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlCuda > :: vectorSaxpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( y. getSize() > 0,
+ cerr << "Vector name is " << y. getName() );
+ tnlAssert( y. getSize() == x. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+#ifdef HAVE_CUDA
+ dim3 blockSize, gridSize;
+ blockSize. x = 512;
+ gridSize. x = x. getSize() / 512 + 1;
+
+ // TODO: fix this
+ /*tnlVectorCUDASaxpyKernel<<< gridSize, blockSize >>>( y. getSize(),
+ alpha,
+ x. getData(),
+ y. getData() );*/
+#else
+ cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
+#endif
+}
+
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlCuda > :: vectorSaxmy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( y. getSize() > 0,
+ cerr << "Vector name is " << y. getName() );
+ tnlAssert( y. getSize() == x. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+#ifdef HAVE_CUDA
+ dim3 blockSize, gridSize;
+ blockSize. x = 512;
+ gridSize. x = x. getSize() / 512 + 1;
+
+ // TODO: fix this
+ /*tnlVectorCUDASaxmyKernel<<< gridSize, blockSize >>>( y. getSize(),
+ alpha,
+ x. getData(),
+ y. getData() );*/
+#else
+ cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
+#endif
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlCuda > :: vectorSaxpsby( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const typename Vector1 :: RealType& beta )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( y. getSize() > 0,
+ cerr << "Vector name is " << y. getName() );
+ tnlAssert( y. getSize() == x. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+#ifdef HAVE_CUDA
+ dim3 blockSize, gridSize;
+ blockSize. x = 512;
+ gridSize. x = x. getSize() / 512 + 1;
+
+ // TODO: fix this
+ /*tnlVectorCUDASaxpsbzKernel<<< gridSize, blockSize >>>( y. getSize(),
+ alpha,
+ x. getData(),
+ beta );*/
+#else
+ cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
+#endif
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlCuda > :: vectorSaxpsbz( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( y. getSize() > 0,
+ cerr << "Vector name is " << y. getName() );
+ tnlAssert( y. getSize() == x. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+#ifdef HAVE_CUDA
+ dim3 blockSize, gridSize;
+ blockSize. x = 512;
+ gridSize. x = x. getSize() / 512 + 1;
+
+ // TODO: fix this
+ /*tnlVectorCUDASaxpsbzKernel<<< gridSize, blockSize >>>( y. getSize(),
+ alpha,
+ x. getData(),
+ beta,
+ z. getData() );*/
+#else
+ cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
+#endif
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlCuda > :: vectorSaxpsbzpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( y. getSize() > 0,
+ cerr << "Vector name is " << y. getName() );
+ tnlAssert( y. getSize() == x. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+#ifdef HAVE_CUDA
+ dim3 blockSize, gridSize;
+ blockSize. x = 512;
+ gridSize. x = x. getSize() / 512 + 1;
+
+ // TODO: fix this
+ /*tnlVectorCUDASaxpsbzpyKernel<<< gridSize, blockSize >>>( y. getSize(),
+ alpha,
+ x. getData(),
+ beta,
+ z. getData() );*/
+#else
+ cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
+#endif
+}
+
+
+#endif /* TNLVECTOROPERATIONSCUDA_IMPL_H_ */
diff --git a/src/implementation/core/vectors/tnlVectorOperationsHost_impl.h b/src/implementation/core/vectors/tnlVectorOperationsHost_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..080a46447924a2fdd08c99046c33083cded39fef
--- /dev/null
+++ b/src/implementation/core/vectors/tnlVectorOperationsHost_impl.h
@@ -0,0 +1,405 @@
+/***************************************************************************
+ tnlVectorOperationsHost_impl.h - description
+ -------------------
+ begin : Nov 8, 2012
+ copyright : (C) 2012 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLVECTOROPERATIONSHOST_IMPL_H_
+#define TNLVECTOROPERATIONSHOST_IMPL_H_
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorMax( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+ Real result = v. getElement( 0 );
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Max( result, v. getElement( i ) );
+ return result;
+}
+
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorMin( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+ Real result = v. getElement( 0 );
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Min( result, v. getElement( i ) );
+ return result;
+}
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorAbsMax( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+ Real result = fabs( v. getElement( 0 ) );
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Max( result, ( Real ) fabs( v. getElement( i ) ) );
+ return result;
+}
+
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorAbsMin( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+ Real result = fabs( v. getElement( 0 ) );
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Min( result, ( Real ) fabs( v. getElement( i ) ) );
+ return result;
+}
+
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorLpNorm( const Vector& v,
+ const typename Vector :: RealType& p )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+ tnlAssert( p > 0.0,
+ cerr << " p = " << p );
+ if( p == 1.0 )
+ {
+ Real result = fabs( v. getElement( 0 ) );
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result += fabs( v. getElement( i ) );
+ return result;
+ }
+ if( p == 2.0 )
+ {
+ Real result = v. getElement( 0 );
+ result *= result;
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ {
+ const Real aux = v. getElement( i );
+ result += aux * aux;
+ }
+ return sqrt( result );
+ }
+ Real result = pow( fabs( v. getElement( 0 ) ), p );
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result += pow( fabs( v. getElement( i ) ), p );
+ return pow( result, 1.0 / p );
+}
+
+template< typename Vector >
+typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorSum( const Vector& v )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+
+ Real result = v. getElement( 0 );
+ const Index n = v. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result += v. getElement( i );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceMax( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+ Real result = v1. getElement( 0 ) - v2. getElement( 0 );
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Max( result, v1. getElement( i ) - v2. getElement( i ) );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceMin( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result = v1. getElement( 0 ) - v2. getElement( 0 );
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Min( result, v1. getElement( i ) - v2. getElement( i ) );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceAbsMax( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Max( result, ( Real ) fabs( v1. getElement( i ) - v2. getElement( i ) ) );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceAbsMin( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result = Min( result, ( Real ) fabs( v1. getElement( i ) - v2. getElement( i ) ) );
+ return result;
+}
+
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceLpNorm( const Vector1& v1,
+ const Vector2& v2,
+ const typename Vector1 :: RealType& p )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( p > 0.0,
+ cerr << " p = " << p );
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ if( p == 1.0 )
+ {
+ Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result += fabs( v1. getElement( i ) - v2. getElement( i ) );
+ return result;
+ }
+ if( p == 2.0 )
+ {
+ Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
+ result *= result;
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ {
+ Real aux = fabs( v1. getElement( i ) - v2. getElement( i ) );
+ result += aux * aux;
+ }
+ return sqrt( result );
+ }
+ Real result = pow( fabs( v1. getElement( 0 ) - v2. getElement( 0 ) ), p );
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result += pow( fabs( v1. getElement( i ) - v2. getElement( i ) ), p );
+ return pow( result, 1.0 / p );
+}
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceSum( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result = v1. getElement( 0 ) - v2. getElement( 0 );
+ const Index n = v1. getSize();
+ for( Index i = 1; i < n; i ++ )
+ result += v1. getElement( i ) - v2. getElement( i );
+ return result;
+}
+
+
+template< typename Vector >
+void tnlVectorOperations< tnlHost > :: vectorScalarMultiplication( Vector& v,
+ const typename Vector :: RealType& alpha )
+{
+ typedef typename Vector :: RealType Real;
+ typedef typename Vector :: IndexType Index;
+
+ tnlAssert( v. getSize() > 0,
+ cerr << "Vector name is " << v. getName() );
+
+ const Index n = v. getSize();
+ for( Index i = 0; i < n; i ++ )
+ v[ i ] *= alpha;
+}
+
+
+template< typename Vector1, typename Vector2 >
+typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorSdot( const Vector1& v1,
+ const Vector2& v2 )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( v1. getSize() > 0,
+ cerr << "Vector name is " << v1. getName() );
+ tnlAssert( v1. getSize() == v2. getSize(),
+ cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
+
+ Real result = 0;
+ const Index n = v1. getSize();
+ for( Index i = 0; i < n; i ++ )
+ result += v1. getElement( i ) * v2. getElement( i );
+ return result;
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlHost > :: vectorSaxpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( x. getSize() > 0,
+ cerr << "Vector name is " << x. getName() );
+ tnlAssert( x. getSize() == y. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+ const Index n = y. getSize();
+ for( Index i = 0; i < n; i ++ )
+ y[ i ] += alpha * x[ i ];
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlHost > :: vectorSaxmy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( x. getSize() > 0,
+ cerr << "Vector name is " << x. getName() );
+ tnlAssert( x. getSize() == y. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+ const Index n = y. getSize();
+ for( Index i = 0; i < n; i ++ )
+ y[ i ] = alpha * x[ i ] - y[ i ];
+}
+
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlHost > :: vectorSaxpsby( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const typename Vector1 :: RealType& beta )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( x. getSize() > 0,
+ cerr << "Vector name is " << x. getName() );
+ tnlAssert( x. getSize() == y. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+
+ const Index n = y. getSize();
+ for( Index i = 0; i < n; i ++ )
+ y[ i ] = alpha * x[ i ] + beta * y[ i ];
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlHost > :: vectorSaxpsbz( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( x. getSize() > 0,
+ cerr << "Vector name is " << x. getName() );
+ tnlAssert( x. getSize() == y. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+ tnlAssert( x. getSize() == z. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << z. getName() );
+
+
+ const Index n = y. getSize();
+ for( Index i = 0; i < n; i ++ )
+ y[ i ] = alpha * x[ i ] + beta * z[ i ];
+}
+
+template< typename Vector1, typename Vector2 >
+void tnlVectorOperations< tnlHost > :: vectorSaxpsbzpy( Vector1& y,
+ const Vector2& x,
+ const typename Vector1 :: RealType& alpha,
+ const Vector2& z,
+ const typename Vector1 :: RealType& beta )
+{
+ typedef typename Vector1 :: RealType Real;
+ typedef typename Vector1 :: IndexType Index;
+
+ tnlAssert( x. getSize() > 0,
+ cerr << "Vector name is " << x. getName() );
+ tnlAssert( x. getSize() == y. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << y. getName() );
+ tnlAssert( x. getSize() == z. getSize(),
+ cerr << "Vector names are " << x. getName() << " and " << z. getName() );
+
+ const Index n = y. getSize();
+ for( Index i = 0; i < n; i ++ )
+ y[ i ] += alpha * x[ i ] + beta * z[ i ];
+}
+
+#endif /* TNLVECTOROPERATIONSHOST_IMPL_H_ */
diff --git a/src/implementation/core/vectors/tnlVector_impl.h b/src/implementation/core/vectors/tnlVector_impl.h
index c8770dd376b49abda9332a7eb2f4cb1ca16ff1a2..c8aeebbd7d47f361188348a7ee08e9f84157dffa 100644
--- a/src/implementation/core/vectors/tnlVector_impl.h
+++ b/src/implementation/core/vectors/tnlVector_impl.h
@@ -18,7 +18,7 @@
#ifndef TNLVECTOR_H_IMPLEMENTATION
#define TNLVECTOR_H_IMPLEMENTATION
-#include <implementation/core/vectors/vector-operations.h>
+#include <core/vectors/tnlVectorOperations.h>
template< typename Real,
typename Device,
@@ -110,7 +110,7 @@ template< typename Real,
typename Index >
Real tnlVector< Real, Device, Index > :: max() const
{
- return getVectorMax( *this );
+ return tnlVectorOperations< Device > :: getVectorMax( *this );
}
template< typename Real,
@@ -118,7 +118,7 @@ template< typename Real,
typename Index >
Real tnlVector< Real, Device, Index > :: min() const
{
- return getVectorMin( *this );
+ return tnlVectorOperations< Device > :: getVectorMin( *this );
}
@@ -127,7 +127,7 @@ template< typename Real,
typename Index >
Real tnlVector< Real, Device, Index > :: absMax() const
{
- return getVectorAbsMax( *this );
+ return tnlVectorOperations< Device > :: getVectorAbsMax( *this );
}
template< typename Real,
@@ -135,7 +135,7 @@ template< typename Real,
typename Index >
Real tnlVector< Real, Device, Index > :: absMin() const
{
- return getVectorAbsMin( *this );
+ return tnlVectorOperations< Device > :: getVectorAbsMin( *this );
}
template< typename Real,
@@ -143,7 +143,7 @@ template< typename Real,
typename Index >
Real tnlVector< Real, Device, Index > :: lpNorm( const Real& p ) const
{
- return getVectorLpNorm( *this, p );
+ return tnlVectorOperations< Device > :: getVectorLpNorm( *this, p );
}
@@ -152,7 +152,7 @@ template< typename Real,
typename Index >
Real tnlVector< Real, Device, Index > :: sum() const
{
- return getVectorSum( *this );
+ return tnlVectorOperations< Device > :: getVectorSum( *this );
}
@@ -162,7 +162,7 @@ template< typename Real,
template< typename Vector >
Real tnlVector< Real, Device, Index > :: differenceMax( const Vector& v ) const
{
- return getVectorDifferenceMax( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceMax( *this, v );
}
@@ -172,7 +172,7 @@ template< typename Real,
template< typename Vector >
Real tnlVector< Real, Device, Index > :: differenceMin( const Vector& v ) const
{
- return getVectorDifferenceMin( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceMin( *this, v );
}
@@ -182,7 +182,7 @@ template< typename Real,
template< typename Vector >
Real tnlVector< Real, Device, Index > :: differenceAbsMax( const Vector& v ) const
{
- return getVectorDifferenceAbsMax( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceAbsMax( *this, v );
}
template< typename Real,
@@ -191,7 +191,7 @@ template< typename Real,
template< typename Vector >
Real tnlVector< Real, Device, Index > :: differenceAbsMin( const Vector& v ) const
{
- return getVectorDifferenceAbsMin( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceAbsMin( *this, v );
}
template< typename Real,
@@ -200,7 +200,7 @@ template< typename Real,
template< typename Vector >
Real tnlVector< Real, Device, Index > :: differenceLpNorm( const Vector& v, const Real& p ) const
{
- return getVectorDifferenceLpNorm( *this, v, p );
+ return tnlVectorOperations< Device > :: getVectorDifferenceLpNorm( *this, v, p );
}
@@ -210,7 +210,7 @@ template< typename Real,
template< typename Vector >
Real tnlVector< Real, Device, Index > :: differenceSum( const Vector& v ) const
{
- return getVectorDifferenceSum( *this, v );
+ return tnlVectorOperations< Device > :: getVectorDifferenceSum( *this, v );
}
@@ -219,7 +219,7 @@ template< typename Real,
typename Index >
void tnlVector< Real, Device, Index > :: scalarMultiplication( const Real& alpha )
{
- vectorScalarMultiplication( *this, alpha );
+ tnlVectorOperations< Device > :: vectorScalarMultiplication( *this, alpha );
}
@@ -229,7 +229,7 @@ template< typename Real,
template< typename Vector >
Real tnlVector< Real, Device, Index > :: sdot( const Vector& v )
{
- return getVectorSdot( *this, v );
+ return tnlVectorOperations< Device > :: getVectorSdot( *this, v );
}
@@ -240,7 +240,7 @@ template< typename Vector >
void tnlVector< Real, Device, Index > :: saxpy( const Real& alpha,
const Vector& x )
{
- vectorSaxpy( *this, x, alpha );
+ tnlVectorOperations< Device > :: vectorSaxpy( *this, x, alpha );
}
template< typename Real,
@@ -250,7 +250,7 @@ template< typename Vector >
void tnlVector< Real, Device, Index > :: saxmy( const Real& alpha,
const Vector& x )
{
- vectorSaxmy( *this, x, alpha );
+ tnlVectorOperations< Device > :: vectorSaxmy( *this, x, alpha );
}
template< typename Real,
@@ -261,7 +261,7 @@ void tnlVector< Real, Device, Index > :: saxpsby( const Real& alpha,
const Vector& x,
const Real& beta )
{
- vectorSaxpsby( *this, x, alpha, beta );
+ tnlVectorOperations< Device > :: vectorSaxpsby( *this, x, alpha, beta );
}
template< typename Real,
@@ -273,7 +273,7 @@ void tnlVector< Real, Device, Index > :: saxpsbz( const Real& alpha,
const Real& beta,
const Vector& z )
{
- vectorSaxpsbz( *this, x, alpha, z, beta );
+ tnlVectorOperations< Device > :: vectorSaxpsbz( *this, x, alpha, z, beta );
}
template< typename Real,
@@ -285,7 +285,7 @@ void tnlVector< Real, Device, Index > :: saxpsbzpy( const Real& alpha,
const Real& beta,
const Vector& z )
{
- vectorSaxpsbzpy( *this, x, alpha, z, beta );
+ tnlVectorOperations< Device > :: vectorSaxpsbzpy( *this, x, alpha, z, beta );
}
diff --git a/src/implementation/core/vectors/vector-operations.h b/src/implementation/core/vectors/vector-operations.h
deleted file mode 100644
index 53e37f6547304fa42f6b22556e50a471f5582c77..0000000000000000000000000000000000000000
--- a/src/implementation/core/vectors/vector-operations.h
+++ /dev/null
@@ -1,1085 +0,0 @@
-/***************************************************************************
- vector-operations.h - description
- -------------------
- begin : Nov 8, 2012
- copyright : (C) 2012 by Tomas Oberhuber
- email : tomas.oberhuber@fjfi.cvut.cz
- ***************************************************************************/
-
-/***************************************************************************
- * *
- * This program is free software; you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation; either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- ***************************************************************************/
-
-#ifndef VECTOROPERATIONS_H_
-#define VECTOROPERATIONS_H_
-
-#include <core/cuda/device-check.h>
-#include <core/cuda/cuda-reduction.h>
-#include <core/cuda/reduction-operations.h>
-
-template< typename Vector >
-typename Vector :: RealType getHostVectorMax( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result = v. getElement( 0 );
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Max( result, v. getElement( i ) );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getCudaVectorMax( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionMax< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v. getSize(),
- v. getData(),
- ( Real* ) 0,
- result );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getVectorMax( const Vector& v )
-{
- tnlAssert( v. getSize() > 0,
- cerr << "Vector name is " << v. getName() );
- typedef typename Vector :: DeviceType Device;
- switch( Device :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorMax( v );
- case tnlCudaDevice:
- return getCudaVectorMax( v );
- }
-}
-
-template< typename Vector >
-typename Vector :: RealType getHostVectorMin( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result = v. getElement( 0 );
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Min( result, v. getElement( i ) );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getCudaVectorMin( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionMin< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v. getSize(),
- v. getData(),
- ( Real* ) 0,
- result );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getVectorMin( const Vector& v )
-{
- tnlAssert( v. getSize() > 0,
- cerr << "Vector name is " << v. getName() );
- typedef typename Vector :: DeviceType Device;
- switch( Device :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorMin( v );
- case tnlCudaDevice:
- return getCudaVectorMin( v );
- }
-}
-
-template< typename Vector >
-typename Vector :: RealType getHostVectorAbsMax( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result = fabs( v. getElement( 0 ) );
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Max( result, ( Real ) fabs( v. getElement( i ) ) );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getCudaVectorAbsMax( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionAbsMax< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v. getSize(),
- v. getData(),
- ( Real* ) 0,
- result );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getVectorAbsMax( const Vector& v )
-{
- tnlAssert( v. getSize() > 0,
- cerr << "Vector name is " << v. getName() );
- typedef typename Vector :: DeviceType Device;
- switch( Device :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorAbsMax( v );
- case tnlCudaDevice:
- return getCudaVectorAbsMax( v );
- }
-}
-
-template< typename Vector >
-typename Vector :: RealType getHostVectorAbsMin( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result = fabs( v. getElement( 0 ) );
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Min( result, ( Real ) fabs( v. getElement( i ) ) );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getCudaVectorAbsMin( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionAbsMin< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v. getSize(),
- v. getData(),
- ( Real* ) 0,
- result );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getVectorAbsMin( const Vector& v )
-{
- tnlAssert( v. getSize() > 0,
- cerr << "Vector name is " << v. getName() );
- typedef typename Vector :: DeviceType Device;
- switch( Device :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorAbsMin( v );
- case tnlCudaDevice:
- return getCudaVectorAbsMin( v );
- }
-}
-
-template< typename Vector >
-typename Vector :: RealType getHostVectorLpNorm( const Vector& v,
- const typename Vector :: RealType& p )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
- if( p == 1.0 )
- {
- Real result = fabs( v. getElement( 0 ) );
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- result += fabs( v. getElement( i ) );
- return result;
- }
- if( p == 2.0 )
- {
- Real result = v. getElement( 0 );
- result *= result;
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- {
- const Real aux = v. getElement( i );
- result += aux * aux;
- }
- return sqrt( result );
- }
- Real result = pow( fabs( v. getElement( 0 ) ), p );
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- result += pow( fabs( v. getElement( i ) ), p );
- return pow( result, 1.0 / p );
-}
-
-template< typename Vector >
-typename Vector :: RealType getCudaVectorLpNorm( const Vector& v,
- const typename Vector :: RealType& p )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
-
- Real result( 0 );
- tnlParallelReductionLpNorm< Real, Index > operation;
- operation. setPower( p );
- reductionOnCudaDevice( operation,
- v. getSize(),
- v. getData(),
- ( Real* ) 0,
- result );
- return pow( result, 1.0 / p );
-}
-
-template< typename Vector >
-typename Vector :: RealType getVectorLpNorm( const Vector& v,
- const typename Vector :: RealType& p )
-{
- tnlAssert( v. getSize() > 0,
- cerr << "Vector name is " << v. getName() );
- tnlAssert( p > 0.0,
- cerr << " p = " << p );
- typedef typename Vector :: DeviceType Device;
- switch( Device :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorLpNorm( v, p );
- case tnlCudaDevice:
- return getCudaVectorLpNorm( v, p );
- }
-}
-template< typename Vector >
-typename Vector :: RealType getHostVectorSum( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
-
- Real result = v. getElement( 0 );
- const Index n = v. getSize();
- for( Index i = 1; i < n; i ++ )
- result += v. getElement( i );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getCudaVectorSum( const Vector& v )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
-
- Real result( 0 );
- tnlParallelReductionSum< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v. getSize(),
- v. getData(),
- ( Real* ) 0,
- result );
- return result;
-}
-
-template< typename Vector >
-typename Vector :: RealType getVectorSum( const Vector& v )
-{
- tnlAssert( v. getSize() > 0,
- cerr << "Vector name is " << v. getName() );
- typedef typename Vector :: DeviceType Device;
- switch( Device :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorSum( v );
- case tnlCudaDevice:
- return getCudaVectorSum( v );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getHostVectorDifferenceMax( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- Real result = v1. getElement( 0 ) - v2. getElement( 0 );
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Max( result, v1. getElement( i ) - v2. getElement( i ) );
- return result;
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getCudaVectorDifferenceMax( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionDiffMax< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v1. getSize(),
- v1. getData(),
- v2. getData(),
- result );
- return result;
-}
-
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getVectorDifferenceMax( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( v1. getSize() > 0,
- cerr << "Vector name is " << v1. getName() );
- tnlAssert( v1. getSize() == v2. getSize(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorDifferenceMax( v1, v2 );
- case tnlCudaDevice:
- return getCudaVectorDifferenceMax( v1, v2 );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getHostVectorDifferenceMin( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- Real result = v1. getElement( 0 ) - v2. getElement( 0 );
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Min( result, v1. getElement( i ) - v2. getElement( i ) );
- return result;
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getCudaVectorDifferenceMin( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionDiffMin< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v1. getSize(),
- v1. getData(),
- v2. getData(),
- result );
- return result;
-}
-
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getVectorDifferenceMin( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( v1. getSize() > 0,
- cerr << "Vector name is " << v1. getName() );
- tnlAssert( v1. getSize() == v2. getSize(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorDifferenceMin( v1, v2 );
- case tnlCudaDevice:
- return getCudaVectorDifferenceMin( v1, v2 );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getHostVectorDifferenceAbsMax( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Max( result, ( Real ) fabs( v1. getElement( i ) - v2. getElement( i ) ) );
- return result;
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getCudaVectorDifferenceAbsMax( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionDiffAbsMax< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v1. getSize(),
- v1. getData(),
- v2. getData(),
- result );
- return result;
-}
-
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getVectorDifferenceAbsMax( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( v1. getSize() > 0,
- cerr << "Vector name is " << v1. getName() );
- tnlAssert( v1. getSize() == v2. getSize(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorDifferenceAbsMax( v1, v2 );
- case tnlCudaDevice:
- return getCudaVectorDifferenceAbsMax( v1, v2 );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getHostVectorDifferenceAbsMin( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- result = Min( result, ( Real ) fabs( v1. getElement( i ) - v2. getElement( i ) ) );
- return result;
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getCudaVectorDifferenceAbsMin( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionDiffAbsMin< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v1. getSize(),
- v1. getData(),
- v2. getData(),
- result );
- return result;
-}
-
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getVectorDifferenceAbsMin( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( v1. getSize() > 0,
- cerr << "Vector name is " << v1. getName() );
- tnlAssert( v1. getSize() == v2. getSize(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorDifferenceAbsMin( v1, v2 );
- case tnlCudaDevice:
- return getCudaVectorDifferenceAbsMin( v1, v2 );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getHostVectorDifferenceLpNorm( const Vector1& v1,
- const Vector2& v2,
- const typename Vector1 :: RealType& p )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- if( p == 1.0 )
- {
- Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- result += fabs( v1. getElement( i ) - v2. getElement( i ) );
- return result;
- }
- if( p == 2.0 )
- {
- Real result = fabs( v1. getElement( 0 ) - v2. getElement( 0 ) );
- result *= result;
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- {
- Real aux = fabs( v1. getElement( i ) - v2. getElement( i ) );
- result += aux * aux;
- }
- return sqrt( result );
- }
- Real result = pow( fabs( v1. getElement( 0 ) - v2. getElement( 0 ) ), p );
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- result += pow( fabs( v1. getElement( i ) - v2. getElement( i ) ), p );
- return pow( result, 1.0 / p );
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getCudaVectorDifferenceLpNorm( const Vector1& v1,
- const Vector2& v2,
- const typename Vector1 :: RealType& p )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionDiffLpNorm< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v1. getSize(),
- v1. getData(),
- v2. getData(),
- result );
- return result;
-}
-
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getVectorDifferenceLpNorm( const Vector1& v1,
- const Vector2& v2,
- const typename Vector1 :: RealType& p )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( p > 0.0,
- cerr << " p = " << p );
- tnlAssert( v1. getSize() > 0,
- cerr << "Vector name is " << v1. getName() );
- tnlAssert( v1. getSize() == v2. getSize(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorDifferenceLpNorm( v1, v2, p );
- case tnlCudaDevice:
- return getCudaVectorDifferenceLpNorm( v1, v2, p );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getHostVectorDifferenceSum( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- Real result = v1. getElement( 0 ) - v2. getElement( 0 );
- const Index n = v1. getSize();
- for( Index i = 1; i < n; i ++ )
- result += v1. getElement( i ) - v2. getElement( i );
- return result;
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getCudaVectorDifferenceSum( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
- Real result( 0 );
- tnlParallelReductionDiffSum< Real, Index > operation;
- reductionOnCudaDevice( operation,
- v1. getSize(),
- v1. getData(),
- v2. getData(),
- result );
- return result;
-}
-
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getVectorDifferenceSum( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( v1. getSize() > 0,
- cerr << "Vector name is " << v1. getName() );
- tnlAssert( v1. getSize() == v2. getSize(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorDifferenceSum( v1, v2 );
- case tnlCudaDevice:
- return getCudaVectorDifferenceSum( v1, v2 );
- }
-}
-
-template< typename Vector >
-void hostVectorScalarMultiplication( Vector& v,
- const typename Vector :: RealType& alpha )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
-
- const Index n = v. getSize();
- for( Index i = 0; i < n; i ++ )
- v[ i ] *= alpha;
-}
-
-template< typename Vector >
-void cudaVectorScalarMultiplication( Vector& v,
- const typename Vector :: RealType& alpha )
-{
- typedef typename Vector :: RealType Real;
- typedef typename Vector :: IndexType Index;
-
-#ifdef HAVE_CUDA
- dim3 blockSize, gridSize;
- blockSize. x = 512;
- gridSize. x = v. getSize() / 512 + 1;
- // TODO: Fix this - the grid size might be limiting for large vectors.
-
- /*
- tnlVectorCUDAScalaMultiplicationKernel<<< gridSize, blockSize >>>( v. getSize(),
- alpha,
- v. getData() );*/
-#else
- cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
-#endif
-
-}
-
-template< typename Vector >
-void vectorScalarMultiplication( Vector& v,
- const typename Vector :: RealType& alpha )
-{
- tnlAssert( v. getSize() > 0,
- cerr << "Vector name is " << v. getName() );
-
- typedef typename Vector :: DeviceType Device;
-
- switch( Device :: getDevice() )
- {
- case tnlHostDevice:
- return hostVectorScalarMultiplication( v, alpha );
- case tnlCudaDevice:
- return cudaVectorScalarMultiplication( v, alpha );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getHostVectorSdot( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- Real result = 0;
- const Index n = v1. getSize();
- for( Index i = 0; i < n; i ++ )
- result += v1. getElement( i ) * v2. getElement( i );
- return result;
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getCudaVectorSdot( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- Real result( 0 );
- /*reductionOnCudaDevice< Real,
- Real,
- Index,
- tnlParallelReductionSdot >
- ( v1. getSize(),
- v1. getData(),
- v2. getData(),
- result,
- ( Real ) 0 );*/
- return result;
-}
-
-template< typename Vector1, typename Vector2 >
-typename Vector1 :: RealType getVectorSdot( const Vector1& v1,
- const Vector2& v2 )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( v1. getSize() > 0,
- cerr << "Vector name is " << v1. getName() );
- tnlAssert( v1. getSize() == v2. getSize(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << v1. getName() << " and " << v2. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return getHostVectorSdot( v1, v2 );
- case tnlCudaDevice:
- return getCudaVectorSdot( v1, v2 );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-void hostVectorSaxpy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- const Index n = y. getSize();
- for( Index i = 0; i < n; i ++ )
- y[ i ] += alpha * x[ i ];
-}
-
-template< typename Vector1, typename Vector2 >
-void cudaVectorSaxpy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
-#ifdef HAVE_CUDA
- dim3 blockSize, gridSize;
- blockSize. x = 512;
- gridSize. x = x. getSize() / 512 + 1;
-
- // TODO: fix this
- /*tnlVectorCUDASaxpyKernel<<< gridSize, blockSize >>>( y. getSize(),
- alpha,
- x. getData(),
- y. getData() );*/
-#else
- cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
-#endif
-}
-
-template< typename Vector1, typename Vector2 >
-void vectorSaxpy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( y. getSize() > 0,
- cerr << "Vector name is " << y. getName() );
- tnlAssert( y. getSize() == x. getSize(),
- cerr << "Vector names are " << x. getName() << " and " << y. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << y. getName() << " and " << x. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return hostVectorSaxpy( y, x, alpha );
- case tnlCudaDevice:
- return cudaVectorSaxpy( y, x, alpha );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-void hostVectorSaxmy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- const Index n = y. getSize();
- for( Index i = 0; i < n; i ++ )
- y[ i ] = alpha * x[ i ] - y[ i ];
-}
-
-template< typename Vector1, typename Vector2 >
-void cudaVectorSaxmy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
-#ifdef HAVE_CUDA
- dim3 blockSize, gridSize;
- blockSize. x = 512;
- gridSize. x = x. getSize() / 512 + 1;
-
- // TODO: fix this
- /*tnlVectorCUDASaxmyKernel<<< gridSize, blockSize >>>( y. getSize(),
- alpha,
- x. getData(),
- y. getData() );*/
-#else
- cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
-#endif
-}
-
-template< typename Vector1, typename Vector2 >
-void vectorSaxmy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( y. getSize() > 0,
- cerr << "Vector name is " << y. getName() );
- tnlAssert( y. getSize() == x. getSize(),
- cerr << "Vector names are " << x. getName() << " and " << y. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << y. getName() << " and " << x. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return hostVectorSaxmy( y, x, alpha );
- case tnlCudaDevice:
- return cudaVectorSaxmy( y, x, alpha );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-void hostVectorSaxpsby( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- const Index n = y. getSize();
- for( Index i = 0; i < n; i ++ )
- y[ i ] = alpha * x[ i ] + beta * y[ i ];
-}
-
-template< typename Vector1, typename Vector2 >
-void cudaVectorSaxpsby( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
-#ifdef HAVE_CUDA
- dim3 blockSize, gridSize;
- blockSize. x = 512;
- gridSize. x = x. getSize() / 512 + 1;
-
- // TODO: fix this
- /*tnlVectorCUDASaxpsbzKernel<<< gridSize, blockSize >>>( y. getSize(),
- alpha,
- x. getData(),
- beta );*/
-#else
- cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
-#endif
-}
-
-
-template< typename Vector1, typename Vector2 >
-void vectorSaxpsby( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( y. getSize() > 0,
- cerr << "Vector name is " << y. getName() );
- tnlAssert( y. getSize() == x. getSize(),
- cerr << "Vector names are " << x. getName() << " and " << y. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << y. getName() << " and " << x. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return hostVectorSaxpsby( y, x, alpha, beta );
- case tnlCudaDevice:
- return cudaVectorSaxpsby( y, x, alpha, beta );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-void hostVectorSaxpsbz( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const Vector2& z,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- const Index n = y. getSize();
- for( Index i = 0; i < n; i ++ )
- y[ i ] = alpha * x[ i ] + beta * z[ i ];
-}
-
-template< typename Vector1, typename Vector2 >
-void cudaVectorSaxpsbz( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const Vector2& z,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
-#ifdef HAVE_CUDA
- dim3 blockSize, gridSize;
- blockSize. x = 512;
- gridSize. x = x. getSize() / 512 + 1;
-
- // TODO: fix this
- /*tnlVectorCUDASaxpsbzKernel<<< gridSize, blockSize >>>( y. getSize(),
- alpha,
- x. getData(),
- beta,
- z. getData() );*/
-#else
- cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
-#endif
-}
-
-
-template< typename Vector1, typename Vector2 >
-void vectorSaxpsbz( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const Vector2& z,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( y. getSize() > 0,
- cerr << "Vector name is " << y. getName() );
- tnlAssert( y. getSize() == x. getSize(),
- cerr << "Vector names are " << x. getName() << " and " << y. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << y. getName() << " and " << x. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return hostVectorSaxpsbz( y, x, alpha, z, beta );
- case tnlCudaDevice:
- return cudaVectorSaxpsbz( y, x, alpha, z, beta );
- }
-}
-
-template< typename Vector1, typename Vector2 >
-void hostVectorSaxpsbzpy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const Vector2& z,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
- const Index n = y. getSize();
- for( Index i = 0; i < n; i ++ )
- y[ i ] += alpha * x[ i ] + beta * z[ i ];
-}
-
-template< typename Vector1, typename Vector2 >
-void cudaVectorSaxpsbzpy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const Vector2& z,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: RealType Real;
- typedef typename Vector1 :: IndexType Index;
-
-#ifdef HAVE_CUDA
- dim3 blockSize, gridSize;
- blockSize. x = 512;
- gridSize. x = x. getSize() / 512 + 1;
-
- // TODO: fix this
- /*tnlVectorCUDASaxpsbzpyKernel<<< gridSize, blockSize >>>( y. getSize(),
- alpha,
- x. getData(),
- beta,
- z. getData() );*/
-#else
- cerr << "I am sorry but CUDA support is missing on this system " << __FILE__ << " line " << __LINE__ << "." << endl;
-#endif
-}
-
-
-template< typename Vector1, typename Vector2 >
-void vectorSaxpsbzpy( Vector1& y,
- const Vector2& x,
- const typename Vector1 :: RealType& alpha,
- const Vector2& z,
- const typename Vector1 :: RealType& beta )
-{
- typedef typename Vector1 :: DeviceType Device1;
- typedef typename Vector2 :: DeviceType Device2;
-
- tnlAssert( y. getSize() > 0,
- cerr << "Vector name is " << y. getName() );
- tnlAssert( y. getSize() == x. getSize(),
- cerr << "Vector names are " << x. getName() << " and " << y. getName() );
- tnlAssert( Device1 :: getDevice() == Device2 :: getDevice(),
- cerr << "Vector names are " << y. getName() << " and " << x. getName() );
-
- switch( Device1 :: getDevice() )
- {
- case tnlHostDevice:
- return hostVectorSaxpsbzpy( y, x, alpha, z, beta );
- case tnlCudaDevice:
- return cudaVectorSaxpsbzpy( y, x, alpha, z, beta );
- }
-}
-
-
-#endif /* VECTOROPERATIONS_H_ */
diff --git a/tests/unit-tests/CMakeLists.txt b/tests/unit-tests/CMakeLists.txt
index c392a8e5bcd5c85df77f17c08552337e064336f7..45ae2e6c57859b3f5c4d6e8b252f84760da35bb7 100755
--- a/tests/unit-tests/CMakeLists.txt
+++ b/tests/unit-tests/CMakeLists.txt
@@ -31,6 +31,7 @@ else()
tnl${mpiExt}${debugExt}-0.1 )
endif()
ADD_TEST( core/arrays/tnlArrayTest${mpiExt}${debugExt} ${EXECUTABLE_OUTPUT_PATH}/tnlArrayTest${mpiExt}${debugExt} )
+ADD_TEST( core/vectors/tnlVectorOperationsTest${mpiExt}${debugExt} ${EXECUTABLE_OUTPUT_PATH}/tnlVectorOperationsTest${mpiExt}${debugExt} )
if( BUILD_CUDA )
ADD_TEST( core/cuda/tnl-device-check-test${mpiExt}${debugExt} ${EXECUTABLE_OUTPUT_PATH}/tnl-device-check-test${mpiExt}${debugExt} )
diff --git a/tests/unit-tests/core/vectors/CMakeLists.txt b/tests/unit-tests/core/vectors/CMakeLists.txt
index c1d5649581c59c9d479e6cb700dd97bad1ef7946..e793b1c3a046f8578b86af4f42ebc210f66b9e4b 100755
--- a/tests/unit-tests/core/vectors/CMakeLists.txt
+++ b/tests/unit-tests/core/vectors/CMakeLists.txt
@@ -1,6 +1,6 @@
set( headers tnlVectorTester.h
tnlSharedVectorTester.h
- tnlCudaVectorOperationsTester.h )
+ tnlVectorOperationsTester.h )
if( BUILD_CUDA )
CUDA_ADD_EXECUTABLE( tnlCudaVectorOperationsTest${mpiExt}${debugExt} ${headers} tnlCudaVectorOperationsTest.cu )
@@ -8,6 +8,13 @@ if( BUILD_CUDA )
tnl${mpiExt}${debugExt}-0.1 )
endif()
+
+ADD_EXECUTABLE( tnlVectorOperationsTest${mpiExt}${debugExt} ${headers} tnlVectorOperationsTest.cpp )
+TARGET_LINK_LIBRARIES( tnlVectorOperationsTest${mpiExt}${debugExt} ${CPPUNIT_LIBRARIES}
+ tnl${mpiExt}${debugExt}-0.1 )
+
+
+
#ADD_EXECUTABLE( tnlVectorTest${mpiExt}${debugExt} ${headers} tnlVectorTest.cpp )
#TARGET_LINK_LIBRARIES( tnlVectorTest${mpiExt}${debugExt} ${CPPUNIT_LIBRARIES}
# tnl${mpiExt}${debugExt}-0.1 )
diff --git a/tests/unit-tests/core/vectors/tnlCudaVectorOperationsTest.cu b/tests/unit-tests/core/vectors/tnlCudaVectorOperationsTest.cu
index 4e9af3b2ef973855255de9ae0789231149d2f656..0223a5a7c3a8d49c7d6b3d85991bb44610043c82 100644
--- a/tests/unit-tests/core/vectors/tnlCudaVectorOperationsTest.cu
+++ b/tests/unit-tests/core/vectors/tnlCudaVectorOperationsTest.cu
@@ -15,12 +15,12 @@
* *
***************************************************************************/
-#include "tnlCudaVectorOperationsTester.h"
+#include "tnlVectorOperationsTester.h"
#include "../../tnlUnitTestStarter.h"
int main( int argc, char* argv[] )
{
- if( ! tnlUnitTestStarter :: run< tnlCudaVectorOperationsTester< double > >() )
+ if( ! tnlUnitTestStarter :: run< tnlVectorOperationsTester< double, tnlCuda > >() )
return EXIT_FAILURE;
return EXIT_SUCCESS;
}
\ No newline at end of file
diff --git a/tests/unit-tests/core/vectors/tnlVectorOperationsTest.cpp b/tests/unit-tests/core/vectors/tnlVectorOperationsTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..08580d1c455ca63d08274c09cde90c47ef406383
--- /dev/null
+++ b/tests/unit-tests/core/vectors/tnlVectorOperationsTest.cpp
@@ -0,0 +1,30 @@
+/***************************************************************************
+ tnlVectorOperationsTest.cpp - description
+ -------------------
+ begin : Jul 15, 2013
+ copyright : (C) 2013 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#include "tnlVectorOperationsTester.h"
+#include "../../tnlUnitTestStarter.h"
+#include <core/tnlHost.h>
+
+int main( int argc, char* argv[] )
+{
+ if( ! tnlUnitTestStarter :: run< tnlVectorOperationsTester< double, tnlHost > >() )
+ return EXIT_FAILURE;
+ return EXIT_SUCCESS;
+}
+
+
+
diff --git a/tests/unit-tests/core/vectors/tnlCudaVectorOperationsTester.h b/tests/unit-tests/core/vectors/tnlVectorOperationsTester.h
similarity index 56%
rename from tests/unit-tests/core/vectors/tnlCudaVectorOperationsTester.h
rename to tests/unit-tests/core/vectors/tnlVectorOperationsTester.h
index 7d306eb5cd875b7b6aeda9288d2ee1c0d3985781..9da836ff19a8dc2fd94e64136afe56198b69a425 100644
--- a/tests/unit-tests/core/vectors/tnlCudaVectorOperationsTester.h
+++ b/tests/unit-tests/core/vectors/tnlVectorOperationsTester.h
@@ -1,5 +1,5 @@
/***************************************************************************
- tnlCudaVectorOperationsTester.h - description
+ tnlVectorOperationsTester.h - description
-------------------
begin : Mar 30, 2013
copyright : (C) 2013 by Tomas Oberhuber
@@ -15,8 +15,8 @@
* *
***************************************************************************/
-#ifndef TNLCUDAVECTOROPERATIONSTESTER_H_
-#define TNLCUDAVECTOROPERATIONSTESTER_H_
+#ifndef TNLVECTOROPERATIONSTESTER_H_
+#define TNLVECTOROPERATIONSTESTER_H_
#include <tnlConfig.h>
@@ -26,68 +26,67 @@
#include <cppunit/TestCaller.h>
#include <cppunit/TestCase.h>
#include <cppunit/Message.h>
-#include <core/cuda/device-check.h>
-#include <implementation/core/memory-operations.h>
+
#include <core/vectors/tnlVector.h>
-#include <implementation/core/vectors/vector-operations.h>
+#include <core/vectors/tnlVectorOperations.h>
-template< typename Type >
-class tnlCudaVectorOperationsTester : public CppUnit :: TestCase
+template< typename Real, typename Device >
+class tnlVectorOperationsTester : public CppUnit :: TestCase
{
public:
- tnlCudaVectorOperationsTester(){};
+ tnlVectorOperationsTester(){};
virtual
- ~tnlCudaVectorOperationsTester(){};
+ ~tnlVectorOperationsTester(){};
static CppUnit :: Test* suite()
{
- CppUnit :: TestSuite* suiteOfTests = new CppUnit :: TestSuite( "tnlCudaVectorOperationsTester" );
+ CppUnit :: TestSuite* suiteOfTests = new CppUnit :: TestSuite( "tnlVectorOperationsTester" );
CppUnit :: TestResult result;
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorMaxTest",
- &tnlCudaVectorOperationsTester :: getVectorMaxTest )
+ &tnlVectorOperationsTester :: getVectorMaxTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorMinTest",
- &tnlCudaVectorOperationsTester :: getVectorMinTest )
+ &tnlVectorOperationsTester :: getVectorMinTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorAbsMaxTest",
- &tnlCudaVectorOperationsTester :: getVectorAbsMaxTest )
+ &tnlVectorOperationsTester :: getVectorAbsMaxTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorAbsMinTest",
- &tnlCudaVectorOperationsTester :: getVectorAbsMinTest )
+ &tnlVectorOperationsTester :: getVectorAbsMinTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorLpNormTest",
- &tnlCudaVectorOperationsTester :: getVectorLpNormTest )
+ &tnlVectorOperationsTester :: getVectorLpNormTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorSumTest",
- &tnlCudaVectorOperationsTester :: getVectorSumTest )
+ &tnlVectorOperationsTester :: getVectorSumTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorDifferenceMaxTest",
- &tnlCudaVectorOperationsTester :: getVectorDifferenceMaxTest )
+ &tnlVectorOperationsTester :: getVectorDifferenceMaxTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorDifferenceMinTest",
- &tnlCudaVectorOperationsTester :: getVectorDifferenceMinTest )
+ &tnlVectorOperationsTester :: getVectorDifferenceMinTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorDifferenceAbsMaxTest",
- &tnlCudaVectorOperationsTester :: getVectorDifferenceAbsMaxTest )
+ &tnlVectorOperationsTester :: getVectorDifferenceAbsMaxTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorDifferenceAbsMinTest",
- &tnlCudaVectorOperationsTester :: getVectorDifferenceAbsMinTest )
+ &tnlVectorOperationsTester :: getVectorDifferenceAbsMinTest )
);
- suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlCudaVectorOperationsTester >(
+ suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlVectorOperationsTester >(
"getVectorDifferenceLpNormTest",
- &tnlCudaVectorOperationsTester :: getVectorDifferenceLpNormTest )
+ &tnlVectorOperationsTester :: getVectorDifferenceLpNormTest )
);
return suiteOfTests;
@@ -141,139 +140,139 @@ class tnlCudaVectorOperationsTester : public CppUnit :: TestCase
void getVectorMaxTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > v;
+ tnlVector< Real, Device > v;
v. setSize( size );
setLinearSequence( v );
- CPPUNIT_ASSERT( getCudaVectorMax( v ) == size - 1 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorMax( v ) == size - 1 );
}
void getVectorMinTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > v;
+ tnlVector< Real, Device > v;
v. setSize( size );
setLinearSequence( v );
- CPPUNIT_ASSERT( getCudaVectorMin( v ) == 0 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorMin( v ) == 0 );
}
void getVectorAbsMaxTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > v;
+ tnlVector< Real, Device > v;
v. setSize( size );
setNegativeLinearSequence( v );
- CPPUNIT_ASSERT( getCudaVectorAbsMax( v ) == size - 1 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorAbsMax( v ) == size - 1 );
}
void getVectorAbsMinTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > v;
+ tnlVector< Real, Device > v;
v. setSize( size );
setNegativeLinearSequence( v );
- CPPUNIT_ASSERT( getCudaVectorAbsMin( v ) == 0 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorAbsMin( v ) == 0 );
}
void getVectorLpNormTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > v;
+ tnlVector< Real, Device > v;
v. setSize( size );
setOnesSequence( v );
- CPPUNIT_ASSERT( getCudaVectorLpNorm( v, 2.0 ) == sqrt( size ) );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorLpNorm( v, 2.0 ) == sqrt( size ) );
}
void getVectorSumTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > v;
+ tnlVector< Real, Device > v;
v. setSize( size );
setOnesSequence( v );
- CPPUNIT_ASSERT( getCudaVectorSum( v ) == size );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorSum( v ) == size );
setLinearSequence( v );
- CPPUNIT_ASSERT( getCudaVectorSum( v ) == ( ( Type ) size ) * ( ( Type ) size - 1 ) / 2 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorSum( v ) == ( ( Real ) size ) * ( ( Real ) size - 1 ) / 2 );
}
void getVectorDifferenceMaxTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > u, v;
+ tnlVector< Real, Device > u, v;
u. setSize( size );
v. setSize( size );
setLinearSequence( u );
setOnesSequence( v );
- CPPUNIT_ASSERT( getCudaVectorDifferenceMax( u, v ) == size - 2 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceMax( u, v ) == size - 2 );
}
void getVectorDifferenceMinTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > u, v;
+ tnlVector< Real, Device > u, v;
u. setSize( size );
v. setSize( size );
setLinearSequence( u );
setOnesSequence( v );
- CPPUNIT_ASSERT( getCudaVectorDifferenceMin( u, v ) == -1 );
- CPPUNIT_ASSERT( getCudaVectorDifferenceMin( v, u ) == -1234565 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceMin( u, v ) == -1 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceMin( v, u ) == -1234565 );
}
void getVectorDifferenceAbsMaxTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > u, v;
+ tnlVector< Real, Device > u, v;
u. setSize( size );
v. setSize( size );
setNegativeLinearSequence( u );
setOnesSequence( v );
- CPPUNIT_ASSERT( getCudaVectorDifferenceAbsMax( u, v ) == size );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceAbsMax( u, v ) == size );
}
void getVectorDifferenceAbsMinTest()
{
const int size( 1234567 );
- tnlVector< Type, tnlCuda > u, v;
+ tnlVector< Real, Device > u, v;
u. setSize( size );
v. setSize( size );
setLinearSequence( u );
setOnesSequence( v );
- CPPUNIT_ASSERT( getCudaVectorDifferenceAbsMin( u, v ) == 0 );
- CPPUNIT_ASSERT( getCudaVectorDifferenceAbsMin( v, u ) == 0 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceAbsMin( u, v ) == 0 );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceAbsMin( v, u ) == 0 );
}
void getVectorDifferenceLpNormTest()
{
const int size( 1024 );
- tnlVector< Type, tnlCuda > u, v;
+ tnlVector< Real, Device > u, v;
u. setSize( size );
v. setSize( size );
u. setValue( 3.0 );
v. setValue( 1.0 );
- cout << getCudaVectorDifferenceLpNorm( u, v, 1.0 ) << " " << 2.0 * size << endl;
- CPPUNIT_ASSERT( getCudaVectorDifferenceLpNorm( u, v, 1.0 ) == 2.0 * size );
- CPPUNIT_ASSERT( getCudaVectorDifferenceLpNorm( u, v, 2.0 ) == sqrt( 4.0 * size ) );
+ cout << tnlVectorOperations< Device > :: getVectorDifferenceLpNorm( u, v, 1.0 ) << " " << 2.0 * size << endl;
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceLpNorm( u, v, 1.0 ) == 2.0 * size );
+ CPPUNIT_ASSERT( tnlVectorOperations< Device > :: getVectorDifferenceLpNorm( u, v, 2.0 ) == sqrt( 4.0 * size ) );
}
};
#else
-template< typename Type >
-class tnlCudaVectorOperationsTester
+template< typename Device >
+class tnlVectorOperationsTester
{};
#endif /* HAVE_CPPUNIT */
-#endif /* TNLCUDAVECTOROPERATIONSTESTER_H_ */
+#endif /* TNLVECTOROPERATIONSTESTER_H_ */