Cross product and triple scalar product for all vectors (dense, tiny, engine)

flens/blas/closures/auxiliary/result.h

@@ -53,10 +53,13 @@ struct Result
 template <typename Op, typename L, typename R>
 struct Result<VectorClosure<Op, L, R> >
 {
-    typedef typename VectorClosure<Op, L, R>::ElementType T;
+    typedef typename Result<L>::Type  LT; // result type of lhs
+    typedef typename Result<R>::Type  RT; // result type of rhs
 
-    typedef DenseVector<Array<T> >  Type;
-    typedef typename Type::NoView   NoView;
+    typedef typename std::remove_reference<LT>::type  LTT; // remove if necesary reference from lhs result type
+
+    typedef typename std::conditional<!IsVectorClosure<LTT>::value && IsVector<LTT>::value, LT, RT>::type  Type;   // the vector type involved in the closure
+    typedef typename NoView<Type>::Type                                                                    NoView;
 };
 
 //-- MatrixClosures ------------------------------------------------------------

flens/blas/closures/level1/axpy.h

@@ -90,6 +90,15 @@ template <typename ALPHA, typename VX, typename SV, typename VY>
     axpy(const ALPHA &alpha,
          const VectorClosure<OpDiv, VX, SV> &x, Vector<VY> &y);
 
+// y += x1 % x2
+template <typename ALPHA, typename VL, typename VR, typename VY>
+    typename RestrictTo<VCDefaultEval<OpCross, VL, VR>::value
+                     && IsVector<VL>::value
+                     && IsVector<VR>::value,
+             void>::Type
+    axpy(const ALPHA &alpha,
+         const VectorClosure<OpCross, VL, VR> &x, Vector<VY> &y);
+
 // y += A*x
 template <typename ALPHA, typename ML, typename VR, typename VY>
     typename RestrictTo<VCDefaultEval<OpMult, ML, VR>::value

flens/blas/closures/level1/axpy.tcc

@@ -429,6 +429,32 @@ axpy(const ALPHA &alpha,
     FLENS_BLASLOG_END;
 }
 
+//------------------------------------------------------------------------------
+//
+// y += x1 % x2
+//
+template <typename ALPHA, typename VL, typename VR, typename VY>
+typename RestrictTo<VCDefaultEval<OpCross, VL, VR>::value
+                 && IsVector<VL>::value
+                 && IsVector<VR>::value,
+         void>::Type
+axpy(const ALPHA &alpha, const VectorClosure<OpCross, VL, VR> &x, Vector<VY> &y)
+{
+    FLENS_BLASLOG_BEGIN_AXPY(alpha, x, y);
+//
+// Computes the result of the cross closure x first and stores it in vx.
+//
+    typedef VectorClosure<OpCross, VL, VR> VC;
+    const typename Result<VC>::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

flens/blas/closures/level1/copy.h

@@ -83,6 +83,15 @@ template <typename VX, typename VY>
              void>::Type
     copy(const VectorClosureOpConj<VX> &x, Vector<VY> &y);
 
+// y = x1 % x2
+template <typename VL, typename VR, typename VY>
+    typename RestrictTo<VCDefaultEval<OpCross, VL, VR>::value
+                     && IsVector<VL>::value
+                     && IsVector<VR>::value,
+             void>::Type
+    copy(const VectorClosure<OpCross, VL, VR> &x, Vector<VY> &y);
+
+
 // y = A*x
 template <typename ML, typename VR, typename VY>
     typename RestrictTo<VCDefaultEval<OpMult, ML, VR>::value
@@ -100,6 +109,7 @@ template <typename VL, typename MR, typename VY>
     copy(const VectorClosure<OpMult, VL, MR> &xA, Vector<VY> &y);
 
 
+
 //
 //  This gets called if everything else fails
 //

flens/blas/closures/level1/copy.tcc

@@ -145,6 +145,23 @@ copy(const VectorClosureOpConj<VX> &x, Vector<VY> &y)
     FLENS_BLASLOG_END;
 }
 
+//------------------------------------------------------------------------------
+//
+//  y = x1 % x2
+//
+
+template <typename VL, typename VR, typename VY>
+typename RestrictTo<VCDefaultEval<OpCross, VL, VR>::value
+                 && IsVector<VL>::value
+                 && IsVector<VR>::value,
+         void>::Type
+copy(const VectorClosure<OpCross, VL, VR> &x, Vector<VY> &y)
+{
+    FLENS_BLASLOG_BEGIN_COPY(x1x2, y);
+    copyCross(x.left(), x.right(), y.impl());
+    FLENS_BLASLOG_END;
+}
+
 //------------------------------------------------------------------------------
 //
 //  y = A*x

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 <cxxblas/cxxblas.h>
+#include <flens/blas/closures/tweaks/defaulteval.h>
+#include <flens/blas/operators/operators.h>
+#include <flens/matrixtypes/matrixtypes.h>
+#include <flens/typedefs.h>
+#include <flens/vectortypes/vectortypes.h>
+
+namespace flens { namespace blas {
+
+//------------------------------------------------------------------------------
+//
+//  y = x1 % x2
+//
+template <typename VX1, typename VX2, typename VY>
+    void
+    copyCross(const VX1 &x1, VX2 &x2, VY &y);
+
+} } // namespace blas, flens
+
+#endif // FLENS_BLAS_CLOSURES_LEVEL1_COPYCROSS_H

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 <flens/auxiliary/auxiliary.h>
+#include <flens/blas/closures/closures.h>
+#include <flens/blas/level1/level1.h>
+#include <flens/blas/level2/level2.h>
+#include <flens/blas/level3/level3.h>
+#include <flens/typedefs.h>
+
+#ifdef FLENS_DEBUG_CLOSURES
+#   include <flens/blas/blaslogon.h>
+#else
+#   include <flens/blas/blaslogoff.h>
+#endif
+
+namespace flens { namespace blas {
+
+//------------------------------------------------------------------------------
+//
+//  y = x1 % x2
+//
+template <typename VX1, typename VX2, typename VY>
+void
+copyCross(const VX1 &x1, VX2 &x2, VY &y)
+{
+    // Temporary are maybe created but it is not so bad ;-)
+    const typename Result<VX1>::Type& lhs = x1;
+    const typename Result<VX2>::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

flens/blas/closures/level1/level1.h

@@ -36,6 +36,7 @@
 #include <flens/blas/closures/level1/axpy.h>
 #include <flens/blas/closures/level1/copy.h>
 #include <flens/blas/closures/level1/copyconj.h>
+#include <flens/blas/closures/level1/copycross.h>
 #include <flens/blas/closures/level1/copyrscal.h>
 #include <flens/blas/closures/level1/copyscal.h>
 #include <flens/blas/closures/level1/copysum.h>

flens/blas/closures/level1/level1.tcc

@@ -36,6 +36,7 @@
 #include <flens/blas/closures/level1/axpy.tcc>
 #include <flens/blas/closures/level1/copy.tcc>
 #include <flens/blas/closures/level1/copyconj.tcc>
+#include <flens/blas/closures/level1/copycross.tcc>
 #include <flens/blas/closures/level1/copyrscal.tcc>
 #include <flens/blas/closures/level1/copyscal.tcc>
 #include <flens/blas/closures/level1/copysum.tcc>

flens/blas/level1/asum.h

@@ -48,6 +48,14 @@ template <typename X>
     const typename ComplexTrait<typename X::ElementType>::PrimitiveType
     asum(const DenseVector<X> &x);
 
+template <typename X, typename T>
+    typename RestrictTo<IsNotComplex<T>::value, void>::Type
+    asum(const TinyVector<X> &x, T &absoluteSum);
+
+template <typename X>
+    const typename ComplexTrait<typename X::ElementType>::PrimitiveType
+    asum(const TinyVector<X> &x);
+
 //-- BLAS Level 1 extensions ---------------------------------------------------
 
 //== GeneralMatrix

flens/blas/level1/asum.tcc

@@ -65,6 +65,30 @@ asum(const DenseVector<X> &x)
     return absoluteSum;
 }
 
+template <typename X, typename T>
+typename RestrictTo<IsNotComplex<T>::value, void>::Type
+asum(const TinyVector<X> &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 <typename X>
+const typename ComplexTrait<typename X::ElementType>::PrimitiveType
+asum(const TinyVector<X> &x)
+{
+    typename ComplexTrait<typename X::ElementType>::PrimitiveType  absoluteSum;
+
+    asum(x, absoluteSum);
+    return absoluteSum;
+}
+
 //-- BLAS Level 1 extensions ---------------------------------------------------
 
 //== GeneralMatrix

flens/blas/level1/dot.h

@@ -64,6 +64,35 @@ template <typename X, typename Y>
                             typename Y::ElementType>::Type
     dotu(const DenseVector<X> &x, const DenseVector<Y> &y);
 
+// TinyVector version
+
+template <typename X, typename Y, typename T>
+    void
+    dot(const TinyVector<X> &x, const TinyVector<Y> &y, T &result);
+
+template <typename X, typename Y, typename T>
+    void
+    dotc(const TinyVector<X> &x, const TinyVector<Y> &y, T &result);
+
+template <typename X, typename Y, typename T>
+    void
+    dotu(const TinyVector<X> &x, const TinyVector<Y> &y, T &result);
+
+template <typename X, typename Y>
+    typename CompatibleType<typename X::ElementType,
+                            typename Y::ElementType>::Type
+    dot(const TinyVector<X> &x, const TinyVector<Y> &y);
+
+template <typename X, typename Y>
+    typename CompatibleType<typename X::ElementType,
+                            typename Y::ElementType>::Type
+    dotc(const TinyVector<X> &x, const TinyVector<Y> &y);
+
+template <typename X, typename Y>
+    typename CompatibleType<typename X::ElementType,
+                            typename Y::ElementType>::Type
+    dotu(const TinyVector<X> &x, const TinyVector<Y> &y);
+
 } } // namespace blas, flens
 
 #endif // FLENS_BLAS_LEVEL1_DOT_H