diff --git a/flens/blas/closures/auxiliary/result.h b/flens/blas/closures/auxiliary/result.h index bfc6eb7c..b6dfd9ca 100644 --- a/flens/blas/closures/auxiliary/result.h +++ b/flens/blas/closures/auxiliary/result.h @@ -53,10 +53,13 @@ struct Result template struct Result > { - typedef typename VectorClosure::ElementType T; + typedef typename Result::Type LT; // result type of lhs + typedef typename Result::Type RT; // result type of rhs - typedef DenseVector > Type; - typedef typename Type::NoView NoView; + typedef typename std::remove_reference::type LTT; // remove if necesary reference from lhs result type + + typedef typename std::conditional::value && IsVector::value, LT, RT>::type Type; // the vector type involved in the closure + typedef typename NoView::Type NoView; }; //-- MatrixClosures ------------------------------------------------------------ diff --git a/flens/blas/closures/level1/axpy.h b/flens/blas/closures/level1/axpy.h index e1220a2e..8ee212a4 100644 --- a/flens/blas/closures/level1/axpy.h +++ b/flens/blas/closures/level1/axpy.h @@ -90,6 +90,15 @@ template axpy(const ALPHA &alpha, const VectorClosure &x, Vector &y); +// y += x1 % x2 +template + typename RestrictTo::value + && IsVector::value + && IsVector::value, + void>::Type + axpy(const ALPHA &alpha, + const VectorClosure &x, Vector &y); + // y += A*x template typename RestrictTo::value diff --git a/flens/blas/closures/level1/axpy.tcc b/flens/blas/closures/level1/axpy.tcc index 8cb34197..93235de6 100644 --- a/flens/blas/closures/level1/axpy.tcc +++ b/flens/blas/closures/level1/axpy.tcc @@ -429,6 +429,32 @@ axpy(const ALPHA &alpha, FLENS_BLASLOG_END; } +//------------------------------------------------------------------------------ +// +// y += x1 % x2 +// +template +typename RestrictTo::value + && IsVector::value + && IsVector::value, + void>::Type +axpy(const ALPHA &alpha, const VectorClosure &x, Vector &y) +{ + FLENS_BLASLOG_BEGIN_AXPY(alpha, x, y); +// +// Computes the result of the cross closure x first and stores it in vx. +// + typedef VectorClosure VC; + const typename Result::Type vx = x; + +// +// Update y with vx, i.e. compute y = y + alpha*vx +// + axpy(alpha, vx.impl(), y.impl()); + + FLENS_BLASLOG_END; +} + //------------------------------------------------------------------------------ // // y += x*A diff --git a/flens/blas/closures/level1/copy.h b/flens/blas/closures/level1/copy.h index 9a15bc18..451d2f3c 100644 --- a/flens/blas/closures/level1/copy.h +++ b/flens/blas/closures/level1/copy.h @@ -83,6 +83,15 @@ template void>::Type copy(const VectorClosureOpConj &x, Vector &y); +// y = x1 % x2 +template + typename RestrictTo::value + && IsVector::value + && IsVector::value, + void>::Type + copy(const VectorClosure &x, Vector &y); + + // y = A*x template typename RestrictTo::value @@ -100,6 +109,7 @@ template copy(const VectorClosure &xA, Vector &y); + // // This gets called if everything else fails // diff --git a/flens/blas/closures/level1/copy.tcc b/flens/blas/closures/level1/copy.tcc index 516304ab..6acf19de 100644 --- a/flens/blas/closures/level1/copy.tcc +++ b/flens/blas/closures/level1/copy.tcc @@ -145,6 +145,23 @@ copy(const VectorClosureOpConj &x, Vector &y) FLENS_BLASLOG_END; } +//------------------------------------------------------------------------------ +// +// y = x1 % x2 +// + +template +typename RestrictTo::value + && IsVector::value + && IsVector::value, + void>::Type +copy(const VectorClosure &x, Vector &y) +{ + FLENS_BLASLOG_BEGIN_COPY(x1x2, y); + copyCross(x.left(), x.right(), y.impl()); + FLENS_BLASLOG_END; +} + //------------------------------------------------------------------------------ // // y = A*x diff --git a/flens/blas/closures/level1/copycross.h b/flens/blas/closures/level1/copycross.h new file mode 100644 index 00000000..0839739d --- /dev/null +++ b/flens/blas/closures/level1/copycross.h @@ -0,0 +1,55 @@ +/* + * Copyright (c) 2016, Thibaud Kloczko + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1) Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2) Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * 3) Neither the name of the FLENS development group nor the names of + * its contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef FLENS_BLAS_CLOSURES_LEVEL1_COPYCROSS_H +#define FLENS_BLAS_CLOSURES_LEVEL1_COPYCROSS_H 1 + +#include +#include +#include +#include +#include +#include + +namespace flens { namespace blas { + +//------------------------------------------------------------------------------ +// +// y = x1 % x2 +// +template + void + copyCross(const VX1 &x1, VX2 &x2, VY &y); + +} } // namespace blas, flens + +#endif // FLENS_BLAS_CLOSURES_LEVEL1_COPYCROSS_H diff --git a/flens/blas/closures/level1/copycross.tcc b/flens/blas/closures/level1/copycross.tcc new file mode 100644 index 00000000..5b2667e7 --- /dev/null +++ b/flens/blas/closures/level1/copycross.tcc @@ -0,0 +1,70 @@ +/* + * Copyright (c) 2016, Thibaud Kloczko + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1) Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2) Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * 3) Neither the name of the FLENS development group nor the names of + * its contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef FLENS_BLAS_CLOSURES_LEVEL1_COPYCROSS_TCC +#define FLENS_BLAS_CLOSURES_LEVEL1_COPYCROSS_TCC 1 + +#include +#include +#include +#include +#include +#include + +#ifdef FLENS_DEBUG_CLOSURES +# include +#else +# include +#endif + +namespace flens { namespace blas { + +//------------------------------------------------------------------------------ +// +// y = x1 % x2 +// +template +void +copyCross(const VX1 &x1, VX2 &x2, VY &y) +{ + // Temporary are maybe created but it is not so bad ;-) + const typename Result::Type& lhs = x1; + const typename Result::Type& rhs = x2; + + y(1) = lhs(2) * rhs(3) - lhs(3) * rhs(2); + y(2) = lhs(3) * rhs(1) - lhs(1) * rhs(3); + y(3) = lhs(1) * rhs(2) - lhs(2) * rhs(1); +} + +} } // namespace blas, flens + +#endif // FLENS_BLAS_CLOSURES_LEVEL1_COPYCROSS_TCC diff --git a/flens/blas/closures/level1/level1.h b/flens/blas/closures/level1/level1.h index 656c8f82..9daf7426 100644 --- a/flens/blas/closures/level1/level1.h +++ b/flens/blas/closures/level1/level1.h @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include diff --git a/flens/blas/closures/level1/level1.tcc b/flens/blas/closures/level1/level1.tcc index ed7141af..ba1355b9 100644 --- a/flens/blas/closures/level1/level1.tcc +++ b/flens/blas/closures/level1/level1.tcc @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include diff --git a/flens/blas/level1/asum.h b/flens/blas/level1/asum.h index a0d8c068..e014d5db 100644 --- a/flens/blas/level1/asum.h +++ b/flens/blas/level1/asum.h @@ -48,6 +48,14 @@ template const typename ComplexTrait::PrimitiveType asum(const DenseVector &x); +template + typename RestrictTo::value, void>::Type + asum(const TinyVector &x, T &absoluteSum); + +template + const typename ComplexTrait::PrimitiveType + asum(const TinyVector &x); + //-- BLAS Level 1 extensions --------------------------------------------------- //== GeneralMatrix diff --git a/flens/blas/level1/asum.tcc b/flens/blas/level1/asum.tcc index 262e60e4..4bc386a4 100644 --- a/flens/blas/level1/asum.tcc +++ b/flens/blas/level1/asum.tcc @@ -65,6 +65,30 @@ asum(const DenseVector &x) return absoluteSum; } +template +typename RestrictTo::value, void>::Type +asum(const TinyVector &x, T &absoluteSum) +{ +# ifdef HAVE_CXXBLAS_ASUM + const int length = X::length; + const int stride = X::stride; + + cxxblas::asum(length, x.data(), stride, absoluteSum); +# else + ASSERT(0); +# endif +} + +template +const typename ComplexTrait::PrimitiveType +asum(const TinyVector &x) +{ + typename ComplexTrait::PrimitiveType absoluteSum; + + asum(x, absoluteSum); + return absoluteSum; +} + //-- BLAS Level 1 extensions --------------------------------------------------- //== GeneralMatrix diff --git a/flens/blas/level1/dot.h b/flens/blas/level1/dot.h index 92b7d09a..d5e57664 100644 --- a/flens/blas/level1/dot.h +++ b/flens/blas/level1/dot.h @@ -64,6 +64,35 @@ template typename Y::ElementType>::Type dotu(const DenseVector &x, const DenseVector &y); +// TinyVector version + +template + void + dot(const TinyVector &x, const TinyVector &y, T &result); + +template + void + dotc(const TinyVector &x, const TinyVector &y, T &result); + +template + void + dotu(const TinyVector &x, const TinyVector &y, T &result); + +template + typename CompatibleType::Type + dot(const TinyVector &x, const TinyVector &y); + +template + typename CompatibleType::Type + dotc(const TinyVector &x, const TinyVector &y); + +template + typename CompatibleType::Type + dotu(const TinyVector &x, const TinyVector &y); + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_DOT_H diff --git a/flens/blas/level1/dot.tcc b/flens/blas/level1/dot.tcc index fa000aea..5fa53a5c 100644 --- a/flens/blas/level1/dot.tcc +++ b/flens/blas/level1/dot.tcc @@ -127,6 +127,103 @@ dotu(const DenseVector &x, const DenseVector &y) return val; } +// -- TinyVector version + +template +void +dot(const TinyVector &x, const TinyVector &y, T &result) +{ + FLENS_BLASLOG_SETTAG("--> "); + FLENS_BLASLOG_BEGIN_DOT(x, y); + + typedef typename X::ElementType TX; + typedef typename Y::ElementType TY; + + const int lengthX = X::length; + const int strideX = X::stride; + const int lengthY = Y::length; + const int strideY = Y::stride; + + ASSERT(lengthX == lengthY); + +# ifdef HAVE_CXXBLAS_DOT + cxxblas::dot(lengthX, + x.data(), strideX, + y.data(), strideY, result); +# else + ASSERT(0); +# endif + + FLENS_BLASLOG_END; + FLENS_BLASLOG_UNSETTAG; +} + +template +void +dotc(const TinyVector &x, const TinyVector &y, T &result) +{ + dotc(x, y, result); +} + +template +void +dotu(const TinyVector &x, const TinyVector &y, T &result) +{ + FLENS_BLASLOG_SETTAG("--> "); + FLENS_BLASLOG_BEGIN_DOTU(x, y); + + typedef typename X::ElementType TX; + typedef typename Y::ElementType TY; + + const int lengthX = X::length; + const int strideX = X::stride; + const int lengthY = Y::length; + const int strideY = Y::stride; + + ASSERT(lengthX == lengthY); + +# ifdef HAVE_CXXBLAS_DOTU + cxxblas::dotu(lengthX, + x.data(), strideX, + y.data(), strideY, result); +# else + ASSERT(0); +# endif + + FLENS_BLASLOG_END; + FLENS_BLASLOG_UNSETTAG; +} + +template +typename CompatibleType::Type +dot(const TinyVector &x, const TinyVector &y) +{ + typedef typename X::ElementType TX; + typedef typename Y::ElementType TY; + + typename CompatibleType::Type val; + dot(x, y, val); + return val; +} + +template +typename CompatibleType::Type +dotc(const TinyVector &x, const TinyVector &y) +{ + return dot(x, y); +} + +template +typename CompatibleType::Type +dotu(const TinyVector &x, const TinyVector &y) +{ + typedef typename X::ElementType TX; + typedef typename Y::ElementType TY; + + typename CompatibleType::Type val; + dotu(x, y, val); + return val; +} } } // namespace blas, flens diff --git a/flens/blas/level1/iamax.h b/flens/blas/level1/iamax.h index cc8d517d..b7be80f7 100644 --- a/flens/blas/level1/iamax.h +++ b/flens/blas/level1/iamax.h @@ -42,6 +42,10 @@ template typename DenseVector::IndexType iamax(const DenseVector &x); +template + typename TinyVector::IndexType + iamax(const TinyVector &x); + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_IAMAX_H diff --git a/flens/blas/level1/iamax.tcc b/flens/blas/level1/iamax.tcc index af98c783..4655d624 100644 --- a/flens/blas/level1/iamax.tcc +++ b/flens/blas/level1/iamax.tcc @@ -61,6 +61,24 @@ iamax(const DenseVector &x) return i + x.indexBase(); } +template +typename TinyVector::IndexType +iamax(const TinyVector &x) +{ + typename TinyVector::IndexType i; + +# ifdef HAVE_CXXBLAS_IAMAX + const int length = X::length; + const int stride = X::stride; + + cxxblas::iamax(length, x.data(), stride, i); +# else + ASSERT(0); +# endif + const int idbase = X::firstIndex; + return i + idbase; +} + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_IAMAX_TCC diff --git a/flens/blas/level1/nrm2.h b/flens/blas/level1/nrm2.h index 7369508c..7c00aaf7 100644 --- a/flens/blas/level1/nrm2.h +++ b/flens/blas/level1/nrm2.h @@ -46,6 +46,14 @@ template typename ComplexTrait::PrimitiveType nrm2(const DenseVector &x); +template + typename RestrictTo::value, void>::Type + nrm2(const TinyVector &x, T &norm); + +template + typename ComplexTrait::PrimitiveType + nrm2(const TinyVector &x); + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_NRM2_H diff --git a/flens/blas/level1/nrm2.tcc b/flens/blas/level1/nrm2.tcc index d0958602..c65d119b 100644 --- a/flens/blas/level1/nrm2.tcc +++ b/flens/blas/level1/nrm2.tcc @@ -67,6 +67,29 @@ nrm2(const DenseVector &x) return norm; } +template +typename RestrictTo::value, void>::Type +nrm2(const TinyVector &x, T &norm) +{ +# ifdef HAVE_CXXBLAS_NRM2 + const int length = X::length; + const int stride = X::stride; + + cxxblas::nrm2(length, x.data(), stride, norm); +# else + ASSERT(0); +# endif +} + +template +typename ComplexTrait::PrimitiveType +nrm2(const TinyVector &x) +{ + typename ComplexTrait::PrimitiveType norm; + nrm2(x, norm); + return norm; +} + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_NRM2_TCC diff --git a/flens/blas/level1/raxpy.h b/flens/blas/level1/raxpy.h index ce4adfb9..f2aae1ce 100644 --- a/flens/blas/level1/raxpy.h +++ b/flens/blas/level1/raxpy.h @@ -47,6 +47,12 @@ template void>::Type raxpy(const ALPHA &alpha, const VX &x, VY &&y); +template + typename RestrictTo::value + && IsTinyVector::value, + void>::Type + raxpy(const ALPHA &alpha, const VX &x, VY &&y); + template typename RestrictTo::value && IsGeMatrix::value, diff --git a/flens/blas/level1/raxpy.tcc b/flens/blas/level1/raxpy.tcc index 5cd323e2..022e6fbb 100644 --- a/flens/blas/level1/raxpy.tcc +++ b/flens/blas/level1/raxpy.tcc @@ -80,6 +80,37 @@ raxpy(const ALPHA &alpha, const VX &x, VY &&y) FLENS_BLASLOG_UNSETTAG; } +//-- tinyvector raxpy +// +template +typename RestrictTo::value + && IsTinyVector::value, + void>::Type +raxpy(const ALPHA &alpha, const VX &x, VY &&y) +{ + FLENS_BLASLOG_SETTAG("--> "); + FLENS_BLASLOG_BEGIN_RAXPY(alpha, x, y); + + typedef typename VX::ElementType TX; + typedef typename RemoveRef::Type VectorY; + typedef typename VectorY::ElementType TY; + + const int lengthX = VX::Engine::length; + const int strideX = VX::Engine::stride; + const int lengthY = VectorY::Engine::length; + const int strideY = VectorY::Engine::stride; + + ASSERT(lengthY == lengthX && lengthY != 0); + +# ifdef HAVE_CXXBLAS_RAXPY + cxxblas::raxpy(lengthX, alpha, x.data(), strideX, y.data(), strideY); +# else + ASSERT(0); +# endif + FLENS_BLASLOG_END; + FLENS_BLASLOG_UNSETTAG; +} + //-- geraxpy // // B += A/alpha diff --git a/flens/blas/level1/rot.h b/flens/blas/level1/rot.h index 0d0f8863..889ac6e2 100644 --- a/flens/blas/level1/rot.h +++ b/flens/blas/level1/rot.h @@ -57,6 +57,12 @@ template void>::Type rot(VX &&x, VY &&y, const TC &c, const TS &s); +template + typename RestrictTo::value + && IsTinyVector::value, + void>::Type + rot(VX &&x, VY &&y, const TC &c, const TS &s); + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_ROT_H diff --git a/flens/blas/level1/rot.tcc b/flens/blas/level1/rot.tcc index 3cbe7175..3d74529c 100644 --- a/flens/blas/level1/rot.tcc +++ b/flens/blas/level1/rot.tcc @@ -82,5 +82,30 @@ rot(VX &&x, VY &&y, const TC &c, const TS &s) # endif } +//-- rot +template +typename RestrictTo::value + && IsTinyVector::value, + void>::Type +rot(VX &&x, VY &&y, const TC &c, const TS &s) +{ + typedef typename std::remove_reference::type VectorX; + typedef typename VectorX::ElementType TX; + typedef typename std::remove_reference::type VectorY; + typedef typename VectorY::ElementType TY; + + const int lengthX = VectorX::Engine::length; + const int strideX = VectorX::Engine::stride; + const int lengthY = VectorY::Engine::length; + const int strideY = VectorY::Engine::stride; + + ASSERT(lengthX == lengthY); +# ifdef HAVE_CXXBLAS_ROT + cxxblas::rot(lengthX, x.data(), strideX, y.data(), strideY, c, s); +# else + ASSERT(0); +# endif +} + } } // namespace blas, flens diff --git a/flens/blas/level1/rotm.h b/flens/blas/level1/rotm.h index 24365886..44afdb06 100644 --- a/flens/blas/level1/rotm.h +++ b/flens/blas/level1/rotm.h @@ -50,6 +50,18 @@ template void>::Type rotm(VX &&x, VY &&y, const DenseVector &p); +//-- rotmg +template + void + rotmg(T &d1, T &d2, T &b1, T &b2, TinyVector &p); + +//-- rotm +template + typename RestrictTo::value + && IsTinyVector::value, + void>::Type + rotm(VX &&x, VY &&y, const TinyVector &p); + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_ROT_H diff --git a/flens/blas/level1/rotm.tcc b/flens/blas/level1/rotm.tcc index 9544bd93..2b3c48f5 100644 --- a/flens/blas/level1/rotm.tcc +++ b/flens/blas/level1/rotm.tcc @@ -97,6 +97,66 @@ rotm(VX &&x, VY &&y, const DenseVector &p) # endif } +//-- rotmg +template +void +rotmg(T &d1, T &d2, T &b1, T &b2, TinyVector &p) +{ +# ifdef HAVE_CXXBLAS_ROTMG + +# ifndef NDEBUG + typedef typename TinyVector::IndexType IndexType; +# endif + + const int length = VP::length; + ASSERT(length == IndexType(5)); + + cxxblas::rotmg(d1, d2, b1, b2, p.data()); +# else + ASSERT(0); + (void)(d1); + (void)(d2); + (void)(b1); + (void)(b2); + (void)(p); +# endif +} + +//-- rotm +template +typename RestrictTo::value + && IsTinyVector::value, + void>::Type +rotm(VX &&x, VY &&y, const TinyVector &p) +{ +# ifdef HAVE_CXXBLAS_ROTM + typedef typename std::remove_reference::type VectorX; + typedef typename VectorX::ElementType TX; + typedef typename std::remove_reference::type VectorY; + typedef typename VectorY::ElementType TY; + + const int lengthX = VectorX::Engine::length; + const int strideX = VectorX::Engine::stride; + const int lengthY = VectorY::Engine::length; + const int strideY = VectorY::Engine::stride; + + const int length = VP::length; + + typedef typename VP::IndexType IndexType; + + ASSERT(length == IndexType(5)); + ASSERT(lengthX == lengthY); + + cxxblas::rotm(lengthX, x.data(), strideX, y.data(), strideY, p.data()); + +# else + ASSERT(0); + (void)(x); + (void)(y); + (void)(p); +# endif +} + } } // namespace blas, flens #endif // FLENS_BLAS_LEVEL1_ROT_TCC diff --git a/flens/blas/operators/opcross.h b/flens/blas/operators/opcross.h new file mode 100644 index 00000000..d3eb7fe0 --- /dev/null +++ b/flens/blas/operators/opcross.h @@ -0,0 +1,65 @@ +/* + * Copyright (c) 2016, Thibaud Kloczko + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1) Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2) Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * 3) Neither the name of the FLENS development group nor the names of + * its contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef FLENS_BLAS_OPERATORS_OPCROSS_H +#define FLENS_BLAS_OPERATORS_OPCROSS_H 1 + +#include + +namespace flens { + +struct OpCross {}; + +//-- cross product ------------------------------------------------------------- +// (x % y) +template + const VectorClosure + operator%(const Vector &x, const Vector &y); + +//-- triple vector products ---------------------------------------------------- +// x * (y % z) +template + typename Promotion::Type>::Type + operator*(const VX &x, const VectorClosure &yz); + +// (x % y) * z +template + typename Promotion::Type, + typename VZ::Impl::ElementType>::Type + operator*(const VectorClosure &xy, const VZ &z); + +} // namespace flens + +#endif // FLENS_BLAS_OPERATORS_OPCROSS_H diff --git a/flens/blas/operators/opcross.tcc b/flens/blas/operators/opcross.tcc new file mode 100644 index 00000000..01ad3900 --- /dev/null +++ b/flens/blas/operators/opcross.tcc @@ -0,0 +1,93 @@ +/* + * Copyright (c) 2016, Thibaud Kloczko + * + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1) Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2) Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * 3) Neither the name of the FLENS development group nor the names of + * its contributors may be used to endorse or promote products derived + * from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef FLENS_BLAS_OPERATORS_OPCROSS_TCC +#define FLENS_BLAS_OPERATORS_OPCROSS_TCC 1 + +#include + +namespace flens { + +//-- cross product ------------------------------------------------------------- +// (x % y) + +template +const VectorClosure +operator%(const Vector &x, const Vector &y) +{ + typedef VectorClosure VC; + return VC(x.impl(), y.impl()); +} + +//-- triple vector products ---------------------------------------------------- +// x * (y % z) +template + typename Promotion::Type>::Type + operator*(const VX &x, const VectorClosure &yz) +{ + typedef typename VX::Impl::ElementType TX; + typedef typename VY::Impl::ElementType TY; + typedef typename VZ::Impl::ElementType TZ; + + typedef typename Promotion::Type>::Type T; + + const typename Result::Type& _x = x; + const typename Result::Type& _y = yz.left(); + const typename Result::Type& _z = yz.right(); + + T res = _x(1) * (_y(2) * _z(3) - _y(3) * _z(2)); + res += _x(2) * (_y(3) * _z(1) - _y(1) * _z(3)); + res += _x(3) * (_y(1) * _z(2) - _y(2) * _z(1)); + + return res; +} + +// (x % y) * z +template + typename Promotion::Type, + typename VZ::Impl::ElementType>::Type + operator*(const VectorClosure &xy, const VZ &z) +{ + typedef VectorClosure VC; + + const typename Result::Type xy_ = xy; + const typename Result::Type& z_ = z; + + return (xy_ * z_); +} + +} // namespace flens + +#endif // FLENS_BLAS_OPERATORS_OPCROSS_TCC diff --git a/flens/blas/operators/operators.h b/flens/blas/operators/operators.h index 9dc5170a..aa4bef33 100644 --- a/flens/blas/operators/operators.h +++ b/flens/blas/operators/operators.h @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include diff --git a/flens/blas/operators/operators.tcc b/flens/blas/operators/operators.tcc index 2bd5d556..e47fe8e4 100644 --- a/flens/blas/operators/operators.tcc +++ b/flens/blas/operators/operators.tcc @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include