Commit 9268e568 authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Fixed VectorOperationsMIC_impl.h 

- Use math function from TNL/Math.h instead of the standard library
- Many cosmetic changes for readability
parent 2fdbf71c

src/TNL/Containers/Algorithms/VectorOperationsMIC_impl.h

+318 −298
Original line number Diff line number Diff line
@@ -13,7 +13,7 @@ 
#pragma once



#include <TNL/Devices/MIC.h>



#include <TNL/Math.h>



namespace TNL {

namespace Containers {
@@ -22,7 +22,9 @@ namespace Algorithms { 
//static const int OpenMPVectorOperationsThreshold = 65536; // TODO: check this threshold



template< typename Vector >

void VectorOperations< Devices::MIC >::addElement( Vector& v,

void

VectorOperations< Devices::MIC >::

addElement( Vector& v,

            const typename Vector::IndexType i,

            const typename Vector::RealType& value )

{
@@ -32,7 +34,9 @@ void VectorOperations< Devices::MIC >::addElement( Vector& v, 
}



template< typename Vector >

void VectorOperations< Devices::MIC >::addElement( Vector& v,

void

VectorOperations< Devices::MIC >::

addElement( Vector& v,

            const typename Vector::IndexType i,

            const typename Vector::RealType& value,

            const typename Vector::RealType& thisElementMultiplicator )
@@ -43,7 +47,9 @@ void VectorOperations< Devices::MIC >::addElement( Vector& v, 
}



template< typename Vector >

typename Vector::RealType VectorOperations< Devices::MIC >::getVectorMax( const Vector& v )

typename Vector::RealType

VectorOperations< Devices::MIC >::

getVectorMax( const Vector& v )

{

   //tady je možnost paralelizace

   typename Vector::RealType result;
@@ -64,7 +70,9 @@ typename Vector::RealType VectorOperations< Devices::MIC >::getVectorMax( const 
}



template< typename Vector >

typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorMin( const Vector& v )

typename Vector::RealType

VectorOperations< Devices::MIC >::

getVectorMin( const Vector& v )

{

   //tady je možnost paralelizace

   typename Vector::RealType result;
@@ -85,7 +93,9 @@ typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorMin( co 
}



template< typename Vector >

typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorAbsMax( const Vector& v )

typename Vector::RealType

VectorOperations< Devices::MIC >::

getVectorAbsMax( const Vector& v )

{

   //tady je možnost paralelizace

   typename Vector::RealType result;
@@ -95,19 +105,20 @@ typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorAbsMax( 


   #pragma offload target(mic) in(vct,size) out(result)

   {

       result=std::fabs(vct.pointer[0]);

      result=TNL::abs(vct.pointer[0]);

      for(typename Vector::IndexType i=1;i<size;i++)

      {

           if(result<std::fabs(vct.pointer[i]))

               result=std::fabs(vct.pointer[i]);

         if(result<TNL::abs(vct.pointer[i]))

            result=TNL::abs(vct.pointer[i]);

      }

   }

   return result;

}





template< typename Vector >

typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorAbsMin( const Vector& v )

typename Vector::RealType

VectorOperations< Devices::MIC >::

getVectorAbsMin( const Vector& v )

{

   //tady je možnost paralelizace

   typename Vector::RealType result;
@@ -117,11 +128,11 @@ typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorAbsMin( 


   #pragma offload target(mic) in(vct,size) out(result)

   {

       result=std::fabs(vct.pointer[0]);

      result=TNL::abs(vct.pointer[0]);

      for(typename Vector::IndexType i=1;i<size;i++)

      {

           if(result>std::fabs(vct.pointer[i]))

               result=std::fabs(vct.pointer[i]);

         if(result>TNL::abs(vct.pointer[i]))

            result=TNL::abs(vct.pointer[i]);

      }

   }

   return result;
@@ -129,7 +140,8 @@ typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorAbsMin( 


template< typename Vector >

typename Vector::RealType

VectorOperations< Devices::MIC >::getVectorL1Norm( const Vector& v )

VectorOperations< Devices::MIC >::

getVectorL1Norm( const Vector& v )

{

   typedef typename Vector::RealType Real;

   typedef typename Vector::IndexType Index;
@@ -144,14 +156,15 @@ VectorOperations< Devices::MIC >::getVectorL1Norm( const Vector& v ) 
   {

      #pragma omp parallel for reduction(+:result)// if( n > OpenMPVectorOperationsThreshold ) // TODO: check this threshold

      for( Index i = 0; i < n; i ++ )

            result += std::fabs( vct.pointer[ i ] );

         result += TNL::abs( vct.pointer[ i ] );

   }

   return result;

}



template< typename Vector >

typename Vector::RealType

VectorOperations< Devices::MIC >::getVectorL2Norm( const Vector& v )

VectorOperations< Devices::MIC >::

getVectorL2Norm( const Vector& v )

{

   typedef typename Vector::RealType Real;

   typedef typename Vector::IndexType Index;
@@ -170,12 +183,13 @@ VectorOperations< Devices::MIC >::getVectorL2Norm( const Vector& v ) 
         result += aux * aux;

      }

   }

   return std::sqrt( result );

   return TNL::sqrt( result );

}



template< typename Vector >

typename Vector::RealType

VectorOperations< Devices::MIC >:: getVectorLpNorm( const Vector& v,

VectorOperations< Devices::MIC >::

getVectorLpNorm( const Vector& v,

                 const typename Vector::RealType& p )

{

   typedef typename Vector::RealType Real;
@@ -198,16 +212,17 @@ VectorOperations< Devices::MIC >:: getVectorLpNorm( const Vector& v, 
      #pragma omp parallel for reduction(+:result) //if( n > OpenMPVectorOperationsThreshold ) // TODO: check this threshold

      for( Index i = 0; i < n; i ++ )

      {

      result += pow( std::fabs( vct.pointer[ i ] ), p );

         result += TNL::pow( TNL::abs( vct.pointer[ i ] ), p );

      }

   }

   return pow( result, 1.0 / p );

   return TNL::pow( result, 1.0 / p );

}



template< typename Vector >

typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorSum( const Vector& v )

typename Vector::RealType

VectorOperations< Devices::MIC >::

getVectorSum( const Vector& v )

{

 

   typedef typename Vector::RealType Real;

   typedef typename Vector::IndexType Index;

   TNL_ASSERT( v. getSize() > 0, );
@@ -226,9 +241,10 @@ typename Vector :: RealType VectorOperations< Devices::MIC > :: getVectorSum( co 
   return result;

}





template< typename Vector1, typename Vector2 >

typename Vector1 :: RealType VectorOperations< Devices::MIC> :: getVectorDifferenceMax( const Vector1& v1,

typename Vector1::RealType

VectorOperations< Devices::MIC>::

getVectorDifferenceMax( const Vector1& v1,

                        const Vector2& v2 )

{

   typedef typename Vector1::RealType Real;
@@ -247,13 +263,15 @@ typename Vector1 :: RealType VectorOperations< Devices::MIC> :: getVectorDiffere 
   {

      result = vct1.pointer[0] - vct2.pointer[0];

      for( Index i = 1; i < n; i ++ )

      result = max( result, vct1.pointer[ i ] - vct2.pointer[ i ] );

         result = TNL::max( result, vct1.pointer[ i ] - vct2.pointer[ i ] );

   }

   return result;

}



template< typename Vector1, typename Vector2 >

typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getVectorDifferenceMin( const Vector1& v1,

typename Vector1::RealType

VectorOperations< Devices::MIC >::

getVectorDifferenceMin( const Vector1& v1,

                        const Vector2& v2 )

{

   typedef typename Vector1::RealType Real;
@@ -272,13 +290,15 @@ typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getVectorDiffer 
   {

      result = vct1.pointer[0] - vct2.pointer[0];

      for( Index i = 1; i < n; i ++ )

      result = min( result, vct1.pointer[ i ] - vct2.pointer[ i ] );

         result = TNL::min( result, vct1.pointer[ i ] - vct2.pointer[ i ] );

   }

   return result;

}



template< typename Vector1, typename Vector2 >

typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getVectorDifferenceAbsMax( const Vector1& v1,

typename Vector1::RealType

VectorOperations< Devices::MIC >::

getVectorDifferenceAbsMax( const Vector1& v1,

                           const Vector2& v2 )

{

   typedef typename Vector1::RealType Real;
@@ -295,15 +315,17 @@ typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getVectorDiffer 


   #pragma offload target(mic) in(n,vct1,vct2) out(result)

   {

   result = std::fabs(vct1.pointer[0] - vct2.pointer[0]);

      result = TNL::abs(vct1.pointer[0] - vct2.pointer[0]);

      for( Index i = 1; i < n; i ++ )

      result = max( result, std::fabs(vct1.pointer[ i ] - vct2.pointer[ i ]) );

         result = TNL::max( result, TNL::abs(vct1.pointer[ i ] - vct2.pointer[ i ]) );

   }

   return result;

}



template< typename Vector1, typename Vector2 >

typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getVectorDifferenceAbsMin( const Vector1& v1,

typename Vector1::RealType

VectorOperations< Devices::MIC >::

getVectorDifferenceAbsMin( const Vector1& v1,

                           const Vector2& v2 )

{

   typedef typename Vector1::RealType Real;
@@ -320,9 +342,9 @@ typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getVectorDiffer 


   #pragma offload target(mic) in(n,vct1,vct2) out(result)

   {

   result = std::fabs(vct1.pointer[0] - vct2.pointer[0]);

      result = TNL::abs(vct1.pointer[0] - vct2.pointer[0]);

      for( Index i = 1; i < n; i ++ )

      result = min( result, std::fabs(vct1.pointer[ i ] - vct2.pointer[ i ]) );

         result = TNL::min( result, TNL::abs(vct1.pointer[ i ] - vct2.pointer[ i ]) );

   }

   return result;

}
@@ -349,10 +371,9 @@ getVectorDifferenceL1Norm( const Vector1& v1, 
   #pragma offload target(mic) in(n,vct1,vct2) inout(result)

   {

      for( Index i = 0; i < n; i ++ )

      result += std::fabs( vct1.pointer[ i ] - vct2.pointer[ i ] );

         result += TNL::abs( vct1.pointer[ i ] - vct2.pointer[ i ] );

   }

   return result;

   

}



template< typename Vector1, typename Vector2 >
@@ -361,7 +382,6 @@ VectorOperations< Devices::MIC >:: 
getVectorDifferenceL2Norm( const Vector1& v1,

                           const Vector2& v2 )

{

    

   typedef typename Vector1::RealType Real;

   typedef typename Vector1::IndexType Index;
@@ -380,13 +400,12 @@ getVectorDifferenceL2Norm( const Vector1& v1, 
   {

      for( Index i = 0; i < n; i ++ )

      {

      Real aux = std::fabs( vct1.pointer[ i ] - vct2.pointer[ i ] );

         Real aux = TNL::abs( vct1.pointer[ i ] - vct2.pointer[ i ] );

         result += aux * aux;

      }

   }



 return std::sqrt( result );

   

   return TNL::sqrt( result );

}
@@ -397,7 +416,6 @@ getVectorDifferenceLpNorm( const Vector1& v1, 
                           const Vector2& v2,

                           const typename Vector1::RealType& p )

{

    

   typedef typename Vector1::RealType Real;

   typedef typename Vector1::IndexType Index;
@@ -423,15 +441,16 @@ getVectorDifferenceLpNorm( const Vector1& v1, 
   {

      for( Index i = 0; i < n; i ++ )

      {

      result += pow( std::fabs( vct1.pointer[ i ] - vct2.pointer[ i ] ), p );

         result += TNL::pow( TNL::abs( vct1.pointer[ i ] - vct2.pointer[ i ] ), p );

      }

   }

   return pow( result, 1.0 / p );

 

   return TNL::pow( result, 1.0 / p );

}



template< typename Vector1, typename Vector2 >

typename Vector1::RealType VectorOperations< Devices::MIC > :: getVectorDifferenceSum( const Vector1& v1,

typename Vector1::RealType

VectorOperations< Devices::MIC >::

getVectorDifferenceSum( const Vector1& v1,

                        const Vector2& v2 )

{

   typedef typename Vector1::RealType Real;
@@ -457,7 +476,9 @@ typename Vector1::RealType VectorOperations< Devices::MIC > :: getVectorDifferen 




template< typename Vector >

void VectorOperations< Devices::MIC > :: vectorScalarMultiplication( Vector& v,

void

VectorOperations< Devices::MIC >::

vectorScalarMultiplication( Vector& v,

                            const typename Vector::RealType& alpha )

{

   typedef typename Vector::RealType Real;
@@ -475,12 +496,13 @@ void VectorOperations< Devices::MIC > :: vectorScalarMultiplication( Vector& v, 
      for( Index i = 0; i < n; i ++ )

         vct.pointer[ i ] *= a;

   }



}





template< typename Vector1, typename Vector2 >

typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getScalarProduct( const Vector1& v1,

typename Vector1::RealType

VectorOperations< Devices::MIC >::

getScalarProduct( const Vector1& v1,

                  const Vector2& v2 )

{

   typedef typename Vector1::RealType Real;
@@ -497,7 +519,6 @@ typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getScalarProduc 


   #pragma offload target(mic) in(vct1,vct2,n) inout(result)

   {



      #pragma omp parallel for reduction(+:result)// if( n > OpenMPVectorOperationsThreshold ) // TODO: check this threshold

      for( Index i = 0; i < n; i++ )

         result += vct1.pointer[ i ] * vct2.pointer[ i ];
@@ -524,7 +545,9 @@ typename Vector1 :: RealType VectorOperations< Devices::MIC > :: getScalarProduc 
}



template< typename Vector1, typename Vector2 >

void VectorOperations< Devices::MIC > :: addVector( Vector1& y,

void

VectorOperations< Devices::MIC >::

addVector( Vector1& y,

           const Vector2& x,

           const typename Vector2::RealType& alpha,

           const typename Vector1::RealType& thisMultiplicator )
@@ -547,7 +570,6 @@ void VectorOperations< Devices::MIC > :: addVector( Vector1& y, 
      for( Index i = 0; i < n; i ++ )

         vct.pointer[ i ] = t * vct.pointer[ i ] + a * vct2.pointer[ i ];

   }

    

}



template< typename Vector1,
@@ -584,15 +606,15 @@ addVectors( Vector1& v, 
      for( Index i = 0; i < n; i ++ )

         vct.pointer[ i ] = t * vct.pointer[ i ] + m1 * vct1.pointer[ i ] + m2 * vct2.pointer[ i ];

   }

    

}



template< typename Vector >

void VectorOperations< Devices::MIC >::computePrefixSum( Vector& v,

void

VectorOperations< Devices::MIC >::

computePrefixSum( Vector& v,

                  typename Vector::IndexType begin,

                  typename Vector::IndexType end )

{

    

   typedef typename Vector::IndexType Index;



   //std::cout << v.getSize()<< "    " << end <<endl;
@@ -609,15 +631,15 @@ void VectorOperations< Devices::MIC >::computePrefixSum( Vector& v, 
      for( Index i = begin + 1; i < end; i++ )

         vct.pointer[ i ] += vct.pointer[ i - 1 ];

   }

    

}



template< typename Vector >

void VectorOperations< Devices::MIC >::computeExclusivePrefixSum( Vector& v,

void

VectorOperations< Devices::MIC >::

computeExclusivePrefixSum( Vector& v,

                           typename Vector::IndexType begin,

                           typename Vector::IndexType end )

{

   

   typedef typename Vector::IndexType Index;

   typedef typename Vector::RealType Real;

   TNL_ASSERT( v.getSize() > 0, );
@@ -631,7 +653,6 @@ void VectorOperations< Devices::MIC >::computeExclusivePrefixSum( Vector& v, 


   #pragma offload target(mic) in(vct,begin,end)

   {

   

      Real aux( vct.pointer[ begin ] );

      vct.pointer[ begin ] = 0.0;

      for( Index i = begin + 1; i < end; i++ )
@@ -643,7 +664,6 @@ void VectorOperations< Devices::MIC >::computeExclusivePrefixSum( Vector& v, 
   }

}





} // namespace Algorithms

} // namespace Containers

} // namespace TNL