diff --git a/src/core/cuda/CMakeLists.txt b/src/core/cuda/CMakeLists.txt index 3530dda77a5aed1f1b2acdbd160af0a26555b74d..b1e0aababd34c9be5196d865e5ce6b29f38adfcd 100755 --- a/src/core/cuda/CMakeLists.txt +++ b/src/core/cuda/CMakeLists.txt @@ -20,6 +20,7 @@ IF( BUILD_CUDA ) ${CURRENT_DIR}/cuda-reduction-abs-max_impl.cu ${CURRENT_DIR}/cuda-reduction-and_impl.cu ${CURRENT_DIR}/cuda-reduction-or_impl.cu + ${CURRENT_DIR}/cuda-reduction-l2-norm_impl.cu ${CURRENT_DIR}/cuda-reduction-lp-norm_impl.cu ${CURRENT_DIR}/cuda-reduction-equalities_impl.cu ${CURRENT_DIR}/cuda-reduction-inequalities_impl.cu @@ -30,6 +31,7 @@ IF( BUILD_CUDA ) ${CURRENT_DIR}/cuda-reduction-diff-abs-sum_impl.cu ${CURRENT_DIR}/cuda-reduction-diff-abs-min_impl.cu ${CURRENT_DIR}/cuda-reduction-diff-abs-max_impl.cu + ${CURRENT_DIR}/cuda-reduction-diff-l2-norm_impl.cu ${CURRENT_DIR}/cuda-reduction-diff-lp-norm_impl.cu ${CURRENT_DIR}/cuda-prefix-sum_impl.cu PARENT_SCOPE ) diff --git a/src/core/cuda/cuda-reduction-diff-l2-norm_impl.cu b/src/core/cuda/cuda-reduction-diff-l2-norm_impl.cu new file mode 100644 index 0000000000000000000000000000000000000000..4f3e95b74945960d11a389c4732fd22c845ca2df --- /dev/null +++ b/src/core/cuda/cuda-reduction-diff-l2-norm_impl.cu @@ -0,0 +1,87 @@ +/*************************************************************************** + cuda-reduction-diff-lp-norm_impl.cu - description + ------------------- + begin : Jan 19, 2014 + copyright : (C) 2014 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 <core/cuda/reduction-operations.h> +#include <core/cuda/cuda-reduction.h> + +#ifdef TEMPLATE_EXPLICIT_INSTANTIATION + +/**** + * Diff L2 Norm + */ +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< float, int > > + ( const tnlParallelReductionDiffL2Norm< float, int >& operation, + const typename tnlParallelReductionDiffL2Norm< float, int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< float, int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< float, int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< float, int> :: ResultType& result ); + +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< double, int > > + ( const tnlParallelReductionDiffL2Norm< double, int>& operation, + const typename tnlParallelReductionDiffL2Norm< double, int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< double, int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< double, int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< double, int> :: ResultType& result ); + +#ifdef INSTANTIATE_LONG_DOUBLE +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< long double, int > > + ( const tnlParallelReductionDiffL2Norm< long double, int>& operation, + const typename tnlParallelReductionDiffL2Norm< long double, int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< long double, int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< long double, int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< long double, int> :: ResultType& result ); +#endif + +#ifdef INSTANTIATE_LONG_INT +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< char, long int > > + ( const tnlParallelReductionDiffL2Norm< char, long int >& operation, + const typename tnlParallelReductionDiffL2Norm< char, long int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< char, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< char, long int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< char, long int > :: ResultType& result ); + +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< int, long int > > + ( const tnlParallelReductionDiffL2Norm< int, long int >& operation, + const typename tnlParallelReductionDiffL2Norm< int, long int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< int, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< int, long int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< int, long int > :: ResultType& result ); + +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< float, long int > > + ( const tnlParallelReductionDiffL2Norm< float, long int >& operation, + const typename tnlParallelReductionDiffL2Norm< float, long int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< float, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< float, long int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< float, long int> :: ResultType& result ); + +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< double, long int > > + ( const tnlParallelReductionDiffL2Norm< double, long int>& operation, + const typename tnlParallelReductionDiffL2Norm< double, long int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< double, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< double, long int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< double, long int> :: ResultType& result ); + +#ifdef INSTANTIATE_LONG_DOUBLE +template bool reductionOnCudaDevice< tnlParallelReductionDiffL2Norm< long double, long int > > + ( const tnlParallelReductionDiffL2Norm< long double, long int>& operation, + const typename tnlParallelReductionDiffL2Norm< long double, long int > :: IndexType size, + const typename tnlParallelReductionDiffL2Norm< long double, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionDiffL2Norm< long double, long int > :: RealType* deviceInput2, + typename tnlParallelReductionDiffL2Norm< long double, long int> :: ResultType& result ); +#endif +#endif +#endif diff --git a/src/core/cuda/cuda-reduction-l2-norm_impl.cu b/src/core/cuda/cuda-reduction-l2-norm_impl.cu new file mode 100644 index 0000000000000000000000000000000000000000..29e6b265c134e4113480f1836e65960f45b0f61f --- /dev/null +++ b/src/core/cuda/cuda-reduction-l2-norm_impl.cu @@ -0,0 +1,80 @@ +/*************************************************************************** + cuda-reduction-lp-norm_impl.cu - description + ------------------- + begin : Jan 19, 2014 + copyright : (C) 2014 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 <core/cuda/reduction-operations.h> +#include <core/cuda/cuda-reduction.h> + +#ifdef TEMPLATE_EXPLICIT_INSTANTIATION + +/**** + * L2 Norm + */ +template bool reductionOnCudaDevice< tnlParallelReductionL2Norm< float, int > > + ( const tnlParallelReductionL2Norm< float, int >& operation, + const typename tnlParallelReductionL2Norm< float, int > :: IndexType size, + const typename tnlParallelReductionL2Norm< float, int > :: RealType* deviceInput1, + const typename tnlParallelReductionL2Norm< float, int > :: RealType* deviceInput2, + typename tnlParallelReductionL2Norm< float, int> :: ResultType& result ); + +template bool reductionOnCudaDevice< tnlParallelReductionL2Norm< double, int > > + ( const tnlParallelReductionL2Norm< double, int>& operation, + const typename tnlParallelReductionL2Norm< double, int > :: IndexType size, + const typename tnlParallelReductionL2Norm< double, int > :: RealType* deviceInput1, + const typename tnlParallelReductionL2Norm< double, int > :: RealType* deviceInput2, + typename tnlParallelReductionL2Norm< double, int> :: ResultType& result ); + +#ifdef INSTANTIATE_LONG_DOUBLE +template bool reductionOnCudaDevice< tnlParallelReductionL2Norm< long double, int > > + ( const tnlParallelReductionL2Norm< long double, int>& operation, + const typename tnlParallelReductionL2Norm< long double, int > :: IndexType size, + const typename tnlParallelReductionL2Norm< long double, int > :: RealType* deviceInput1, + const typename tnlParallelReductionL2Norm< long double, int > :: RealType* deviceInput2, + typename tnlParallelReductionL2Norm< long double, int> :: ResultType& result ); +#endif + +#ifdef INSTANTIATE_LONG_INT +template bool reductionOnCudaDevice< tnlParallelReductionL2Norm< int, long int > > + ( const tnlParallelReductionL2Norm< int, long int >& operation, + const typename tnlParallelReductionL2Norm< int, long int > :: IndexType size, + const typename tnlParallelReductionL2Norm< int, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionL2Norm< int, long int > :: RealType* deviceInput2, + typename tnlParallelReductionL2Norm< int, long int> :: ResultType& result ); + +template bool reductionOnCudaDevice< tnlParallelReductionL2Norm< float, long int > > + ( const tnlParallelReductionL2Norm< float, long int >& operation, + const typename tnlParallelReductionL2Norm< float, long int > :: IndexType size, + const typename tnlParallelReductionL2Norm< float, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionL2Norm< float, long int > :: RealType* deviceInput2, + typename tnlParallelReductionL2Norm< float, long int> :: ResultType& result ); + +template bool reductionOnCudaDevice< tnlParallelReductionL2Norm< double, long int > > + ( const tnlParallelReductionL2Norm< double, long int>& operation, + const typename tnlParallelReductionL2Norm< double, long int > :: IndexType size, + const typename tnlParallelReductionL2Norm< double, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionL2Norm< double, long int > :: RealType* deviceInput2, + typename tnlParallelReductionL2Norm< double, long int> :: ResultType& result ); + +#ifdef INSTANTIATE_LONG_DOUBLE +template bool reductionOnCudaDevice< tnlParallelReductionL2Norm< long double, long int > > + ( const tnlParallelReductionL2Norm< long double, long int>& operation, + const typename tnlParallelReductionL2Norm< long double, long int > :: IndexType size, + const typename tnlParallelReductionL2Norm< long double, long int > :: RealType* deviceInput1, + const typename tnlParallelReductionL2Norm< long double, long int > :: RealType* deviceInput2, + typename tnlParallelReductionL2Norm< long double, long int> :: ResultType& result ); +#endif +#endif +#endif \ No newline at end of file diff --git a/src/core/cuda/reduction-operations.h b/src/core/cuda/reduction-operations.h index 93de5076d8fad0fd24c41c60d5f829cadd5b27c9..8733bcd88120154888be10c1f4742e38e5e06c55 100644 --- a/src/core/cuda/reduction-operations.h +++ b/src/core/cuda/reduction-operations.h @@ -501,6 +501,35 @@ class tnlParallelReductionAbsMax : public tnlParallelReductionMax< Real, Index > } }; +template< typename Real, typename Index > +class tnlParallelReductionL2Norm : public tnlParallelReductionSum< Real, Index > +{ + public: + + typedef Real RealType; + typedef Index IndexType; + typedef Real ResultType; + typedef tnlParallelReductionSum< Real, Index > LaterReductionOperation; + + ResultType reduceOnHost( const IndexType idx, + const ResultType& current, + const RealType* data1, + const RealType* data2 ) const + { + return current + data1[ idx ] * data1[ idx ]; + }; + + __cuda_callable__ ResultType initialValue() const { return ( ResultType ) 0; }; + + __cuda_callable__ void cudaFirstReduction( ResultType& result, + const IndexType index, + const RealType* data1, + const RealType* data2 ) const + { + result += data1[ index ] * data1[ index ]; + } +}; + template< typename Real, typename Index > class tnlParallelReductionLpNorm : public tnlParallelReductionSum< Real, Index > @@ -801,6 +830,37 @@ class tnlParallelReductionDiffAbsMax : public tnlParallelReductionMax< Real, Ind } }; +template< typename Real, typename Index > +class tnlParallelReductionDiffL2Norm : public tnlParallelReductionSum< Real, Index > +{ + public: + + typedef Real RealType; + typedef Index IndexType; + typedef Real ResultType; + typedef tnlParallelReductionSum< Real, Index > LaterReductionOperation; + + ResultType reduceOnHost( const IndexType idx, + const ResultType& current, + const RealType* data1, + const RealType* data2 ) const + { + const RealType aux( data2[ idx ] - data1[ idx ] ); + return current + aux * aux; + }; + + __cuda_callable__ ResultType initialValue() const { return ( ResultType ) 0; }; + + __cuda_callable__ void cudaFirstReduction( ResultType& result, + const IndexType index, + const RealType* data1, + const RealType* data2 ) const + { + const RealType aux( data2[ index ] - data1[ index ] ); + result += aux * aux; + } +}; + template< typename Real, typename Index > class tnlParallelReductionDiffLpNorm : public tnlParallelReductionSum< Real, Index > { diff --git a/src/core/vectors/tnlVectorOperations.h b/src/core/vectors/tnlVectorOperations.h index 3438b5250ff5a39dfc3529c5cf452938ff59cdbc..afa7b2e96249ae383cd9f982b98e303c2d559c20 100644 --- a/src/core/vectors/tnlVectorOperations.h +++ b/src/core/vectors/tnlVectorOperations.h @@ -54,6 +54,12 @@ class tnlVectorOperations< tnlHost > template< typename Vector > static typename Vector::RealType getVectorAbsMin( const Vector& v ); + template< typename Vector > + static typename Vector::RealType getVectorL1Norm( const Vector& v ); + + template< typename Vector > + static typename Vector::RealType getVectorL2Norm( const Vector& v ); + template< typename Vector > static typename Vector::RealType getVectorLpNorm( const Vector& v, const typename Vector::RealType& p ); @@ -77,6 +83,14 @@ class tnlVectorOperations< tnlHost > static typename Vector1::RealType getVectorDifferenceAbsMin( const Vector1& v1, const Vector2& v2 ); + template< typename Vector1, typename Vector2 > + static typename Vector1::RealType getVectorDifferenceL1Norm( const Vector1& v1, + const Vector2& v2 ); + + template< typename Vector1, typename Vector2 > + static typename Vector1::RealType getVectorDifferenceL2Norm( const Vector1& v1, + const Vector2& v2 ); + template< typename Vector1, typename Vector2 > static typename Vector1::RealType getVectorDifferenceLpNorm( const Vector1& v1, const Vector2& v2, @@ -85,6 +99,8 @@ class tnlVectorOperations< tnlHost > 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 ); @@ -146,11 +162,17 @@ class tnlVectorOperations< tnlCuda > template< typename Vector > static typename Vector::RealType getVectorAbsMin( const Vector& v ); - + + template< typename Vector > + static typename Vector::RealType getVectorL1Norm( const Vector& v ); + + template< typename Vector > + static typename Vector::RealType getVectorL2Norm( 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 ); @@ -168,8 +190,16 @@ class tnlVectorOperations< tnlCuda > template< typename Vector1, typename Vector2 > static typename Vector1::RealType getVectorDifferenceAbsMin( const Vector1& v1, - const Vector2& v2 ); + const Vector2& v2 ); + + template< typename Vector1, typename Vector2 > + static typename Vector1::RealType getVectorDifferenceL1Norm( const Vector1& v1, + const Vector2& v2 ); + template< typename Vector1, typename Vector2 > + static typename Vector1::RealType getVectorDifferenceL2Norm( const Vector1& v1, + const Vector2& v2 ); + template< typename Vector1, typename Vector2 > static typename Vector1::RealType getVectorDifferenceLpNorm( const Vector1& v1, const Vector2& v2, @@ -178,6 +208,7 @@ class tnlVectorOperations< tnlCuda > 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 ); diff --git a/src/core/vectors/tnlVectorOperationsCuda_impl.h b/src/core/vectors/tnlVectorOperationsCuda_impl.h index ccd191e63fd521f48f63d44039525e9d844071cd..d026378a8567300063a01a221fa7091fcba278f7 100644 --- a/src/core/vectors/tnlVectorOperationsCuda_impl.h +++ b/src/core/vectors/tnlVectorOperationsCuda_impl.h @@ -112,8 +112,51 @@ typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorAbsMin( c } template< typename Vector > -typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorLpNorm( const Vector& v, - const typename Vector :: RealType& p ) +typename Vector::RealType +tnlVectorOperations< tnlCuda >:: +getVectorL1Norm( const Vector& v ) +{ + typedef typename Vector :: RealType Real; + typedef typename Vector :: IndexType Index; + + tnlAssert( v. getSize() > 0, ); + + Real result( 0 ); + tnlParallelReductionAbsSum< Real, Index > operation; + reductionOnCudaDevice( operation, + v. getSize(), + v. getData(), + ( Real* ) 0, + result ); + return result; +} + +template< typename Vector > +typename Vector::RealType +tnlVectorOperations< tnlCuda >:: +getVectorL2Norm( const Vector& v ) +{ + typedef typename Vector :: RealType Real; + typedef typename Vector :: IndexType Index; + + tnlAssert( v. getSize() > 0, ); + + Real result( 0 ); + tnlParallelReductionL2Norm< Real, Index > operation; + reductionOnCudaDevice( operation, + v. getSize(), + v. getData(), + ( Real* ) 0, + result ); + return sqrt( 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; @@ -121,7 +164,11 @@ typename Vector :: RealType tnlVectorOperations< tnlCuda > :: getVectorLpNorm( c tnlAssert( v. getSize() > 0, ); tnlAssert( p > 0.0, cerr << " p = " << p ); - + + if( p == 1 ) + return getVectorL1Norm( v ); + if( p == 2 ) + return getVectorL2Norm( v ); Real result( 0 ); tnlParallelReductionLpNorm< Real, Index > operation; operation. setPower( p ); @@ -232,6 +279,51 @@ typename Vector1 :: RealType tnlVectorOperations< tnlCuda > :: getVectorDifferen return result; } +template< typename Vector1, typename Vector2 > +typename Vector1::RealType +tnlVectorOperations< tnlCuda >:: +getVectorDifferenceL1Norm( const Vector1& v1, + const Vector2& v2 ) +{ + typedef typename Vector1 :: RealType Real; + typedef typename Vector1 :: IndexType Index; + + tnlAssert( v1. getSize() > 0, ); + tnlAssert( v1. getSize() == v2. getSize(), ); + + Real result( 0 ); + tnlParallelReductionDiffAbsSum< Real, Index > operation; + reductionOnCudaDevice( operation, + v1. getSize(), + v1. getData(), + v2. getData(), + result ); + return result; +} + +template< typename Vector1, typename Vector2 > +typename Vector1::RealType +tnlVectorOperations< tnlCuda >:: +getVectorDifferenceL2Norm( const Vector1& v1, + const Vector2& v2 ) +{ + typedef typename Vector1 :: RealType Real; + typedef typename Vector1 :: IndexType Index; + + tnlAssert( v1. getSize() > 0, ); + tnlAssert( v1. getSize() == v2. getSize(), ); + + Real result( 0 ); + tnlParallelReductionDiffL2Norm< Real, Index > operation; + reductionOnCudaDevice( operation, + v1. getSize(), + v1. getData(), + v2. getData(), + result ); + return sqrt( result ); +} + + template< typename Vector1, typename Vector2 > typename Vector1::RealType tnlVectorOperations< tnlCuda >:: diff --git a/src/core/vectors/tnlVectorOperationsHost_impl.h b/src/core/vectors/tnlVectorOperationsHost_impl.h index 25de5f300e536a261b0528ac1ea40e719eda07e5..67ed621ecfa06381c55a016109e0e83ddfcdd951 100644 --- a/src/core/vectors/tnlVectorOperationsHost_impl.h +++ b/src/core/vectors/tnlVectorOperationsHost_impl.h @@ -90,10 +90,46 @@ typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorAbsMin( c return result; } +template< typename Vector > +typename Vector::RealType +tnlVectorOperations< tnlHost >:: +getVectorL1Norm( const Vector& v ) +{ + typedef typename Vector :: RealType Real; + typedef typename Vector :: IndexType Index; + tnlAssert( v. getSize() > 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; +} template< typename Vector > -typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorLpNorm( const Vector& v, - const typename Vector :: RealType& p ) +typename Vector::RealType +tnlVectorOperations< tnlHost >:: +getVectorL2Norm( const Vector& v ) +{ + typedef typename Vector :: RealType Real; + typedef typename Vector :: IndexType Index; + tnlAssert( v. getSize() > 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 ); +} + +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; @@ -101,25 +137,10 @@ typename Vector :: RealType tnlVectorOperations< tnlHost > :: getVectorLpNorm( c 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; - } + return getVectorL1Norm( v ); 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 ); - } + return getVectorL2Norm( v ); + Real result = pow( fabs( v. getElement( 0 ) ), p ); const Index n = v. getSize(); for( Index i = 1; i < n; i ++ ) @@ -207,11 +228,55 @@ typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferen return result; } +template< typename Vector1, typename Vector2 > +typename Vector1::RealType +tnlVectorOperations< tnlHost >:: +getVectorDifferenceL1Norm( const Vector1& v1, + const Vector2& v2 ) +{ + typedef typename Vector1 :: RealType Real; + typedef typename Vector1 :: IndexType Index; + + tnlAssert( v1. getSize() > 0, ); + tnlAssert( v1. getSize() == v2. getSize(), ); + + 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; +} + +template< typename Vector1, typename Vector2 > +typename Vector1::RealType +tnlVectorOperations< tnlHost >:: +getVectorDifferenceL2Norm( const Vector1& v1, + const Vector2& v2 ) +{ + typedef typename Vector1 :: RealType Real; + typedef typename Vector1 :: IndexType Index; + + tnlAssert( v1. getSize() > 0, ); + tnlAssert( v1. getSize() == v2. getSize(), ); + + 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 ); +} + template< typename Vector1, typename Vector2 > -typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceLpNorm( const Vector1& v1, - const Vector2& v2, - const typename Vector1 :: RealType& p ) +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; @@ -222,25 +287,10 @@ typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferen tnlAssert( v1. getSize() == v2. getSize(), ); 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; - } + return getVectorDifferenceL1Norm( v1, v2 ); 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 ); - } + return getVectorDifferenceL2Norm( v1, v2 ); + Real result = pow( fabs( v1. getElement( 0 ) - v2. getElement( 0 ) ), p ); const Index n = v1. getSize(); for( Index i = 1; i < n; i ++ ) @@ -249,8 +299,8 @@ typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferen } template< typename Vector1, typename Vector2 > -typename Vector1 :: RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceSum( const Vector1& v1, - const Vector2& v2 ) +typename Vector1::RealType tnlVectorOperations< tnlHost > :: getVectorDifferenceSum( const Vector1& v1, + const Vector2& v2 ) { typedef typename Vector1 :: RealType Real; typedef typename Vector1 :: IndexType Index;