HPCC-16282 Add BLAS wrappers as standard library attributes and plugins

Signed-off-by: johnholt <john.d.holt@lexisnexis.com>

cmake_modules/FindCBLAS.cmake

@@ -0,0 +1,33 @@
+################################################################################
+#    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+################################################################################
+
+if(NOT CBLAS_FOUND)
+    if(WIN32)
+        set(cblas_lib "cblas")
+    else()
+        set(cblas_lib cblas tatlas satlas)
+    endif()
+
+    find_path(CBLAS_INCLUDE_DIR NAMES cblas.h)
+    find_library(CBLAS_LIBRARIES NAMES ${cblas_lib} PATHS /usr/lib/atlas /usr/lib64/atlas)
+
+    include(FindPackageHandleStandardArgs)
+    find_package_handle_standard_args(CBLAS
+        DEFAULT_MSG
+        CBLAS_LIBRARIES 
+        CBLAS_INCLUDE_DIR)
+    mark_as_advanced(CBLAS_INCLUDE_DIR CBLAS_LIBRARIES)
+endif()

cmake_modules/commonSetup.cmake

@@ -64,6 +64,7 @@ IF ("${COMMONSETUP_DONE}" STREQUAL "")
   option(Boost_USE_STATIC_LIBS "Use boost_regex static library for RPM BUILD" OFF)
   option(USE_OPENSSL "Configure use of OpenSSL" ON)
   option(USE_ZLIB "Configure use of zlib" ON)
+  option(USE_CBLAS "Configure use of cblas" ON)
   if (WIN32)
     option(USE_GIT "Configure use of GIT (Hooks)" OFF)
   else()
@@ -718,6 +719,15 @@ IF ("${COMMONSETUP_DONE}" STREQUAL "")
         endif()
       endif(USE_LIBXML2)
 
+      if(USE_CBLAS)
+        find_package(CBLAS)
+        if(CBLAS_FOUND)
+            add_definitions(-D_USE_CBLAS)
+        else()
+            message(FATAL_ERROR "CBLAS requested but package not found")
+        endif()
+      endif(USE_CBLAS)
+
       if(USE_ZLIB)
         find_package(ZLIB)
         if (ZLIB_FOUND)

ecllibrary/std/BLAS.ecl

@@ -0,0 +1,284 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+IMPORT lib_eclblas AS LIB_ECLBLAS;
+
+EXPORT BLAS := MODULE
+  // Types for the Block Basic Linear Algebra Sub-programs support
+  EXPORT Types := MODULE
+    EXPORT value_t      := LIB_ECLBLAS.blas_value_t;     //REAL8;
+    EXPORT dimension_t  := LIB_ECLBLAS.blas_dimension_t; //UNSIGNED4;
+    EXPORT matrix_t     := LIB_ECLBLAS.blas_matrix_t;    //SET OF REAL8
+    EXPORT Triangle     := LIB_ECLBLAS.blas_Triangle;    //ENUM(UNSIGNED1, Upper=1, Lower=2)
+    EXPORT Diagonal     := LIB_ECLBLAS.blas_Diagonal;    //ENUM(UNSIGNED1, UnitTri=1, NotUnitTri=2)
+    EXPORT Side         := LIB_ECLBLAS.blas_Side;        //ENUM(UNSIGNED1, Ax=1, xA=2)
+  END;
+
+  /**
+   * Function prototype for Apply2Cell.
+   * @param v the value
+   * @param r the row ordinal
+   * @param c the column ordinal
+   * @return the updated value
+   */
+  EXPORT Types.value_t ICellFunc(Types.value_t v,
+                                 Types.dimension_t r,
+                                 Types.dimension_t c) := v;
+
+  /**
+   * Iterate matrix and apply function to each cell
+   *@param m number of rows
+   *@param n number of columns
+   *@param x matrix
+   *@param f function to apply
+   *@return updated matrix
+   */
+  EXPORT Types.matrix_t Apply2Cells(Types.dimension_t m,
+                                    Types.dimension_t n,
+                                    Types.matrix_t x,
+                                    ICellFunc f) := FUNCTION
+    Cell := {Types.value_t v};
+    Cell applyFunc(Cell v, UNSIGNED pos) := TRANSFORM
+      r := ((pos-1)  %  m) + 1;
+      c := ((pos-1) DIV m) + 1;
+      SELF.v := f(v.v, r, c);
+    END;
+    dIn := DATASET(x, Cell);
+    dOut := PROJECT(dIn, applyFunc(LEFT, COUNTER));
+    RETURN SET(dOut, v);
+  END;
+
+  /**
+   * Absolute sum, the 1 norm of a vector.
+   *@param m the number of entries
+   *@param x the column major matrix holding the vector
+   *@param incx the increment for x, 1 in the case of an actual vector
+   *@param skipped default is zero, the number of entries stepped over
+   * to get to the first entry
+   *@return the sum of the absolute values
+   */
+  EXPORT Types.value_t
+      dasum(Types.dimension_t m, Types.matrix_t x,
+            Types.dimension_t incx, Types.dimension_t skipped=0)
+      := LIB_ECLBLAS.ECLBLAS.dasum(m, x, incx, skipped);
+
+/**
+ * alpha*X + Y
+ * @param N number of elements in vector
+ * @param alpha the scalar multiplier
+ * @param X the column major matrix holding the vector X
+ * @param incX the increment or stride for the vector
+ * @param Y the column major matrix holding the vector Y
+ * @param incY the increment or stride of Y
+ * @param x_skipped number of entries skipped to get to the first X
+ * @param y_skipped number of entries skipped to get to the first Y
+ * @return the updated matrix
+ */
+EXPORT Types.matrix_t
+       daxpy(Types.dimension_t N, Types.value_t alpha, Types.matrix_t X,
+             Types.dimension_t incX, Types.matrix_t Y, Types.dimension_t incY,
+             Types.dimension_t x_skipped=0, Types.dimension_t y_skipped=0)
+      := LIB_ECLBLAS.ECLBLAS.daxpy(N, alpha, X, incX, Y, incY, x_skipped, y_skipped);
+
+  /**
+   * alpha*op(A) op(B) + beta*C where op() is transpose
+   * @param transposeA true when transpose of A is used
+   * @param transposeB true when transpose of B is used
+   * @param M number of rows in product
+   * @param N number of columns in product
+   * @param K number of columns/rows for the multiplier/multiplicand
+   * @param alpha scalar used on A
+   * @param A matrix A
+   * @param B matrix B
+   * @param beta scalar for matrix C
+   * @param C matrix C or empty
+   */
+  EXPORT Types.matrix_t
+         dgemm(BOOLEAN transposeA, BOOLEAN transposeB,
+               Types.dimension_t M, Types.dimension_t N, Types.dimension_t K,
+               Types.value_t alpha, Types.matrix_t A, Types.matrix_t B,
+               Types.value_t beta=0.0, Types.matrix_t C=[])
+     := LIB_ECLBLAS.ECLBLAS.dgemm(transposeA, transposeB, M, N, K, alpha, A, B, beta, C);
+
+  /**
+   * Compute LU Factorization of matrix A.
+   * @param m number of rows of A
+   * @param n number of columns of A
+   * @return composite matrix of factors, lower triangle has an
+   *         implied diagonal of ones.  Upper triangle has the diagonal of the
+   *         composite.
+   */
+  EXPORT Types.matrix_t
+         dgetf2(Types.dimension_t m, Types.dimension_t n, Types.matrix_t a)
+      := LIB_ECLBLAS.ECLBLAS.dgetf2(m, n, a);
+
+  /**
+   * DPOTF2 computes the Cholesky factorization of a real symmetric
+   * positive definite matrix A.
+   *The factorization has the form
+   * A = U**T * U ,  if UPLO = 'U', or
+   * A = L  * L**T,  if UPLO = 'L',
+   * where U is an upper triangular matrix and L is lower triangular.
+   * This is the unblocked version of the algorithm, calling Level 2 BLAS.
+   * @param tri indocate whether upper or lower triangle is used
+   * @param r number of rows/columns in the square matrix
+   * @param A the square matrix
+   * @param clear clears the unused triangle
+   * @return the triangular matrix requested.
+   */
+
+  EXPORT Types.matrix_t
+         dpotf2(Types.Triangle tri, Types.dimension_t r, Types.matrix_t A,
+                         BOOLEAN clear=TRUE)
+      := LIB_ECLBLAS.ECLBLAS.dpotf2(tri, r, A, clear);
+
+  /**
+   * Scale a vector alpha
+   * @param N number of elements in the vector
+   * @param alpha the scaling factor
+   * @param X the column major matrix holding the vector
+   * @param incX the stride to get to the next element in the vector
+   * @param skipped the number of elements skipped to get to the first element
+   * @return the updated matrix
+   */
+  EXPORT Types.matrix_t
+         dscal(Types.dimension_t N, Types.value_t alpha, Types.matrix_t X,
+               Types.dimension_t incX, Types.dimension_t skipped=0)
+      := LIB_ECLBLAS.ECLBLAS.dscal(N, alpha, X, incX, skipped);
+
+  /**
+   * Implements symmetric rank update C <- alpha A**T A + beta C or
+   * c <- alpha A A**T + beta C.  C is N x N.
+   * @param tri update upper or lower triangle
+   * @param transposeA Transpose the A matrix to be NxK
+   * @param N number of rows
+   * @param K number of columns in the update matrix or transpose
+   * @param alpha the alpha scalar
+   * @param A the update matrix, either NxK or KxN
+   * @param beta the beta scalar
+   * @param C the matrix to update
+   * @param clear clear the triangle that is not updated.  BLAS assumes
+   * that symmetric matrices have only one of the triangles and this
+   * option lets you make that true.
+   */
+  EXPORT Types.matrix_t
+         dsyrk(Types.Triangle tri, BOOLEAN transposeA,
+               Types.dimension_t N, Types.dimension_t K,
+               Types.value_t alpha, Types.matrix_t A,
+               Types.value_t beta, Types.matrix_t C, BOOLEAN clear=FALSE)
+     := LIB_ECLBLAS.ECLBLAS.dsyrk(tri, transposeA, N, K, alpha, A, beta, C, clear);
+
+  /**
+   * Triangular matrix solver.  op(A) X = alpha B or X op(A) = alpha B
+   * where op is Transpose, X and B is MxN
+   * @param side side for A, Side.Ax is op(A) X = alpha B
+   * @param tri Says whether A is Upper or Lower triangle
+   * @param transposeA is op(A) the transpose of A
+   * @param diag is the diagonal an implied unit diagonal or supplied
+   * @param M number of rows
+   * @param N number of columns
+   * @param lda the leading dimension of the A matrix, either M or N
+   * @param alpha the scalar multiplier for B
+   * @param A a triangular matrix
+   * @param B the matrix of values for the solve
+   * @return the matrix of coefficients to get B.
+   */
+  EXPORT Types.matrix_t
+         dtrsm(Types.Side side, Types.Triangle tri,
+               BOOLEAN transposeA, Types.Diagonal diag,
+               Types.dimension_t M, Types.dimension_t N, Types.dimension_t lda,
+               Types.value_t alpha, Types.matrix_t A, Types.matrix_t B)
+      := LIB_ECLBLAS.ECLBLAS.dtrsm(side, tri, transposeA, diag, M, N, lda, alpha, A, B);
+
+  /**
+   * Extract the diagonal of he matrix
+   * @param m number of rows
+   * @param n number of columns
+   * @param x matrix from which to extract the diagonal
+   * @return diagonal matrix
+   */
+  EXPORT Types.matrix_t extract_diag(Types.dimension_t m,
+                                     Types.dimension_t n,
+                                     Types.matrix_t x) := FUNCTION
+    cell := {Types.value_t v};
+    cell ext(cell v, UNSIGNED pos) := TRANSFORM
+      r := ((pos-1) % m) + 1;
+      c := ((pos-1) DIV m) + 1;
+      SELF.v := IF(r=c AND r<=m AND c<=n, v.v, SKIP);
+    END;
+    diag := SET(PROJECT(DATASET(x, cell), ext(LEFT, COUNTER)), v);
+    RETURN diag;
+  END;
+
+  /**
+   * Extract the upper or lower triangle.  Diagonal can be actual or implied
+   * unit diagonal.
+   * @param m number of rows
+   * @param n number of columns
+   * @param tri Upper or Lower specifier, Triangle.Lower or Triangle.Upper
+   * @param dt Use Diagonal.NotUnitTri or Diagonal.UnitTri
+   * @param a Matrix, usually a composite from factoring
+   * @return the triangle
+   */
+  EXPORT Types.matrix_t
+         Extract_Tri(Types.dimension_t m, Types.dimension_t n,
+                     Types.Triangle tri, Types.Diagonal dt,
+                     Types.matrix_t a)
+       := LIB_ECLBLAS.ECLBLAS.Extract_Tri(m, n, tri, dt, a);
+
+  /**
+   * Generate a diagonal matrix.
+   * @param m number of diagonal entries
+   * @param v option value, defaults to 1
+   * @param X optional input of diagonal values, multiplied by v.
+   * @return a diagonal matrix
+   */
+  EXPORT Types.matrix_t
+         make_diag(Types.dimension_t m, Types.value_t v=1.0,
+                   Types.matrix_t X=[]) := LIB_ECLBLAS.ECLBLAS.make_diag(m, v, X);
+
+  // make_vec helpers
+  Cell := RECORD
+    Types.value_t v;
+  END;
+  Cell makeCell(Types.value_t v) := TRANSFORM
+    SELF.v := v;
+  END;
+  vec_dataset(Types.dimension_t m, Types.value_t v) := DATASET(m, makeCell(v));
+
+  /**
+   * Make a vector of dimension m
+   * @param m number of elements
+   * @param v the values, defaults to 1
+   * @return the vector
+   */
+  EXPORT Types.matrix_t
+      make_vector(Types.dimension_t m,
+                  Types.value_t v=1.0) := SET(vec_dataset(m, v), v);
+
+  /**
+   * The trace of the input matrix
+   * @param m number of rows
+   * @param n number of columns
+   * @param x the matrix
+   * @return the trace (sum of the diagonal entries)
+   */
+  EXPORT Types.value_t
+      trace(Types.dimension_t m,
+            Types.dimension_t n,
+            Types.matrix_t x)
+      := SUM(extract_diag(m,n,x));
+END;

ecllibrary/std/CMakeLists.txt

@@ -19,6 +19,7 @@ add_subdirectory(system)
 set(
     SRCS
     Audit.ecl
+    BLAS.ecl
     BundleBase.ecl
     Date.ecl
     File.ecl

ecllibrary/teststd/BLAS/Test_Apply2Cells.ecl

@@ -0,0 +1,23 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types AS Types;
+value_t := Types.value_t;
+dimension_t := Types.dimension_t;
+matrix_t := Types.matrix_t;
+
+matrix_t init1 := [1.0, 2.0, -3.0, 4.0, 5.0, -6.0, 7.0, 8.0, 9.0];
+value_t half(value_t v, dimension_t r, dimension_t c) := FUNCTION
+ RETURN v/2;
+END;
+value_t zero_d(value_t v, dimension_t r, dimension_t c) := FUNCTION
+  RETURN IF(r=c, 0, v);
+END;
+
+Test1_mat := BLAS.Apply2Cells(3, 3, init1, half);
+Test2_mat := BLAS.Apply2Cells(3, 3, init1, zero_d);
+
+EXPORT Test_Apply2Cells := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=22.5);
+    EXPORT Test02 := ASSERT(BLAS.dasum(9, Test2_mat, 1)=30);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_dasum.ecl

@@ -0,0 +1,11 @@
+IMPORT Std.BLAS AS BLAS;
+SET OF REAL8 init1 := [1.0, 2.0, -3.0, 4.0, 5.0, -6.0, 7.0, 8.0, 9.0];
+
+EXPORT Test_dasum := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(4, init1, 1)=10);
+    EXPORT Test02 := ASSERT(BLAS.dasum(4, init1, 2)=16);
+    EXPORT Test03 := ASSERT(BLAS.dasum(3, init1, 3, 2)=18);
+    EXPORT Test04 := ASSERT(BLAS.dasum(3, init1, 1, 3)=15);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_daxpy.ecl

@@ -0,0 +1,14 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types AS Types;
+Types.matrix_t init1 := [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0];
+Types.matrix_t init2 := [9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0];
+
+Test1_mat := BLAS.daxpy(9, 1.0, init1, 1, init2, 1);
+Test2_mat := BLAS.daxpy(3, -1.0, init1, 3, init2, 3, 2, 2);
+
+EXPORT Test_daxpy := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=90);
+    EXPORT Test02 := ASSERT(BLAS.dasum(9, Test2_mat, 1)=47);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_dgemm.ecl

@@ -0,0 +1,20 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types AS Types;
+Types.matrix_t init1 := [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0];
+Types.matrix_t init2 := [9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0];
+Types.matrix_t init3 := [1.0, 2.0, 3.0];
+Types.matrix_t init4 := [2.0, 2.0, 2.0];
+
+Test1_mat := BLAS.dgemm(FALSE, TRUE, 3, 3, 1, 1.0, init3, init4);
+Test2_mat := BLAS.dgemm(TRUE, FALSE, 1, 1, 3, 1.0, init3, init4);
+Test3_mat := BLAS.dgemm(FALSE, TRUE, 3, 1, 3, -1.0, init1, init4);
+Test4_mat := BLAS.dgemm(FALSE, TRUE, 3, 1, 3, 1.0, init1, init4, 5, init3);
+
+EXPORT Test_dgemm := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=36);
+    EXPORT Test02 := ASSERT(BLAS.dasum(1, Test2_mat, 1)=12);
+    EXPORT Test03 := ASSERT(BLAS.dasum(3, Test3_mat, 1)=90);
+    EXPORT Test04 := ASSERT(BLAS.dasum(3, Test4_mat, 1)=120);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_dgetf2.ecl

@@ -0,0 +1,13 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types AS Types;
+Types.matrix_t init_lower := [1, 2, 3, 0, 1, 4, 0, 0, 1];
+Types.matrix_t init_upper := [2, 0, 0, 3, 4, 0, 9, 16, 25];
+
+input := BLAS.dgemm(FALSE, FALSE, 3, 3, 3, 1.0, init_lower, init_upper);
+Test1_mat := BLAS.dgetf2(3, 3, input);
+
+EXPORT Test_dgetf2 := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=68);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_dpotf2.ecl

@@ -0,0 +1,16 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types AS Types;
+Types.matrix_t init1 := [4, 6, 8, 6, 13, 18, 8, 18, 29];
+Types.matrix_t init2 := [4, 0, 0, 0, 9, 0, 0, 0, 16];
+
+Test1_mat := BLAS.dpotf2(Types.Triangle.lower, 3, init1);
+Test2_mat := BLAS.dpotf2(Types.Triangle.upper, 3, init1, FALSE);
+Test3_mat := BLAS.dpotf2(Types.Triangle.upper, 3, init2);
+
+EXPORT Test_dpotf2 := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=16);
+    EXPORT Test02 := ASSERT(BLAS.dasum(9, Test2_mat, 1)=48);
+    EXPORT Test03 := ASSERT(BLAS.dasum(9, Test3_mat, 1)=9);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_dscal.ecl

@@ -0,0 +1,13 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT Std.BLAS.Types as Types;
+Types.matrix_t init1 := [1, 1, 1, 2, 2, 2, 3, 3, 3];
+Test1_mat := BLAS.dscal(9, 2.0, init1, 1); //all 9 as 1 vector
+Test2_mat := BLAS.dscal(3, 2.0, init1, 3);  // row 2 as 3x3
+Test3_mat := BLAS.dscal(3, 2.0, init1, 1, 3); // column 2
+EXPORT Test_dscal := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=36);
+    EXPORT Test02 := ASSERT(BLAS.dasum(9, Test2_mat, 1)=24);
+    EXPORT Test03 := ASSERT(BLAS.dasum(9, Test3_mat, 1)=24);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_dsyrk.ecl

@@ -0,0 +1,19 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types AS Types;
+Types.matrix_t initC := [1, 1, 1, 2, 2, 2, 3, 3, 3];
+Types.matrix_t initA := [1, 1, 1];
+Types.matrix_t mat_A := [1, 2, 2, 3, 3, 4];
+
+Test1_mat := BLAS.dsyrk(Types.Triangle.upper, FALSE, 3, 1, 1, initA, 1, initC, TRUE);
+Test2_mat := BLAS.dsyrk(Types.Triangle.lower, TRUE, 3, 1, 1, initA, 1, initC, TRUE);
+Test3_mat := BLAS.dsyrk(Types.Triangle.Upper, FALSE, 3, 2, 1, mat_a, 1, Test1_mat, FALSE);
+Test4_mat := BLAS.dsyrk(Types.Triangle.Lower, TRUE, 3, 2, 1, mat_a, 1, Test2_mat, FALSE);
+
+EXPORT Test_dsyrk := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=20);
+    EXPORT Test02 := ASSERT(BLAS.dasum(9, Test2_mat, 1)=16);
+    EXPORT Test03 := ASSERT(BLAS.dasum(9, Test3_mat, 1)=104);
+    EXPORT Test04 := ASSERT(BLAS.dasum(9, Test4_mat, 1)=96);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_dtrsm.ecl

@@ -0,0 +1,33 @@
+IMPORT Std.BLAS;
+IMPORT BLAS.Types;
+Side := Types.Side;
+Diagonal := Types.Diagonal;
+Triangle := Types.Triangle;
+Types.matrix_t left_a0 := [2, 3, 4, 0, 2, 3, 0, 0, 2];
+Types.matrix_t right_a0 := [2, 0, 0, 3, 2, 0, 4, 3, 2];
+Types.matrix_t composite_a0 := [2, 3, 4, 3, 2, 3, 4, 3, 2];
+Types.matrix_t composite_a1 := [4, 1.5, 2, 6, 4, 1.5, 8, 6, 4];
+Types.matrix_t mat_b := [4, 6, 8, 6, 13, 18, 8, 18, 29];
+
+Types.matrix_t ident := [1, 0, 0, 0, 1, 0, 0, 0, 1];
+
+Test1_mat := BLAS.dtrsm(Side.Ax, Triangle.Lower, FALSE, Diagonal.NotUnitTri,
+                        3, 3, 3, 1.0, left_a0, mat_b);
+Test2_mat := BLAS.dtrsm(Side.xA, Triangle.Upper, FALSE, Diagonal.NotUnitTri,
+                        3, 3, 3, 1.0, right_a0, mat_b);
+Test3_mat := BLAS.dtrsm(Side.Ax, Triangle.Upper, TRUE, Diagonal.NotUnitTri,
+                        3, 3, 3, 1.0, right_a0, mat_b);
+Test4_mat := BLAS.dtrsm(Side.Ax, Triangle.Lower, FALSE, Diagonal.NotUnitTri,
+                        3, 3, 3, 1.0, composite_a0, mat_b);
+Test5_mat := BLAS.dtrsm(Side.Ax, Triangle.Lower, FALSE, Diagonal.UnitTri,
+                        3, 3, 3, 1.0, composite_a1, mat_b);
+
+EXPORT Test_dtrsm := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(BLAS.dasum(9, Test1_mat, 1)=16);
+    EXPORT Test02 := ASSERT(BLAS.dasum(9, Test2_mat, 1)=16);
+    EXPORT Test03 := ASSERT(BLAS.dasum(9, Test3_mat, 1)=16);
+    EXPORT Test04 := ASSERT(BLAS.dasum(9, Test4_mat, 1)=16);
+    EXPORT Test05 := ASSERT(BLAS.dasum(9, Test5_mat, 1)=32);
+  END;
+END;

ecllibrary/teststd/BLAS/Test_helpers.ecl

@@ -0,0 +1,22 @@
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types;
+
+Types.matrix_t x := [1.0, 2.0, 3.0, 2.0, 2.0, 2.0, 4.0, 4.0, 4.0];
+Types.matrix_t init1 := [1.0, 2.0, 3.0, 4.0];
+
+Test_d0 := BLAS.extract_diag(3, 3, x);
+Test_norm1 := BLAS.dasum(3, Test_d0, 1);
+Test_trace := BLAS.trace(3, 3, x);
+Test_Diag1 := BLAS.make_diag(4, -1.0);
+Test_Diag2 := BLAS.make_diag(4, 7.0, init1);
+Test_vector := BLAS.make_vector(5, 4);
+
+EXPORT Test_helpers := MODULE
+  EXPORT TestRuntime := MODULE
+    EXPORT Test01 := ASSERT(Test_norm1=6);
+    EXPORT Test02 := ASSERT(Test_trace=7);
+    EXPORT Test03 := ASSERT(BLAS.dasum(9, Test_Diag1, 1)=4);
+    EXPORT Test04 := ASSERT(BLAS.dasum(9, Test_Diag2, 1)=70);
+    EXPORT Test05 := ASSERT(BLAS.dasum(5, Test_vector, 1)=20);
+  END;
+END;

plugins/CMakeLists.txt

@@ -16,6 +16,7 @@
 add_subdirectory (auditlib)
 add_subdirectory (debugservices)
 add_subdirectory (dmetaphone)
+add_subdirectory (eclblas)
 add_subdirectory (fileservices)
 add_subdirectory (logging)
 add_subdirectory (parselib)

plugins/eclblas/CMakeLists.txt

@@ -0,0 +1,67 @@
+################################################################################
+#    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+################################################################################
+
+
+# Component: eclblas
+
+#####################################################
+# Description:
+# ------------
+#    Cmake Input File for eclblas
+#####################################################
+
+
+project(eclblas)
+
+if(USE_CBLAS)
+    set(SRCS
+        dasum.cpp
+        daxpy.cpp
+        dgemm.cpp
+        dgetf2.cpp
+        dpotf2.cpp
+        dscal.cpp
+        dsyrk.cpp
+        dtrsm.cpp
+        eclblas.cpp
+        Extract_Tri.cpp
+        make_diag.cpp)
+
+    include_directories(
+        ./../../system/include
+        ./../../rtl/eclrtl
+        ./../../rtl/include
+        ./../../common/deftype
+        ./../../system/jlib
+        ${CBLAS_INCLUDE_DIR})
+
+    ADD_DEFINITIONS(-D_USRDLL -DECLBLAS_EXPORTS)
+
+    HPCC_ADD_LIBRARY(eclblas SHARED ${SRCS})
+    install(TARGETS eclblas DESTINATION plugins )
+    target_link_libraries(eclblas
+        eclrtl
+        ${CBLAS_LIBRARIES})
+else()
+    message(WARNING "Not building eclblas library for standard library due to lacking libcblas")
+endif(USE_CBLAS)
+
+if(PLATFORM OR CLIENTTOOLS_ONLY)
+    install(
+        FILES ${CMAKE_CURRENT_SOURCE_DIR}/lib_eclblas.ecllib
+        DESTINATION plugins
+        COMPONENT Runtime)
+endif()

plugins/eclblas/Extract_Tri.cpp

@@ -0,0 +1,46 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+//Extract the upper triangular matrix or the lower triangular matrix
+//from a composite.  Composites are produced by some factorizations.
+//
+#include "eclblas.hpp"
+
+ECLBLAS_CALL void Extract_Tri(bool & __isAllResult, size32_t & __lenResult,
+                              void * & __result, uint32_t m, uint32_t n, uint8_t tri,
+                              uint8_t dt, bool isAllA, size32_t lenA, const void * a){
+  int cells = m * n;
+  __isAllResult = false;
+  __lenResult = lenA;
+  double *new_a = (double*) rtlMalloc(lenA);
+  unsigned int r=0;    //row
+  unsigned int c=0;    //column
+  for (int i=0; i<cells; i++) {
+    double a_element = ((const double*)a)[i];
+    if (r==c) new_a[i] = (dt==UNIT_TRI) ? 1.0  : a_element;
+    else if (r > c) { // lower part
+      new_a[i] = (tri==UPPER) ? 0.0  : a_element;
+    } else {          // upper part
+      new_a[i] = (tri==UPPER) ? a_element  : 0.0;
+    }
+    r++;
+    if (r==m) {
+      r=0;
+      c++;
+    }
+  }
+  __result = (void*) new_a;
+}

plugins/eclblas/README.md

@@ -0,0 +1,53 @@
+ECL BLAS Plugin
+================
+
+This is the ECL plugin to utilize the Basic Linear Algebra Sysem (BLAS).
+It utilizes the C wrappers for BLAS
+
+Installation and Dependencies
+----------------------------
+
+To build the <plugin name> plugin with the HPCC-Platform, <dependency> is required.
+```
+sudo apt-get install atlas-dev
+```
+
+Getting started
+---------------
+
+The libcblas.so library can be found in several packages including ATLAS
+and LAPACK.
+
+The Actual Plugin
+-----------------
+
+The eclblas plugin features are exposed via the ECL Standard library
+Std.BLAS module.
+
+
+###An ECL Example
+```c
+IMPORT Std.BLAS AS BLAS;
+IMPORT BLAS.Types AS Types;
+Types.matrix_t init1 := [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0];
+Types.matrix_t init2 := [9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0];
+Types.matrix_t init3 := [1.0, 2.0, 3.0];
+Types.matrix_t init4 := [2.0, 2.0, 2.0];
+
+Test1_mat := BLAS.dgemm(FALSE, TRUE, 3, 3, 1, 1.0, init3, init4);
+Test2_mat := BLAS.dgemm(TRUE, FALSE, 1, 1, 3, 1.0, init3, init4);
+Test3_mat := BLAS.dgemm(FALSE, TRUE, 3, 1, 3, -1.0, init1, init4);
+Test4_mat := BLAS.dgemm(FALSE, TRUE, 3, 1, 3, 1.0, init1, init4, 5, init3);
+OUTPUT(Test1_mat);
+OUTPUT(Test2_mat);
+OUTPUT(Test3_mat);
+OUTPUT(test4_mat);
+```
+Yields the various matrix products.
+
+Behavior and Implementation Details
+------------------------------------
+
+All of the arrays are column major.  The underlying Fortran code
+expects column major, so we do not need the library to interpret
+the data as row major (usual C convention).

plugins/eclblas/dasum.cpp

@@ -0,0 +1,28 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+
+
+// Absolute sum.  the L1-norm of a vector
+
+#include "eclblas.hpp"
+
+ECLBLAS_CALL double dasum(uint32_t m, bool isAllX, size32_t lenX, const void * x,
+                          uint32_t incx, uint32_t skipped) {
+  const double* X = ((const double*)x) + skipped;
+  double rslt = cblas_dasum(m, X, incx);
+  return rslt;
+}

plugins/eclblas/daxpy.cpp

@@ -0,0 +1,35 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+
+// Vector add, alpha X  +  Y
+#include "eclblas.hpp"
+
+
+ECLBLAS_CALL void daxpy(bool & __isAllResult, size32_t & __lenResult,
+                        void * & __result, uint32_t n, double alpha,
+                        bool isAllX, size32_t lenX, const void * x, uint32_t incx,
+                        bool isAllY, size32_t lenY, const void* y, uint32_t incy,
+                        uint32_t x_skipped, uint32_t y_skipped) {
+  __isAllResult = false;
+  __lenResult = lenY;
+  const double* X = ((const double*)x) + x_skipped;
+  double *result = (double*) rtlMalloc(__lenResult);
+  memcpy(result, y,lenY);
+  double* Y = result + y_skipped;
+  cblas_daxpy(n, alpha, X, incx, Y, incy);
+  __result = (void*) result;
+}

plugins/eclblas/dgemm.cpp

@@ -0,0 +1,45 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+
+
+// matrix-matrix multiply.  alpha A B + beta C with flags to
+//transpose A and B.  beta defaults to zero, and C to empty
+#include "eclblas.hpp"
+
+ECLBLAS_CALL void dgemm(bool & __isAllResult, size32_t & __lenResult,
+                        void * & __result, bool transposeA, bool transposeB,
+                        uint32_t m, uint32_t n, uint32_t k,
+                        double alpha, bool isAllA, size32_t lenA, const void* A,
+                        bool isAllB, size32_t lenB, const void* B, double beta,
+                        bool isAllC, size32_t lenC, const void* C) {
+  unsigned int lda = transposeA==0 ? m  : k;
+  unsigned int ldb = transposeB==0 ? k  : n;
+  unsigned int ldc = m;
+  __isAllResult = false;
+  __lenResult = m * n * sizeof(double);
+  double *result = (double*) rtlMalloc(__lenResult);
+  // populate if provided
+  for(uint32_t i=0; i<m*n; i++) result[i] = (__lenResult==lenC) ?((double*)C)[i] :0.0;
+  cblas_dgemm(CblasColMajor,
+              transposeA ? CblasTrans : CblasNoTrans,
+              transposeB ? CblasTrans : CblasNoTrans,
+              m, n, k, alpha,
+              (const double *) A, lda,
+              (const double *) B, ldb,
+              beta, result, ldc);
+  __result = (void *) result;
+  }

plugins/eclblas/dgetf2.cpp

@@ -0,0 +1,61 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+//DGETF2 computes the LU factorization of a matrix A.  Similar to LAPACK routine
+//of same name.  Result matrix holds both Upper and Lower triangular matrix, with
+//lower matrix diagonal implied since it is a unit triangular matrix.
+//This version does not permute the rows.
+//This version does not support a sub-matrix, hence no LDA argument.
+//This routine would be better if dlamch were available to determine safe min
+
+#include <math.h>
+#include "eclblas.hpp"
+
+ECLBLAS_CALL void dgetf2(bool & __isAllResult, size32_t & __lenResult,
+                         void * & __result, uint32_t m, uint32_t n,
+                         bool isAllA, size32_t lenA, const void* a) {
+  //double sfmin = dlamch('S');   // get safe minimum
+  unsigned int cells = m*n;
+  __isAllResult = false;
+  __lenResult = cells * sizeof(double);
+  double *new_a = (double*) rtlMalloc(__lenResult);
+  memcpy(new_a, a, __lenResult);
+  double akk;
+  unsigned int i, k;
+  unsigned int diag, vpos, wpos, mpos;
+  unsigned int sq_dim = (m < n) ? m  : n;
+  for (k=0; k<sq_dim; k++) {
+    diag = (k*m) + k;     // diag cell
+    vpos = diag + 1;      // top cell of v vector
+    wpos = diag + m;      // left cell of w vector
+    mpos = diag + m + 1;  //upper left of sub-matrix to update
+    akk = new_a[diag];
+    if (akk == 0.0) {
+      rtlFree(new_a);
+      rtlFail(0, "Permute required"); // need to permute
+    }
+    //Ideally, akk should be tested against sfmin, and dscal used
+    // to update the vector for the L cells.
+    for (i=vpos; i<vpos+m-k-1; i++) new_a[i] = new_a[i]/akk;
+    //Update sub-matrix
+    if (k < sq_dim - 1) {
+      cblas_dger(CblasColMajor,
+                 m-k-1, n-k-1, -1.0,  // sub-matrix dimensions
+                 (new_a+vpos), 1, (new_a+wpos), m, (new_a+mpos), m);
+    }
+  }
+  __result = (void*) new_a;
+}

plugins/eclblas/dpotf2.cpp

@@ -0,0 +1,76 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+//DPOTF2 computes the Cholesky factorization of a real symmetric
+// positive definite matrix A.
+//The factorization has the form
+// A = U**T * U ,  if UPLO = 'U', or
+// A = L  * L**T,  if UPLO = 'L',
+// where U is an upper triangular matrix and L is lower triangular.
+// This is the unblocked version of the algorithm, calling Level 2 BLAS.
+//
+
+#include "eclblas.hpp"
+#include <math.h>
+
+ECLBLAS_CALL void dpotf2(bool & __isAllResult, size32_t & __lenResult,
+                         void * & __result, uint8_t tri, uint32_t r,
+                         bool isAllA, size32_t lenA, const void * A,
+                         bool clear) {
+  unsigned int cells = r*r;
+  __isAllResult = false;
+  __lenResult = cells * sizeof(double);
+  double *new_a = (double*) rtlMalloc(__lenResult);
+  memcpy(new_a, A, __lenResult);
+  double ajj;
+  // x and y refer to the embedded vectors for the multiply, not an axis
+  unsigned int diag, a_pos, x_pos, y_pos;
+  unsigned int col_step = r;  // between columns
+  unsigned int row_step = 1;  // between rows
+  unsigned int x_step = (tri==UPPER_TRIANGLE)  ? row_step  : col_step;
+  unsigned int y_step = (tri==UPPER_TRIANGLE)  ? col_step  : row_step;
+  for (unsigned int j=0; j<r; j++) {
+    diag = (j * r) + j;    // diagonal
+    x_pos = j * ((tri==UPPER_TRIANGLE) ? col_step  : row_step);
+    a_pos = (j+1) * ((tri==UPPER_TRIANGLE) ? col_step  : row_step);
+    y_pos = diag + y_step;
+    // ddot.value <- x'*y
+    ajj = new_a[diag] - cblas_ddot(j, (new_a+x_pos), x_step, (new_a+x_pos), x_step);
+    //if ajj is 0, negative or NaN, then error
+    if (ajj <= 0.0) {
+      rtlFree(new_a);
+      rtlFail(0, "Not a positive definite matrix");
+    }
+    ajj = sqrt(ajj);
+    new_a[diag] = ajj;
+    if ( j < r-1) {
+      // y <- alpha*op(A)*x + beta*y
+      cblas_dgemv(CblasColMajor,
+                  (tri==UPPER_TRIANGLE)  ? CblasTrans  : CblasNoTrans,
+                  (tri==UPPER_TRIANGLE)  ? j           : r-1-j,    // M
+                  (tri==UPPER_TRIANGLE)  ? r-1-j       : j,        // N
+                   -1.0,                          // alpha
+                   (new_a+a_pos), r,              //A
+                   (new_a+x_pos), x_step,         //X
+                   1.0, (new_a+y_pos), y_step);   // beta and Y
+      // x <- alpha * x
+      cblas_dscal(r-1-j, 1.0/ajj, (new_a+y_pos), y_step);
+    }
+    // clear lower or upper part if clear flag set
+    for(unsigned int k=1; clear && k<r-j; k++) new_a[(k*x_step)+diag] = 0.0;
+  }
+  __result = (void*) new_a;
+}

plugins/eclblas/dscal.cpp

@@ -0,0 +1,31 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+// scale a vector
+
+#include "eclblas.hpp"
+
+ECLBLAS_CALL void dscal(bool & __isAllResult, size32_t & __lenResult,
+                        void * & __result, uint32_t n, double alpha,
+                        bool isAllX, size32_t lenX, const void * x,
+                        uint32_t incx, uint32_t skipped) {
+  double *result = (double*) rtlMalloc(lenX);
+  memcpy(result, x, lenX);
+  cblas_dscal(n, alpha, result+skipped, incx);
+  __result = (void*) result;
+  __isAllResult = false;
+  __lenResult = lenX;
+}

plugins/eclblas/dsyrk.cpp

@@ -0,0 +1,47 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+// Implements symmetric rank update.  C <- alpha A**T * A  + beta C
+//or C <- alpha A * A**T  + beta C.  Triangular parameters says whether
+//the update is upper or lower.  C is N by N
+#include "eclblas.hpp"
+
+ECLBLAS_CALL void dsyrk(bool & __isAllResult, size32_t & __lenResult,
+                        void * &__result, uint8_t tri, bool transposeA,
+                        uint32_t n, uint32_t k, double alpha, bool isAllA,
+                        size32_t lenA, const void * a, double beta,
+                        bool isAllC, size32_t lenC, const void * c,
+                        bool clear) {
+  __isAllResult = false;
+  __lenResult = lenC;
+  double *new_c = (double*) rtlMalloc(lenC);
+  if (clear) {
+    unsigned int pos = 0;
+    for(unsigned int i=0; i<n; i++) {
+      pos = i*n;  // pos is head of column
+      for (unsigned int j=0; j<n; j++) {
+        new_c[pos+j] = tri==UPPER ? i>=j ? ((double*)c)[pos+j]  : 0.0
+                                  : i<=j ? ((double*)c)[pos+j]  : 0.0;
+      }
+    }
+  } else memcpy(new_c, c, __lenResult);
+  unsigned int lda = (transposeA)  ? k  : n;
+  cblas_dsyrk(CblasColMajor,
+              tri==UPPER  ? CblasUpper  : CblasLower,
+              transposeA ? CblasTrans : CblasNoTrans,
+              n, k, alpha, (const double *)a, lda, beta, new_c, n);
+  __result = (void*) new_c;
+}

plugins/eclblas/dtrsm.cpp

@@ -0,0 +1,40 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+// Triangular matrix solver.
+//  op( A )*X = alpha*B,   or   X*op( A ) = alpha*B; B is m x n
+//
+#include "eclblas.hpp"
+
+ECLBLAS_CALL void dtrsm(bool & __isAllResult, size32_t & __lenResult,
+                        void * & __result, uint8_t side, uint8_t tri,
+                        bool transposeA, uint8_t diag, uint32_t m,
+                        uint32_t n, uint32_t lda, double alpha, bool isAllA,
+                        size32_t lenA, const void * a, bool isAllB, size32_t lenB,
+                        const void * b) {
+  unsigned int ldb = m;
+  __isAllResult = false;
+  __lenResult = lenB;
+  double *new_b = (double*) rtlMalloc(lenB);
+  memcpy(new_b, b, __lenResult);
+  cblas_dtrsm(CblasColMajor,
+              side==AX ?  CblasLeft  : CblasRight,
+              tri==UPPER  ? CblasUpper  : CblasLower,
+              transposeA ? CblasTrans : CblasNoTrans,
+              diag==UNIT ? CblasUnit : CblasNonUnit,
+              m, n, alpha, (const double *)a, lda, new_b, ldb);
+  __result = (void*) new_b;
+}

plugins/eclblas/eclblas.cpp

@@ -0,0 +1,36 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+
+#include "platform.h"
+#include "eclrtl.hpp"
+#include "jstring.hpp"
+#include "eclblas.hpp"
+
+#define ECLBLAS_PLUGIN_VERSION "eclblas plugin 1.0.0"
+
+ECLBLAS_PLUGIN_API bool getECLPluginDefinition(ECLPluginDefinitionBlock *pb)
+{
+    if (pb->size != sizeof(ECLPluginDefinitionBlock))
+        return false;
+    pb->magicVersion = PLUGIN_VERSION;
+    pb->version = ECLBLAS_PLUGIN_VERSION;
+    pb->moduleName = "lib_eclblas";
+    pb->ECL = NULL;
+    pb->flags = PLUGIN_IMPLICIT_MODULE;
+    pb->description = "ECL plugin library for BLAS\n";
+    return true;
+}

plugins/eclblas/eclblas.hpp

@@ -0,0 +1,105 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+
+#ifndef ECLBLAS_INCL
+#define ECLBLAS_INCL
+
+#ifdef _WIN32
+#define ECLBLAS_CALL _cdecl
+#ifdef ECLBLAS_EXPORTS
+#define ECLBLAS_PLUGIN_API __declspec(dllexport)
+#else
+#define ECLBLAS_PLUGIN_API __declspec(dllimport)
+#endif
+#else
+#define ECLBLAS_CALL
+#define ECLBLAS_PLUGIN_API
+#endif
+
+#include "platform.h"
+#include "hqlplugins.hpp"
+#include "eclinclude4.hpp"
+#include "eclrtl.hpp"
+#include "eclhelper.hpp"
+
+
+// Defines for triangle, diagonal, Side
+#define UPPER_TRIANGLE 1
+#define UPPER 1
+#define AX 1
+#define UNIT 1
+#define UNIT_TRI 1
+
+
+extern "C" {
+#include <cblas.h>
+}
+
+extern "C" {
+  ECLBLAS_PLUGIN_API bool getECLPluginDefinition(ECLPluginDefinitionBlock *pb);
+
+  ECLBLAS_CALL double dasum(uint32_t m, bool isAllX, size32_t lenX, const void * x,
+                            uint32_t incx, uint32_t skipped);
+
+  ECLBLAS_CALL void daxpy(bool & __isAllResult, size32_t & __lenResult,
+                          void * & __result, uint32_t n, double alpha,
+                          bool isAllX, size32_t lenX, const void * x, uint32_t incx,
+                          bool isAllY, size32_t lenY, const void * y, uint32_t incy,
+                          uint32_t x_skipped, uint32_t y_skipped);
+
+  ECLBLAS_CALL void dgemm(bool & __isAll_Result, size32_t & __lenResult,
+                          void * & __result, bool transposeA, bool transposeB,
+                          uint32_t m, uint32_t n, uint32_t k,
+                          double alpha, bool isAllA, size32_t lenA, const void* A,
+                          bool isAllB, size32_t lenB, const void* B, double beta,
+                          bool isAllC, size32_t lenC, const void* C);
+
+  ECLBLAS_CALL void dgetf2(bool & __isAllResult, size32_t & __lenResult,
+                           void * & result, uint32_t m, uint32_t n,
+                           bool isAllA, size32_t lenA, const void* a);
+
+  ECLBLAS_CALL void dpotf2(bool & __isAllResult, size32_t & __lenResult,
+                           void * & __result, uint8_t tri, uint32_t r,
+                           bool isAllA, size32_t lenA, const void * A,
+                           bool clear);
+  ECLBLAS_CALL void dscal(bool & __isAllResult, size32_t & __lenResult,
+                          void * & __result, uint32_t n, double alpha,
+                          bool isAllX, size32_t lenX, const void * X,
+                          uint32_t incx, uint32_t skipped);
+  ECLBLAS_CALL void dsyrk(bool & __isAllResult, size32_t & __lenResult,
+                          void * &__result, uint8_t tri, bool transposeA,
+                          uint32_t N, uint32_t k, double alpha, bool isAllA,
+                          size32_t lenA, const void * a, double beta,
+                          bool isAllC, size32_t lenC, const void * c,
+                          bool clear);
+  ECLBLAS_CALL void dtrsm(bool & __isAllResult, size32_t & __lenResult,
+                        void * & __result, uint8_t side, uint8_t tri,
+                        bool transposeA, uint8_t diag, uint32_t m,
+                        uint32_t n, uint32_t lda, double alpha, bool isAllA,
+                        size32_t lenA, const void * a, bool isAllB, size32_t lenB,
+                        const void * b);
+  ECLBLAS_CALL void Extract_Tri(bool & __isAllResult, size32_t & __lenResult,
+                                void * & __result, uint32_t m, uint32_t n,
+                                uint8_t dt, bool isAllA, size32_t lenA,
+                                const void * a);
+  ECLBLAS_CALL void make_diag(bool & __isAllResult, size32_t & __lenResult,
+                              void * & __result, size32_t m, double v,
+                              bool isAllX, size32_t lenX, const void * x);
+}
+
+
+#endif

plugins/eclblas/lib_eclblas.ecllib

@@ -0,0 +1,54 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the License);
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an AS IS BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+
+EXPORT blas_dimension_t  := UNSIGNED4;
+EXPORT blas_value_t      := REAL8;
+EXPORT blas_matrix_t     := SET OF REAL8;
+EXPORT blas_Triangle     := ENUM(UNSIGNED1, Upper=1, Lower=2);
+EXPORT blas_Diagonal     := ENUM(UNSIGNED1, UnitTri=1, NotUnitTri=2);
+EXPORT blas_Side         := ENUM(UNSIGNED1, Ax=1, xA=2);
+
+
+EXPORT eclblas := SERVICE : plugin('eclblasplugin'), library('eclblas')
+  REAL8 dasum(blas_dimension_t m, const blas_matrix_t x, blas_dimension_t incx,
+                blas_dimension_t skipped=0) : C, PURE;
+  SET OF REAL8 daxpy(blas_dimension_t N, blas_value_t alpha, const blas_matrix_t X,
+                 blas_dimension_t incX, blas_matrix_t Y, blas_dimension_t incY,
+                 blas_dimension_t x_skipped=0, blas_dimension_t y_skipped=0)
+                 : C, PURE;
+  SET OF REAL8 dgemm(BOOLEAN transposeA, BOOLEAN transposeB,
+                 blas_dimension_t M, blas_dimension_t N, blas_dimension_t K,
+                 blas_value_t alpha, const blas_matrix_t A, const blas_matrix_t B,
+                 blas_value_t beta=0.0, const blas_matrix_t C=[]) : C, PURE;
+  SET OF REAL8 dgetf2(blas_dimension_t m, blas_dimension_t n, const blas_matrix_t a)
+                 : C, PURE;
+  SET OF REAL8 dpotf2(blas_Triangle tri, blas_dimension_t r, blas_matrix_t A,
+                  BOOLEAN clear=TRUE) : C, PURE;
+  SET OF REAL8 dscal(blas_dimension_t N, blas_value_t alpha, const blas_matrix_t X,
+                 blas_dimension_t incX, blas_dimension_t skipped=0) : C, PURE;
+  SET OF REAL8 dsyrk(blas_Triangle tri, BOOLEAN transposeA, blas_dimension_t N,
+                 blas_dimension_t K, blas_value_t alpha, const blas_matrix_t A,
+                 blas_value_t beta, const blas_matrix_t C, BOOLEAN clear=FALSE)
+                 : C, PURE;
+  SET OF REAL8 dtrsm(blas_Side side, blas_Triangle tri,
+                 BOOLEAN transposeA, blas_Diagonal diag,
+                 blas_dimension_t M, blas_dimension_t N,  blas_dimension_t lda,
+                 blas_value_t alpha, const blas_matrix_t A, const blas_matrix_t B) : C, PURE;
+  SET OF REAL8 Extract_Tri(blas_dimension_t m, blas_dimension_t n, blas_Triangle tri,
+                       blas_Diagonal dt, const blas_matrix_t a) : C, PURE;
+  SET OF REAL8 make_diag(blas_dimension_t m, blas_value_t v=1.0, const blas_matrix_t X=[]) : C, PURE;
+END;

plugins/eclblas/make_diag.cpp

@@ -0,0 +1,41 @@
+/*##############################################################################
+
+    HPCC SYSTEMS software Copyright (C) 2016 HPCC Systems®.
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+############################################################################## */
+//Make a diagonal matrix from a vector or a single value
+//If the vector is present the diagonal is the product
+#include "eclblas.hpp"
+
+ECLBLAS_CALL void make_diag(bool & __isAllResult, size32_t & __lenResult,
+                            void * & __result, uint32_t m, double v,
+                            bool isAllX, size32_t lenX, const void * x) {
+  int cells = m * m;
+  __isAllResult = false;
+  __lenResult = cells * sizeof(double);
+  if (lenX > 0 && m*sizeof(double) != lenX) rtlFail(0,"Bad X vector length");
+  double *diag = (double*) rtlMalloc(__lenResult);
+  const double *in_x = (const double*)x;
+  unsigned int r = 0;   // row
+  unsigned int c = 0;   //row and column
+  for (int i=0; i<cells; i++) {
+    diag[i] = (r==c)?  (lenX!=0 ? v*in_x[r] : v) : 0.0;
+    r++;
+    if (r==m) {
+      r = 0;
+      c++;
+    }
+  }
+  __result = (void*) diag;
+}

plugins/eclblas/sourcedoc.xml

@@ -0,0 +1,26 @@
+<?xml version="1.0" encoding="utf-8"?>
+<!--
+################################################################################
+#    HPCC SYSTEMS software Copyright (C) 2012 HPCC Systems®.
+#
+#    Licensed under the Apache License, Version 2.0 (the "License");
+#    you may not use this file except in compliance with the License.
+#    You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#    Unless required by applicable law or agreed to in writing, software
+#    distributed under the License is distributed on an "AS IS" BASIS,
+#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#    See the License for the specific language governing permissions and
+#    limitations under the License.
+################################################################################
+-->
+<!DOCTYPE section PUBLIC "-//OASIS//DTD DocBook XML V4.3//EN" "http://www.oasis-open.org/docbook/xml/4.3/docbookx.dtd">
+<section>
+    <title>plugins/eclblas</title>
+
+    <para>
+        The plugins/eclblas directory contains the sources for the plugins/blas library.
+    </para>
+</section>

plugins/sourcedoc.xml

@@ -25,6 +25,7 @@
     </para>
     <xi:include xmlns:xi="http://www.w3.org/2001/XInclude" href="auditlib/sourcedoc.xml"/>
     <xi:include xmlns:xi="http://www.w3.org/2001/XInclude" href="debugservices/sourcedoc.xml"/>
+    <xi:include xmlns:xi="http://www.w3.org/2001/XInclude" href="eclblas/sourcedoc.xml"/>
     <xi:include xmlns:xi="http://www.w3.org/2001/XInclude" href="examplelib/sourcedoc.xml"/>
     <xi:include xmlns:xi="http://www.w3.org/2001/XInclude" href="fileservices/sourcedoc.xml"/>
     <xi:include xmlns:xi="http://www.w3.org/2001/XInclude" href="logging/sourcedoc.xml"/>