arrayfire · shehzan10 · Aug 25, 2015 · Jul 4, 2015 · Jul 4, 2015 · Jul 4, 2015
diff --git a/CMakeModules/build_boost_compute.cmake b/CMakeModules/build_boost_compute.cmake
@@ -1,6 +1,6 @@
-SET(VER cf5e40ee265f005b144c4e0dda945734188e8542)
+SET(VER 79aa8f9086fdf6ef6db78e889de0273b0eb7bd19)
 SET(URL https://github.com/boostorg/compute/archive/${VER}.tar.gz)
-SET(MD5 a8d0bc4d70f3b7d65fea78327a4e067e)
+SET(MD5 dba3318cbdac912dddce71f2a38ffa43)
 
 SET(thirdPartyDir "${CMAKE_BINARY_DIR}/third_party")
 SET(srcDir "${thirdPartyDir}/compute-${VER}")

diff --git a/docs/details/lapack.dox b/docs/details/lapack.dox
@@ -121,6 +121,28 @@ When memory is a concern, users can perform Cholesky decomposition in place as s
 
 \snippet test/cholesky_dense.cpp ex_chol_inplace
 
+=======================================================================
+
+\defgroup lapack_factor_func_svd svd
+
+\ingroup lapack_factor_mat
+
+\brief Perform Singular Value Decomposition
+
+This function factorizes a  matrix **A** into two unitary matrices **U** and **Vt**, and a diagonal matrix **S** such that
+
+     \f$A = U * S * Vt\f$
+
+If **A** has **M** rows and **N** columns, **U** is of the size **M x M** , **V** is of size **N x N**, and **S** is of size **M x N**
+
+The arrayfire function only returns the non zero diagonal elements of **S**. To reconstruct the original matrix **A** from the individual factors, the following code snuppet can be used:
+
+
+\snippet test/svd_dense.cpp ex_svd_reg
+
+When memory is a concern, and **A** is dispensible, \ref svdInPlace can be used
+
+
 =======================================================================
 
 \defgroup lapack_solve_func_gen solve

diff --git a/examples/lin_algebra/svd.cpp b/examples/lin_algebra/svd.cpp
@@ -0,0 +1,55 @@
+/*******************************************************
+ * Copyright (c) 2014, ArrayFire
+ * All rights reserved.
+ *
+ * This file is distributed under 3-clause BSD license.
+ * The complete license agreement can be obtained at:
+ * http://arrayfire.com/licenses/BSD-3-Clause
+ ********************************************************/
+
+#include <arrayfire.h>
+#include <cstdio>
+#include <cstdlib>
+
+using namespace af;
+
+int main(int argc, char* argv[])
+{
+    try {
+        // Select a device and display arrayfire info
+        int device = argc > 1 ? atoi(argv[1]) : 0;
+        af::setDevice(device);
+        af::info();
+
+        float h_buffer[] = {1, 4, 2, 5, 3, 6 };  // host array
+        array in(2, 3, h_buffer);                        // copy host data to device
+
+        array u;
+        array s_vec;
+        array vt;
+        svd(u, s_vec, vt, in);
+
+        array s_mat = diag(s_vec, 0, false);
+        array in_recon = matmul(u, s_mat, vt(seq(2), span));
+
+        af_print(in);
+        af_print(s_vec);
+        af_print(u);
+        af_print(s_mat);
+        af_print(vt);
+        af_print(in_recon);
+
+    } catch (af::exception& e) {
+        fprintf(stderr, "%s\n", e.what());
+        throw;
+    }
+
+#ifdef WIN32  // pause in Windows
+    if (!(argc == 2 && argv[1][0] == '-')) {
+        printf("hit [enter]...");
+        fflush(stdout);
+        getchar();
+    }
+#endif
+    return 0;
+}
diff --git a/include/af/lapack.h b/include/af/lapack.h
@@ -14,6 +14,33 @@
 #ifdef __cplusplus
 namespace af
 {
+#if AF_API_VERSION >= 31
+    /**
+       C++ Interface for SVD decomposition
+
+       \param[out] u is the output array containing U
+       \param[out] s is the output array containing the diagonal values of sigma, (singular values of the input matrix))
+       \param[out] vt is the output array containing V^H
+       \param[in] in is the input matrix
+
+       \ingroup lapack_factor_func_svd
+    */
+    AFAPI void svd(array &u, array &s, array &vt, const array &in);
+#endif
+
+#if AF_API_VERSION >= 31
+    /**
+       C++ Interface for SVD decomposition
+
+       \param[out] u is the output array containing U
+       \param[out] s is the output array containing the diagonal values of sigma, (singular values of the input matrix))
+       \param[out] vt is the output array containing V^H
+       \param[inout] in is the input matrix and will contain random data after this operation
+
+       \ingroup lapack_factor_func_svd
+    */
+    AFAPI void svdInPlace(array &u, array &s, array &vt, array &in);
+#endif
 
     /**
        C++ Interface for LU decomposition in packed format
@@ -217,6 +244,34 @@ namespace af
 extern "C" {
 #endif
 
+#if AF_API_VERSION >= 31
+    /**
+       C Interface for SVD decomposition
+
+       \param[out] u is the output array containing U
+       \param[out] s is the output array containing the diagonal values of sigma, (singular values of the input matrix))
+       \param[out] vt is the output array containing V^H
+       \param[in] in is the input matrix
+
+       \ingroup lapack_factor_func_svd
+    */
+    AFAPI af_err af_svd(af_array *u, af_array *s, af_array *vt, const af_array in);
+#endif
+
+#if AF_API_VERSION >= 31
+    /**
+       C Interface for SVD decomposition
+
+       \param[out] u is the output array containing U
+       \param[out] s is the output array containing the diagonal values of sigma, (singular values of the input matrix))
+       \param[out] vt is the output array containing V^H
+       \param[inout] in is the input matrix that will contain random data after this operation
+
+       \ingroup lapack_factor_func_svd
+    */
+    AFAPI af_err af_svd_inplace(af_array *u, af_array *s, af_array *vt, af_array in);
+#endif
+
     /**
        C Interface for LU decomposition
 

diff --git a/src/api/c/svd.cpp b/src/api/c/svd.cpp
@@ -0,0 +1,128 @@
+/*******************************************************
+ * Copyright (c) 2014, ArrayFire
+ * All rights reserved.
+ *
+ * This file is distributed under 3-clause BSD license.
+ * The complete license agreement can be obtained at:
+ * http://arrayfire.com/licenses/BSD-3-Clause
+ ********************************************************/
+
+#include <af/dim4.hpp>
+#include <af/array.h>
+#include <af/lapack.h>
+
+#include <af/util.h>
+#include <af/defines.h>
+#include <err_common.hpp>
+#include <backend.hpp>
+#include <Array.hpp>
+#include <handle.hpp>
+#include <svd.hpp>
+
+using namespace detail;
+
+template <typename T>
+static inline void svd(af_array *s, af_array *u, af_array *vt, const af_array in)
+{
+    ArrayInfo info = getInfo(in);  // ArrayInfo is the base class which
+    af::dim4 dims = info.dims();
+    int M = dims[0];
+    int N = dims[1];
+
+    typedef typename af::dtype_traits<T>::base_type Tr;
+
+    //Allocate output arrays
+    Array<Tr> sA  = createEmptyArray<Tr>(af::dim4(min(M, N)));
+    Array<T > uA  = createEmptyArray<T >(af::dim4(M, M));
+    Array<T > vtA = createEmptyArray<T >(af::dim4(N, N));
+
+    svd<T, Tr>(sA, uA, vtA, getArray<T>(in));
+
+    *s = getHandle(sA);
+    *u = getHandle(uA);
+    *vt = getHandle(vtA);
+}
+
+template <typename T>
+static inline void svdInPlace(af_array *s, af_array *u, af_array *vt, af_array in)
+{
+    ArrayInfo info = getInfo(in);  // ArrayInfo is the base class which
+    af::dim4 dims = info.dims();
+    int M = dims[0];
+    int N = dims[1];
+
+    typedef typename af::dtype_traits<T>::base_type Tr;
+
+    //Allocate output arrays
+    Array<Tr> sA  = createEmptyArray<Tr>(af::dim4(min(M, N)));
+    Array<T > uA  = createEmptyArray<T >(af::dim4(M, M));
+    Array<T > vtA = createEmptyArray<T >(af::dim4(N, N));
+
+    svdInPlace<T, Tr>(sA, uA, vtA, getWritableArray<T>(in));
+
+    *s = getHandle(sA);
+    *u = getHandle(uA);
+    *vt = getHandle(vtA);
+}
+
+af_err af_svd(af_array *u, af_array *s, af_array *vt, const af_array in)
+{
+    try {
+        ArrayInfo info = getInfo(in);
+        af::dim4 dims = info.dims();
+
+        ARG_ASSERT(3, (dims.ndims() >= 0 && dims.ndims() <= 3));
+        af_dtype type = info.getType();
+
+        switch (type) {
+        case f64:
+            svd<double>(s, u, vt, in);
+            break;
+        case f32:
+            svd<float>(s, u, vt, in);
+            break;
+        case c64:
+            svd<cdouble>(s, u, vt, in);
+            break;
+        case c32:
+            svd<cfloat>(s, u, vt, in);
+            break;
+        default:
+            TYPE_ERROR(1, type);
+        }
+    }
+    CATCHALL;
+    return AF_SUCCESS;
+}
+
+af_err af_svd_inplace(af_array *u, af_array *s, af_array *vt, af_array in)
+{
+    try {
+        ArrayInfo info = getInfo(in);
+        af::dim4 dims = info.dims();
+
+        DIM_ASSERT(3, dims[0] <= dims[1]);
+        ARG_ASSERT(3, (dims.ndims() >= 0 && dims.ndims() <= 3));
+
+        af_dtype type = info.getType();
+
+        switch (type) {
+        case f64:
+            svdInPlace<double>(s, u, vt, in);
+            break;
+        case f32:
+            svdInPlace<float>(s, u, vt, in);
+            break;
+        case c64:
+            svdInPlace<cdouble>(s, u, vt, in);
+            break;
+        case c32:
+            svdInPlace<cfloat>(s, u, vt, in);
+            break;
+        default:
+            TYPE_ERROR(1, type);
+        }
+    }
+    CATCHALL;
+    return AF_SUCCESS;
+}
diff --git a/src/api/cpp/lapack.cpp b/src/api/cpp/lapack.cpp
@@ -13,6 +13,24 @@
 
 namespace af
 {
+    void svd(array &u, array &s, array &vt, const array &in)
+    {
+        af_array sl = 0, ul = 0, vtl = 0;
+        AF_THROW(af_svd(&ul, &sl, &vtl, in.get()));
+        s = array(sl);
+        u = array(ul);
+        vt = array(vtl);
+    }
+
+    void svdInPlace(array &u, array &s, array &vt, array &in)
+    {
+        af_array sl = 0, ul = 0, vtl = 0;
+        AF_THROW(af_svd_inplace(&ul, &sl, &vtl, in.get()));
+        s = array(sl);
+        u = array(ul);
+        vt = array(vtl);
+    }
+
     void lu(array &out, array &pivot, const array &in, const bool is_lapack_piv)
     {
         out = in.copy();

diff --git a/src/backend/cblas.cpp b/src/backend/cblas.cpp
@@ -23,34 +23,52 @@ static char transChar(CBLAS_TRANSPOSE Trans)
     }
 }
 
-#define GEMM_F77(X, TS, TV, TY)                                                     \
-void cblas_##X##gemm(                                                               \
-       const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA,                       \
-       const CBLAS_TRANSPOSE TransB, const int M, const int N,                      \
-       const int K, const TS alpha, const TV *A,                                    \
-       const int lda, const TV *B, const int ldb,                                   \
-       const TS beta, TV *C, const int ldc)                                         \
-{                                                                                   \
-    char aT = transChar(TransA);                                                    \
-    char bT = transChar(TransB);                                                    \
-    X##gemm_(&aT, &bT, &M, &N, &K,                                                  \
-            (const TY *)ADDR(alpha), (const TY *)A, &lda,                           \
-            (const TY *)B, &ldb,                                                    \
-            (const TY *)ADDR(beta), (TY *)C, &ldc);                                 \
-}                                                                                   \
-void cblas_##X##gemv(                                                               \
-        const CBLAS_ORDER order, const CBLAS_TRANSPOSE TransA,                      \
-        const int M, const int N,                                                   \
-        const TS alpha, const TV *A, const int lda,                                 \
-        const TV *X, const int incX, const TS beta,                                 \
-        TV *Y, const int incY)                                                      \
-{                                                                                   \
-    char aT = transChar(TransA);                                                    \
-    X##gemv_(&aT, &M, &N,                                                           \
-            (const TY *)ADDR(alpha), (const TY *)A, &lda,                           \
-            (const TY *)X, &incX,                                                   \
-            (const TY *)ADDR(beta), (TY *)Y, &incY);                                \
-}                                                                                   \
+#define GEMM_F77(X, TS, TV, TY)                                 \
+    void cblas_##X##gemm(                                       \
+        const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA,  \
+        const CBLAS_TRANSPOSE TransB, const int M, const int N, \
+        const int K, const TS alpha, const TV *A,               \
+        const int lda, const TV *B, const int ldb,              \
+        const TS beta, TV *C, const int ldc)                    \
+    {                                                           \
+        char aT = transChar(TransA);                            \
+        char bT = transChar(TransB);                            \
+        X##gemm_(&aT, &bT, &M, &N, &K,                          \
+                 (const TY *)ADDR(alpha), (const TY *)A, &lda,  \
+                 (const TY *)B, &ldb,                           \
+                 (const TY *)ADDR(beta), (TY *)C, &ldc);        \
+    }                                                           \
+    void cblas_##X##gemv(                                       \
+        const CBLAS_ORDER order, const CBLAS_TRANSPOSE TransA,  \
+        const int M, const int N,                               \
+        const TS alpha, const TV *A, const int lda,             \
+        const TV *X, const int incX, const TS beta,             \
+        TV *Y, const int incY)                                  \
+    {                                                           \
+        char aT = transChar(TransA);                            \
+        X##gemv_(&aT, &M, &N,                                   \
+                 (const TY *)ADDR(alpha), (const TY *)A, &lda,  \
+                 (const TY *)X, &incX,                          \
+                 (const TY *)ADDR(beta), (TY *)Y, &incY);       \
+    }                                                           \
+    void cblas_##X##axpy(                                       \
+        const int N, const TS alpha,                            \
+        const TV *X, const int incX,                            \
+        TV *Y, const int incY)                                  \
+    {                                                           \
+        X##axpy_(&N,                                            \
+                 (const TY *)ADDR(alpha),                       \
+                 (const TY *)X, &incX,                          \
+                 (TY *)Y, &incY);                               \
+    }                                                           \
+    void cblas_##X##scal(                                       \
+        const int N, const TS alpha,                            \
+        TV *X, const int incX)                                  \
+    {                                                           \
+        X##scal_(&N,                                            \
+                 (const TY *)ADDR(alpha),                       \
+                 (TY *)X, &incX);                               \
+    }                                                           \
 
 #define ADDR(val) &val
 GEMM_F77(s, float, float, float)