add Matrix expressions (#177)

* add Matrix expressions

* clang-format

* clang-format

* clang-format

* mv norm2(matrix) from Algebra to AlgebraExpr

Author: gratienj

_clang-format

@@ -2,7 +2,7 @@ Language: Cpp
 # BasedOnStyle:  WebKit
 AccessModifierOffset: -1
 AlignAfterOpenBracket: true
-AlignConsecutiveAssignments: None
+AlignConsecutiveAssignments: false
 AlignEscapedNewlines: DontAlign
 AlignOperands: false
 AlignTrailingComments: false
@@ -17,7 +17,7 @@ AlwaysBreakBeforeMultilineStrings: false
 AlwaysBreakTemplateDeclarations: No
 BinPackArguments: true
 BinPackParameters: true
-BreakBeforeBinaryOperators: None
+BreakBeforeBinaryOperators: false
 
 BreakBeforeBraces: Custom
 BraceWrapping:

src/core/alien/core/backend/IInternalLinearAlgebraExprT.h

@@ -83,6 +83,13 @@ class IInternalLinearAlgebraExpr
    */
   virtual Real norm2(const Vector& x) const = 0;
 
+  /*!
+   * \brief Compute L2 (Frobenous) norm of a matrix
+   * \param[in] x The matrix on which norm2 is computed
+   * \returns The norm2 of the matrix
+   */
+  virtual Real norm2(const Matrix& x) const = 0;
+
   /*!
    * \brief Compute a matrix vector product
    *
@@ -125,6 +132,23 @@ class IInternalLinearAlgebraExpr
    */
   virtual void copy(const Vector& x, Vector& r) const = 0;
 
+  /*!
+   * \brief Copy a matrix in another one
+   *
+   * \param[in] x The vector to copy
+   * \param[in,out] r The copied vector
+   */
+  virtual void copy(const Matrix& a, Matrix& r) const = 0;
+
+  /*!
+   * \brief Add a matrix to another one
+   *
+   * \param[in] a The matrix to copy
+   *
+   * \param[in,out] r The copied vector
+   */
+  virtual void add(const Matrix& a, Matrix& r) const = 0;
+
   /*!
    * \brief Compute the dot product of two vectors
    * \param[in] x The first vector
@@ -140,6 +164,13 @@ class IInternalLinearAlgebraExpr
    */
   virtual void scal(Real alpha, Vector& x) const = 0;
 
+  /*!
+   * \brief Scale a matrix by a factor
+   * \param[in] alpha The real value to scale with
+   * \param[in,out] x The vector to be scaled
+   */
+  virtual void scal(Real alpha, Matrix& a) const = 0;
+
   /*!
    * \brief Extract the diagonal of a matrix in a vector
    * \param[in] a The matrix to extract the diagonal

src/core/alien/core/backend/IInternalLinearAlgebraT.h

@@ -77,7 +77,7 @@ class IInternalLinearAlgebra
 
   /*!
    * \brief Compute L2 norm of a vector
-   * \param[in] x The vector on which norm0 is computed
+   * \param[in] x The vector on which norm2 is computed
    * \returns The norm2 of the vector
    */
   virtual Real norm2(const Vector& x) const = 0;

src/core/alien/core/backend/LinearAlgebra.h

@@ -93,7 +93,7 @@ class LinearAlgebra : public ILinearAlgebra
 
   /*!
    * \brief Compute L2 norm of a vector
-   * \param[in] x The vector on which norm0 is computed
+   * \param[in] x The vector on which norm2 is computed
    * \returns The norm2 of the vector
    */
   Real norm2(const IVector& x) const;

src/core/alien/core/backend/LinearAlgebraExpr.h

@@ -97,6 +97,13 @@ class LinearAlgebraExpr
    */
   Real norm2(const IVector& x) const;
 
+  /*!
+   * \brief Compute L2 (Frobenous) norm of a matrix
+   * \param[in] x The matrix on which norm2 is computed
+   * \returns The norm2 of the matrix
+   */
+  Real norm2(const IMatrix& x) const;
+
   /*!
    * \brief Compute a matrix vector product
    *
@@ -139,6 +146,24 @@ class LinearAlgebraExpr
    */
   void copy(const IVector& x, IVector& r) const;
 
+  /*!
+   * \brief Copy a matrix in another one
+   *
+   * \param[in] x The matrix to copy
+   * \param[in,out] r The copied matrix
+   */
+  void copy(const IMatrix& x, IMatrix& r) const;
+
+  /*!
+   * \brief Add two  matrices A and B
+   *
+   *
+   * \param[in] a The matrix A
+   * \param[in] b The matrix B
+   * \param[in,out] c The resulting matrix
+   */
+  void add(const IMatrix& a, IMatrix& b) const;
+
   /*!
    * \brief Compute the dot product of two vectors
    * \param[in] x The first vector
@@ -154,6 +179,13 @@ class LinearAlgebraExpr
    */
   void scal(Real alpha, IVector& x) const;
 
+  /*!
+   * \brief Scale a matrix by a factor
+   * \param[in] alpha The real value to scale with
+   * \param[in,out] x The vector to be scaled
+   */
+  void scal(Real alpha, IMatrix& a) const;
+
   /*!
    * \brief Extract the diagonal of a matrix in a vector
    * \param[in] a The matrix to extract the diagonal

src/core/alien/core/backend/LinearAlgebraExprT.h

@@ -77,6 +77,17 @@ LinearAlgebraExpr<Tag, TagV>::norm2(const IVector& x) const
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 
+template <class Tag, class TagV>
+Arccore::Real
+LinearAlgebraExpr<Tag, TagV>::norm2(const IMatrix& x) const
+{
+  const auto& mx = x.impl()->get<Tag>();
+  return m_algebra->norm2(mx);
+}
+
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+
 template <class Tag, class TagV>
 void LinearAlgebraExpr<Tag, TagV>::mult(const IMatrix& a, const IVector& x, IVector& r) const
 {
@@ -123,6 +134,24 @@ void LinearAlgebraExpr<Tag, TagV>::copy(const IVector& x, IVector& r) const
   m_algebra->copy(vx, vr);
 }
 
+template <class Tag, class TagV>
+void LinearAlgebraExpr<Tag, TagV>::copy(const IMatrix& x, IMatrix& r) const
+{
+  const auto& mx = x.impl()->get<TagV>();
+  auto& mr = r.impl()->get<TagV>(true);
+  ALIEN_ASSERT((mx.rowSpace() == mr.rowSpace()), ("Incompatible spaces"));
+  m_algebra->copy(mx, mr);
+}
+
+template <class Tag, class TagV>
+void LinearAlgebraExpr<Tag, TagV>::add(const IMatrix& a, IMatrix& b) const
+{
+  const auto& ma = a.impl()->get<Tag>();
+  auto& mb = b.impl()->get<Tag>(true);
+  ALIEN_ASSERT((ma.rowSpace() == mb.rowSpace()), ("Incompatible spaces"));
+  m_algebra->add(ma, mb);
+}
+
 /*---------------------------------------------------------------------------*/
 
 template <class Tag, class TagV>
@@ -166,6 +195,13 @@ void LinearAlgebraExpr<Tag, TagV>::scal(Real alpha, IVector& x) const
   m_algebra->scal(alpha, vx);
 }
 
+template <class Tag, class TagV>
+void LinearAlgebraExpr<Tag, TagV>::scal(Real alpha, IMatrix& A) const
+{
+  auto& ma = A.impl()->get<Tag>(true);
+  m_algebra->scal(alpha, ma);
+}
+
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 

src/core/alien/kernels/simple_csr/SimpleCSRMatrix.h

@@ -282,6 +282,23 @@ class SimpleCSRMatrix : public IMatrixImpl
     return matrix;
   }
 
+  void copy(SimpleCSRMatrix const& matrix)
+  {
+    m_is_parallel = matrix.m_is_parallel;
+    m_local_size = matrix.m_local_size;
+    m_local_offset = matrix.m_local_offset;
+    m_global_size = matrix.m_global_size;
+    m_ghost_size = matrix.m_ghost_size;
+    m_send_policy = matrix.m_send_policy;
+    m_recv_policy = matrix.m_recv_policy;
+    m_nproc = matrix.m_nproc;
+    m_myrank = matrix.m_myrank;
+    m_parallel_mng = matrix.m_parallel_mng;
+    m_trace = matrix.m_trace;
+    m_matrix.copy(matrix.m_matrix);
+    m_matrix_dist_info.copy(matrix.m_matrix_dist_info);
+  }
+
   void notifyChanges()
   {
     m_matrix.notifyChanges();

src/core/alien/kernels/simple_csr/algebra/CBLASMPIKernel.h

@@ -98,6 +98,19 @@ class CBLASMPIKernel
     return value;
   }
 
+  template <typename Distribution, typename VectorT>
+  static typename VectorT::ValueType nrm1(Distribution const& dist, const VectorT& x)
+  {
+    typedef typename VectorT::ValueType ValueType;
+    typename VectorT::ValueType value = cblas::nrm1(x.scalarizedLocalSize(),
+                                                    (ValueType*)x.getDataPtr(), 1);
+    if (dist.isParallel()) {
+      value = Arccore::MessagePassing::mpAllReduce(
+      dist.parallelMng(), Arccore::MessagePassing::ReduceSum, value);
+    }
+    return value;
+  }
+
   template <typename Distribution, typename VectorT>
   static typename VectorT::ValueType nrm2(Distribution const& dist, const VectorT& x)
   {
@@ -110,6 +123,19 @@ class CBLASMPIKernel
     }
     return std::sqrt(value);
   }
+
+  template <typename Distribution, typename MatrixT>
+  static typename MatrixT::ValueType matrix_nrm2(Distribution const& dist, const MatrixT& x)
+  {
+    typedef typename MatrixT::ValueType ValueType;
+    typename MatrixT::ValueType value = cblas::dot(x.getProfile().getNnz(),
+                                                   (ValueType*)x.data(), 1, (ValueType*)x.data(), 1);
+    if (dist.isParallel()) {
+      value = Arccore::MessagePassing::mpAllReduce(
+      dist.parallelMng(), Arccore::MessagePassing::ReduceSum, value);
+    }
+    return std::sqrt(value);
+  }
 };
 
 } // namespace Alien

src/core/alien/kernels/simple_csr/algebra/SimpleCSRInternalLinearAlgebra.cc

@@ -96,11 +96,13 @@ Real SimpleCSRInternalLinearAlgebra::norm0(const CSRVector& vx ALIEN_UNUSED_PARA
 
 /*---------------------------------------------------------------------------*/
 
-Real SimpleCSRInternalLinearAlgebra::norm1(const CSRVector& vx ALIEN_UNUSED_PARAM) const
+Real SimpleCSRInternalLinearAlgebra::norm1(const CSRVector& vx) const
 {
-  // return CBLASMPIKernel::nrm1(x.space().structInfo(),vx);
-  throw NotImplementedException(
-  A_FUNCINFO, "SimpleCSRLinearAlgebra::norm1 not implemented");
+#ifdef ALIEN_USE_PERF_TIMER
+  SentryType s(m_timer, "CSR-NORM1");
+#endif
+
+  return CBLASMPIKernel::nrm1(vx.distribution(), vx);
 }
 
 /*---------------------------------------------------------------------------*/
@@ -290,6 +292,15 @@ Real SimpleCSRInternalLinearAlgebraExpr::norm2(const CSRVector& vx) const
 
 /*---------------------------------------------------------------------------*/
 
+Real SimpleCSRInternalLinearAlgebraExpr::norm2(const CSRMatrix& mx) const
+{
+#ifdef ALIEN_USE_PERF_TIMER
+  SentryType s(m_timer, "MATRIX-CSR-NORM2");
+#endif
+  return CBLASMPIKernel::matrix_nrm2(mx.distribution(), mx);
+}
+/*---------------------------------------------------------------------------*/
+
 void SimpleCSRInternalLinearAlgebraExpr::mult(
 const CSRMatrix& ma, const CSRVector& vx, CSRVector& vr) const
 {
@@ -342,6 +353,32 @@ void SimpleCSRInternalLinearAlgebraExpr::copy(const CSRVector& vx, CSRVector& vr
   CBLASMPIKernel::copy(vx.distribution(), vx, vr);
 }
 
+void SimpleCSRInternalLinearAlgebraExpr::copy(const CSRMatrix& ma, CSRMatrix& mr) const
+{
+#ifdef ALIEN_USE_PERF_TIMER
+  SentryType s(m_timer, "MATRIX-CSR-COPY");
+#endif
+  mr.copy(ma);
+}
+
+void SimpleCSRInternalLinearAlgebraExpr::add(const CSRMatrix& ma, CSRMatrix& mr) const
+{
+#ifdef ALIEN_USE_PERF_TIMER
+  SentryType s(m_timer, "MATRIX-CSR-ADD");
+#endif
+
+  cblas::axpy(mr.getProfile().getNnz(), 1, (CSRMatrix::ValueType*)ma.data(), 1, mr.data(), 1);
+}
+
+void SimpleCSRInternalLinearAlgebraExpr::scal(Real alpha, CSRMatrix& mr) const
+{
+#ifdef ALIEN_USE_PERF_TIMER
+  SentryType s(m_timer, "MATRIX-CSR-SCAL");
+#endif
+
+  cblas::scal(mr.getProfile().getNnz(), alpha, mr.data(), 1);
+}
+
 /*---------------------------------------------------------------------------*/
 
 Real SimpleCSRInternalLinearAlgebraExpr::dot(

src/core/alien/kernels/simple_csr/algebra/SimpleCSRInternalLinearAlgebra.h

@@ -184,6 +184,7 @@ class SimpleCSRInternalLinearAlgebraExpr
   Real norm0(const Vector& x) const;
   Real norm1(const Vector& x) const;
   Real norm2(const Vector& x) const;
+  Real norm2(const Matrix& x) const;
   void mult(const Matrix& a, const Vector& x, Vector& r) const;
   void axpy(Real alpha, const Vector& x, Vector& r) const;
   void aypx(Real alpha, Vector& y, const Vector& x) const;
@@ -194,7 +195,7 @@ class SimpleCSRInternalLinearAlgebraExpr
   void reciprocal(Vector& x) const;
   void pointwiseMult(const Vector& x, const Vector& y, Vector& w) const;
 
-  // IInternalLinearAlgebra interface.
+  // IInternalLinearAlgebraExpr interface.
 
   void mult(const Matrix& a, const UniqueArray<Real>& x, UniqueArray<Real>& r) const;
   void axpy(Real alpha, UniqueArray<Real> const& x, UniqueArray<Real>& r) const;
@@ -205,6 +206,10 @@ class SimpleCSRInternalLinearAlgebraExpr
 
   void scal(Real alpha, UniqueArray<Real>& x) const;
 
+  void copy(const Matrix& a, Matrix& r) const;
+  void add(const Matrix& a, Matrix& r) const;
+  void scal(Real alpha, Matrix& r) const;
+
  private:
   // No member.
 };

src/core/alien/kernels/simple_csr/algebra/alien_cblas.h

@@ -18,6 +18,7 @@ void cblas_dcopy(const int n, const double* x, const int incx, double* y,
                  const int incy);
 void cblas_sscal(const int n, const float alpha, float* x, const int incx);
 void cblas_dscal(const int n, const double alpha, double* x, const int incx);
+double cblas_dasum(const int n, const double* x, const int incx);
 double cblas_dnrm2(const int n, const double* x, const int incx);
 
 #ifdef __cplusplus
@@ -68,6 +69,11 @@ void scal(const int n, const double alpha, double* x, const int incx)
   return cblas_dscal(n, alpha, x, incx);
 }
 
+double nrm1(const int n, const double* x, const int incx)
+{
+  return cblas_dasum(n, x, incx);
+}
+
 double nrm2(const int n, const double* x, const int incx)
 {
   return cblas_dnrm2(n, x, incx);

src/refsemantic/alien/ref/data/scalar/Matrix.h

@@ -54,6 +54,9 @@ class ALIEN_REFSEMANTIC_EXPORT Matrix final : public IMatrix
   Matrix(const Matrix&) = delete;
   Matrix& operator=(const Matrix&) = delete;
 
+  template <typename E>
+  Matrix& operator=(const E&);
+
  public:
   // Pour les expressions
   void visit(ICopyOnWriteMatrix&) const;

src/refsemantic/alien/ref/data/scalar/Vector.h

@@ -53,7 +53,6 @@ class ALIEN_REFSEMANTIC_EXPORT Vector final : public IVector
   Vector(const Vector&) = delete;
   Vector& operator=(const Vector&) = delete;
 
-  // FIXME: not implemented !
   template <typename E>
   Vector& operator=(const E&);