cavis/libnd4j/include/helpers/cpu/svd.cpp

/*******************************************************************************
 * Copyright (c) 2015-2018 Skymind, Inc.
 *
 * This program and the accompanying materials are made available under the
 * terms of the Apache License, Version 2.0 which is available at
 * https://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.
 *
 * SPDX-License-Identifier: Apache-2.0
 ******************************************************************************/

//
// @author Yurii Shyrma (iuriish@yahoo.com), created on 03.01.2018
//

#include <helpers/svd.h>
#include <helpers/jacobiSVD.h>
#include <helpers/biDiagonalUp.h>
#include <array/ResultSet.h>
#include <array/NDArrayFactory.h>


namespace sd {
namespace ops {
namespace helpers {


//////////////////////////////////////////////////////////////////////////
template <typename T>
SVD<T>::SVD(const NDArray& matrix, const int switchSize, const bool calcU, const bool calcV, const bool fullUV ) {

    if(matrix.rankOf() != 2 || matrix.isScalar())
        throw std::runtime_error("ops::helpers::SVD constructor: input array must be 2D matrix !");

    const int rows = matrix.sizeAt(0);
    const int cols = matrix.sizeAt(1);

    if(cols > rows) {

        _transp = true;
        _diagSize = rows;
    }
    else {

        _transp = false;
        _diagSize = cols;
    }

    _switchSize = switchSize;
    _calcU = calcU;
    _calcV = calcV;
    _fullUV = fullUV;

    if (_transp)
        math::nd4j_swap<bool>(_calcU, _calcV);

    _s = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize, 1}, matrix.getContext());
    _m = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize + 1, _diagSize}, matrix.getContext());
    _m.assign(0.);

    if (_calcU)
        _u = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize + 1, _diagSize + 1}, matrix.getContext());
    else
        _u = NDArrayFactory::create<T>(matrix.ordering(), {2, _diagSize + 1}, matrix.getContext());
    _u.assign(0.);

    if (_calcV) {
        _v = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize, _diagSize}, matrix.getContext());
        _v.assign(0.);
    }

    evalData(matrix);
}

//////////////////////////////////////////////////////////////////////////
template <typename T>
SVD<T>::SVD(const NDArray& matrix, const int switchSize, const bool calcU, const bool calcV, const bool fullUV, const char t) {

    if(matrix.rankOf() != 2 || matrix.isScalar())
        throw std::runtime_error("ops::helpers::SVD constructor: input array must be 2D matrix !");

    const int rows = matrix.sizeAt(0);
    const int cols = matrix.sizeAt(1);

    if(cols > rows) {

        _transp = true;
        _diagSize = rows;
    }
    else {

        _transp = false;
        _diagSize = cols;
    }

    _switchSize = switchSize;
    _calcU = calcU;
    _calcV = calcV;
    _fullUV = fullUV;

    if (_transp)
        math::nd4j_swap<bool>(_calcU, _calcV);

    _s = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize, 1}, matrix.getContext());
    _m = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize + 1, _diagSize}, matrix.getContext());
    _m.assign(0.f);

    if (_calcU)
        _u = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize + 1, _diagSize + 1}, matrix.getContext());
    else
        _u = NDArrayFactory::create<T>(matrix.ordering(), {2, _diagSize + 1}, matrix.getContext());
    _u.assign(0.);

    if (_calcV) {
        _v = NDArrayFactory::create<T>(matrix.ordering(), {_diagSize, _diagSize}, matrix.getContext());
        _v.assign(0.);
    }
}


//////////////////////////////////////////////////////////////////////////
template <typename T>
void SVD<T>::deflation1(int col1, int shift, int ind, int size) {

    if(ind <= 0)
        throw std::runtime_error("ops::helpers::SVD::deflation1 method: input int must satisfy condition ind > 0 !");

    int first = col1 + shift;
    T cos = _m.e<T>(first, first);
    T sin = _m.e<T>(first+ind, first);
    T denom = math::nd4j_sqrt<T, T>(cos*cos + sin*sin);

    if (denom == (T)0.) {

        _m.p(first+ind, first+ind, 0.f);
        return;
    }

    cos /= denom;
    sin /= denom;

    _m.p(first,first, denom);
    _m.p(first+ind, first, 0.f);
    _m.p(first+ind, first+ind, 0.f);

    auto rotation = NDArrayFactory::create<T>(_m.ordering(), {2, 2},  _m.getContext());
    rotation.p(0, 0, cos);
    rotation.p(0, 1, -sin);
    rotation.p(1, 0, sin);
    rotation.p(1, 1, cos);

    if (_calcU) {
        auto temp = _u({col1,col1+size+1, 0,0}, true);
        JacobiSVD<T>::mulRotationOnRight(col1, col1+ind, temp, rotation);
    }
    else
        JacobiSVD<T>::mulRotationOnRight(col1, col1+ind, _u, rotation);
}

//////////////////////////////////////////////////////////////////////////
template <typename T>
void SVD<T>::deflation2(int col1U , int col1M, int row1W, int col1W, int ind1, int ind2, int size) {

    if(ind1 >= ind2)
        throw std::runtime_error("ops::helpers::SVD::deflation2 method: input intes must satisfy condition ind1 < ind2 !");

    if(size <= 0)
        throw std::runtime_error("ops::helpers::SVD::deflation2 method: input size must satisfy condition size > 0 !");

    T cos = _m.e<T>(col1M+ind1, col1M);
    T sin = _m.e<T>(col1M+ind2, col1M);
    T denom = math::nd4j_sqrt<T,T>(cos*cos + sin*sin);

    if (denom == (T)0.)  {

      _m.p(col1M + ind1, col1M + ind1, _m.e<T>(col1M + ind2, col1M + ind2));
      return;
    }

    cos /= denom;
    sin /= denom;
    _m.p(col1M + ind1, col1M, denom);
    _m.p(col1M + ind2, col1M + ind2, _m.e<T>(col1M + ind1, col1M + ind1));
    _m.p(col1M + ind2, col1M, 0.f);

    auto rotation = NDArrayFactory::create<T>(_m.ordering(), {2, 2}, _m.getContext());
    rotation.p(0,0, cos);
    rotation.p(1,1, cos);

    rotation.p(0,1, -sin);
    rotation.p(1,0, sin);

    if (_calcU) {
        auto temp = _u({col1U,col1U+size+1, 0,0}, true);
        JacobiSVD<T>::mulRotationOnRight(col1U+ind1, col1U+ind2, temp, rotation);
    }
    else
        JacobiSVD<T>::mulRotationOnRight(col1U+ind1, col1U+ind2, _u, rotation);

    if (_calcV)  {
        auto temp = _v({row1W,row1W+size, 0,0}, true);
        JacobiSVD<T>::mulRotationOnRight(col1W+ind1, col1W+ind2, temp, rotation);
    }
}

//////////////////////////////////////////////////////////////////////////
// has effect on block from (col1+shift, col1+shift) to (col2+shift, col2+shift) inclusively
template <typename T>
void SVD<T>::deflation(int col1, int col2, int ind, int row1W, int col1W, int shift)
{

    const int len = col2 + 1 - col1;

    auto colVec0 = new NDArray(_m({col1+shift,col1+shift+len, col1+shift,col1+shift+1}, true));

    auto diagInterval = _m({col1+shift, col1+shift+len, col1+shift,col1+shift+len}, true).diagonal('c');

    const T almostZero = DataTypeUtils::min<T>();
    T maxElem;
    if(len == 1)
        maxElem = math::nd4j_abs<T>(diagInterval.template e<T>(0));
    else
        maxElem = diagInterval({1,-1, 0,0}, true).reduceNumber(reduce::AMax).template e<T>(0);
    T maxElem0 = colVec0->reduceNumber(reduce::AMax).template e<T>(0);

    T eps = math::nd4j_max<T>(almostZero, DataTypeUtils::eps<T>() * maxElem);
    T epsBig = (T)8. * DataTypeUtils::eps<T>() * math::nd4j_max<T>(maxElem0, maxElem);

    if(diagInterval.template e<T>(0) < epsBig)
        diagInterval.p(Nd4jLong(0), epsBig);

    for(int i=1; i < len; ++i)
        if(math::nd4j_abs<T>(colVec0->template e<T>(i)) < eps)
            colVec0->p(i, 0.f);

    for(int i=1; i < len; i++)
        if(diagInterval.template e<T>(i) < epsBig) {
            deflation1(col1, shift, i, len);
            for(int i = 0; i < len; ++i)
                diagInterval.p(i, _m.e<T>(col1+shift+i,col1+shift+i));
        }

    {

        bool totDefl = true;
        for(int i=1; i < len; i++)
            if(colVec0->template e<T>(i) >= almostZero) {
                totDefl = false;
                break;
            }

        int* permut = nullptr;
        ALLOCATE(permut, _m.getContext()->getWorkspace(), 3*_diagSize, int);
        {
            permut[0] = 0;
            int p = 1;

            for(int i=1; i<len; ++i)
                if(math::nd4j_abs<T>(diagInterval.template e<T>(i)) < almostZero)
                    permut[p++] = i;

            int k = 1, m = ind+1;

            for( ; p < len; ++p) {
                if(k > ind)
                    permut[p] = m++;
                else if(m >= len)
                    permut[p] = k++;
                else if(diagInterval.template e<T>(k) < diagInterval.template e<T>(m))
                    permut[p] = m++;
                else
                    permut[p] = k++;
            }
        }

        if(totDefl) {
            for(int i=1; i<len; ++i) {
                int ki = permut[i];
                if(math::nd4j_abs<T>(diagInterval.template e<T>(ki)) < almostZero || diagInterval.template e<T>(0) < diagInterval.template e<T>(ki))
                    permut[i-1] = permut[i];
                else {
                    permut[i-1] = 0;
                    break;
                }
            }
        }

        int *tInd = permut + len;
        int *tCol = permut + 2*len;

        for(int m = 0; m < len; m++) {
            tCol[m] = m;
            tInd[m] = m;
        }

        for(int i = totDefl ? 0 : 1; i < len; i++) {

            const int ki = permut[len - (totDefl ? i+1 : i)];
            const int jac = tCol[ki];

            T _e0 = diagInterval.template e<T>(jac);
            //math::nd4j_swap<T>(diagInterval)(i), (*diagInterval)(jac));
            diagInterval.p(jac, diagInterval.template e<T>(i));
            diagInterval.p(i, _e0);

            if(i!=0 && jac!=0) {
                _e0 = colVec0->template e<T>(jac);
                //math::nd4j_swap<T>((*colVec0)(i), (*colVec0)(jac));
                colVec0->p(jac, colVec0->template e<T>(i));
                colVec0->p(i, _e0);
            }

            if (_calcU) {
                auto temp1 = _u({col1,col1+len+1, col1+i,  col1+i+1}, true);
                auto temp2 = _u({col1,col1+len+1, col1+jac,col1+jac+1}, true);
                auto temp3 = temp1;
                temp1.assign(temp2);
                temp2.assign(temp3);
            }
            else {
                auto temp1 = _u({0,2, col1+i,   col1+i+1}, true);
                auto temp2 = _u({0,2, col1+jac, col1+jac+1}, true);
                auto temp3 = temp1;
                temp1.assign(temp2);
                temp2.assign(temp3);
            }

            if(_calcV) {
                auto temp1 = _v({row1W,row1W+len, col1W+i,   col1W+i+1}, true);
                auto temp2 = _v({row1W,row1W+len, col1W+jac, col1W+jac+1}, true);
                auto temp3 = temp1;
                temp1.assign(temp2);
                temp2.assign(temp3);
            }

            const int tI = tInd[i];
            tCol[tI] = jac;
            tCol[ki] = i;
            tInd[jac] = tI;
            tInd[i] = ki;
        }

        RELEASE(permut, _m.getContext()->getWorkspace());
    }

    {
        int i = len-1;

        while(i > 0 && (math::nd4j_abs<T>(diagInterval.template e<T>(i)) < almostZero || math::nd4j_abs<T>(colVec0->template e<T>(i)) < almostZero))
            --i;

        for(; i > 1; --i) {
            if( (diagInterval.template e<T>(i) - diagInterval.template e<T>(i-1)) < DataTypeUtils::eps<T>()*maxElem ) {
                if (math::nd4j_abs<T>(diagInterval.template e<T>(i) - diagInterval.template e<T>(i-1)) >= epsBig)
                    throw std::runtime_error("ops::helpers::SVD::deflation: diagonal elements are not properly sorted !");
                deflation2(col1, col1 + shift, row1W, col1W, i-1, i, len);
            }
        }
    }

    delete colVec0;
}


//////////////////////////////////////////////////////////////////////////
template <typename T>
T SVD<T>::secularEq(const T diff, const NDArray& col0, const NDArray& diag, const NDArray& permut, const NDArray& diagShifted, const T shift) {

    auto len = permut.lengthOf();
    T res = 1.;
    T item;
    for(Nd4jLong i=0; i<len; ++i) {
        auto j = permut.e<int>(i);
        item = col0.e<T>(j) / ((diagShifted.e<T>(j) - diff) * (diag.e<T>(j) + shift + diff));
        res += item * col0.e<T>(j);
    }

    return res;
}

//////////////////////////////////////////////////////////////////////////
template <typename T>
void SVD<T>::calcSingVals(const NDArray& col0, const NDArray& diag, const NDArray& permut, NDArray& singVals, NDArray& shifts, NDArray& mus) {

    auto len = col0.lengthOf();
    auto curLen = len;

    while(curLen > 1 && col0.e<T>(curLen-1) == (T)0.f)
        --curLen;

    for (Nd4jLong k = 0; k < len; ++k)  {

        if (col0.e<T>(k) == (T)0.f || curLen==1) {

            singVals.p(k, k==0 ? col0.e<T>(0) : diag.e<T>(k));
            mus.p(k, 0.f);
            shifts.p(k, k==0 ? col0.e<T>(0) : diag.e<T>(k));
            continue;
        }

        T left = diag.e<T>(k);
        T right;

        if(k==curLen-1)
            right = diag.e<T>(curLen-1) + col0.reduceNumber(reduce::Norm2).e<T>(0);
        else {

            int l = k+1;
            while(col0.e<T>(l) == (T)0.f) {
                ++l;
                if(l >= curLen)
                    throw std::runtime_error("ops::helpers::SVD::calcSingVals method: l >= curLen !");
            }

            right = diag.e<T>(l);
        }

        T mid = left + (right - left) / (T)2.;
        T fMid = secularEq(mid, col0, diag, permut, diag, 0.);
        T shift = (k == curLen-1 || fMid > (T)0.) ? left : right;

        auto diagShifted = diag - shift;

        T muPrev, muCur;
        if (shift == left) {
            muPrev = (right - left) * 0.1;
            if (k == curLen-1)
                muCur = right - left;
            else
                muCur = (right - left) * 0.5;
        }
        else {
            muPrev = -(right - left) * 0.1;
            muCur  = -(right - left) * 0.5;
        }

        T fPrev = secularEq(muPrev, col0, diag, permut, diagShifted, shift);
        T fCur = secularEq(muCur, col0, diag, permut, diagShifted, shift);

        if (math::nd4j_abs<T>(fPrev) < math::nd4j_abs<T>(fCur)) {
            math::nd4j_swap<T>(fPrev, fCur);
            math::nd4j_swap<T>(muPrev, muCur);
        }

        bool useBisection = fPrev * fCur > (T)0.;
        while (fCur != (T).0 &&
               math::nd4j_abs<T>(muCur - muPrev) > (T)8. * DataTypeUtils::eps<T>() * math::nd4j_max<T>(math::nd4j_abs<T>(muCur), math::nd4j_abs<T>(muPrev))
               && math::nd4j_abs<T>(fCur - fPrev) > DataTypeUtils::eps<T>() && !useBisection) {

            T a = (fCur - fPrev) / ((T)1./muCur - (T)1./muPrev);
            T jac = fCur - a / muCur;
            T muZero = -a/jac;
            T fZero = secularEq(muZero, col0, diag, permut, diagShifted, shift);

            muPrev = muCur;
            fPrev = fCur;
            muCur = muZero;
            fCur = fZero;

            if (shift == left  && (muCur < (T)0. || muCur > right - left))
                useBisection = true;
            if (shift == right && (muCur < -(right - left) || muCur > (T)0.))
                useBisection = true;
            if (math::nd4j_abs<T>(fCur) > math::nd4j_abs<T>(fPrev) &&  math::nd4j_abs<T>(fCur - fPrev) > (T)16. * DataTypeUtils::eps<T>())
                useBisection = true;
        }


        if (useBisection) {

            T leftShifted, rightShifted;
            if (shift == left) {
                leftShifted = DataTypeUtils::min<T>();
                rightShifted = (k==curLen-1) ? right : ((right - left) * (T)0.6);
            }
            else {

                leftShifted = -(right - left) * (T)0.6;
                rightShifted = -DataTypeUtils::min<T>();
            }

            T fLeft  = secularEq(leftShifted,  col0, diag, permut, diagShifted, shift);
            T fRight = secularEq(rightShifted, col0, diag, permut, diagShifted, shift);
            // if(fLeft * fRight >= (T)0.)
                // throw "ops::helpers::SVD::calcSingVals method: fLeft * fRight >= (T)0. !";

            while (rightShifted - leftShifted > (T)2.f * DataTypeUtils::eps<T>() * math::nd4j_max<T>(math::nd4j_abs<T>(leftShifted), math::nd4j_abs<T>(rightShifted))) {

                T midShifted = (leftShifted + rightShifted) / (T)2.;
                fMid = secularEq(midShifted, col0, diag, permut, diagShifted, shift);
                if (fLeft * fMid < (T)0.)
                    rightShifted = midShifted;
                else {
                    leftShifted = midShifted;
                    fLeft = fMid;
                }
            }
            muCur = (leftShifted + rightShifted) / (T)2.;
        }
        singVals.p(k, shift + muCur);
        shifts.p(k, shift);
        mus.p(k, muCur);
    }

}


//////////////////////////////////////////////////////////////////////////
template <typename T>
void SVD<T>::perturb(const NDArray& col0, const NDArray& diag, const NDArray& permut, const NDArray& singVals,  const NDArray& shifts, const NDArray& mus, NDArray& zhat) {

    int n = col0.lengthOf();
    int m = permut.lengthOf();
    if(m==0) {
        zhat.assign(0.);
        return;
    }

    int last = permut.e<int>(m-1);

    for (int k = 0; k < n; ++k) {

        if (col0.e<T>(k) == (T)0.f)
            zhat.p(k, (T)0.f);
        else {
            T dk   = diag.e<T>(k);
            T prod = (singVals.e<T>(last) + dk) * (mus.e<T>(last) + (shifts.e<T>(last) - dk));

            for(int l = 0; l<m; ++l) {
                int i = permut.e<int>(l);
                if(i!=k) {
                    int j = i<k ? i : permut.e<int>(l-1);
                    prod *= ((singVals.e<T>(j)+dk) / ((diag.e<T>(i)+dk))) * ((mus.e<T>(j)+(shifts.e<T>(j)-dk)) / ((diag.e<T>(i)-dk)));
                }
            }
            T tmp = math::nd4j_sqrt<T,T>(prod);
            zhat.p(k, col0.e<T>(k) > (T)0.f ? tmp : -tmp);
        }
    }
}


//////////////////////////////////////////////////////////////////////////
template <typename T>
void SVD<T>::calcSingVecs(const NDArray& zhat, const NDArray& diag, const NDArray& perm, const NDArray& singVals,
                             const NDArray& shifts, const NDArray& mus, NDArray& U, NDArray& V) {

    int n = zhat.lengthOf();
    int m = perm.lengthOf();

    for (int k = 0; k < n; ++k) {

        auto colU = new NDArray(U({0,0, k,k+1}, true));
        *colU = 0.;
        NDArray* colV = nullptr;

        if (_calcV) {
            colV = new NDArray(V({0,0, k,k+1}, true));
            *colV = 0.;
        }

        if (zhat.e<T>(k) == (T)0.f) {
            colU->p(k, 1.f);

            if (_calcV)
                colV->p(k, 1.f);
        }
        else {

            for(int l = 0; l < m; ++l) {
                int i = perm.e<int>(l);
                U.p(i,k, zhat.e<T>(i)/(((diag.e<T>(i) - shifts.e<T>(k)) - mus.e<T>(k)) )/( (diag.e<T>(i) + singVals.e<T>(k))));
            }
            U.p(n,k, 0.f);
            *colU /= colU->reduceNumber(reduce::Norm2);

            if (_calcV) {

                for(int l = 1; l < m; ++l){
                    int i = perm.e<T>(l);
                    V.p(i,k, diag.e<T>(i) * zhat.e<T>(i) / (((diag.e<T>(i) - shifts.e<T>(k)) - mus.e<T>(k)) )/( (diag.e<T>(i) + singVals.e<T>(k))));
                }
                V.p(0,k, -1.f);
                *colV /= colV->reduceNumber(reduce::Norm2);
            }
        }
        delete colU;
        if (_calcV)
            delete colV;
    }

    auto colU = U({0,0, n,n+1}, true);
    colU = 0.;
    colU.p(n, 1.);
}


//////////////////////////////////////////////////////////////////////////
template <typename T>
void SVD<T>::calcBlockSVD(int col1, int size, NDArray& U, NDArray& singVals, NDArray& V) {

    const T almostZero = DataTypeUtils::min<T>();
    auto col0 = _m({col1, col1+size, col1, col1+1}, true);
    auto diag = static_cast<const NDArray&>(_m({col1, col1+size, col1, col1+size}, true).diagonal('c'));

    diag.p(Nd4jLong(0), T(0));
    singVals = NDArrayFactory::create<T>(_m.ordering(), {size, 1}, _m.getContext());
    U = NDArrayFactory::create<T>(_u.ordering(), {size+1, size+1}, _u.getContext());
    if (_calcV)
        V = NDArrayFactory::create<T>(_v.ordering(), {size, size}, _v.getContext());

    int curSize = size;
    while(curSize > 1 && diag.template e<T>(curSize-1) == (T)0.f)
        --curSize;

    int m = 0;
    std::vector<T> indices;
    for(int k = 0; k < curSize; ++k)
        if(math::nd4j_abs<T>(col0.template e<T>(k)) > almostZero)
            indices.push_back((T)k);

    auto permut = NDArrayFactory::create<T>(_m.ordering(), {1, (int)indices.size()}, indices, _m.getContext());
    auto shifts = NDArrayFactory::create<T>(_m.ordering(), {size, 1}, _m.getContext());
    auto mus    = NDArrayFactory::create<T>(_m.ordering(), {size, 1}, _m.getContext());
    auto zhat   = NDArrayFactory::create<T>(_m.ordering(), {size, 1}, _m.getContext());

    calcSingVals(col0, diag, permut, singVals, shifts, mus);
    perturb(col0, diag, permut, singVals, shifts, mus, zhat);
    calcSingVecs(zhat, diag, permut, singVals, shifts, mus, U, V);

    for(int i=0; i<curSize-1; ++i) {

        if(singVals.e<T>(i) > singVals.e<T>(i+1)) {
            T _e0 = singVals.e<T>(i);
            T _e1 = singVals.e<T>(i+1);
            //math::nd4j_swap<T>(singVals(i),singVals(i+1));
            singVals.p(i, _e1);
            singVals.p(i+1, _e0);

            auto temp1 = U({0,0, i,i+1}, true);
            auto temp2 = U({0,0, i+1,i+2}, true);
            auto temp3 = temp1;
            temp1.assign(temp2);
            temp2.assign(temp3);

            if(_calcV) {
                auto temp1 = V({0,0, i,i+1}, true);
                auto temp2 = V({0,0, i+1,i+2}, true);
                auto temp3 = temp1;
                temp1.assign(temp2);
                temp2.assign(temp3);
            }
        }
    }

    auto temp1 = singVals({0,curSize, 0,0}, true);
    for (int e = 0; e < curSize / 2; ++e) {
        T tmp = temp1.e<T>(e);
        temp1.p(e, temp1.e<T>(curSize-1-e));
        temp1.p(curSize-1-e, tmp);
    }

    auto temp2 = U({0,0, 0,curSize}, true);
    for(int i = 0; i < curSize/2; ++i) {
        auto temp3 = temp2({0,0, i,i+1}, true);
        auto temp4 = temp2({0,0, curSize-1-i,curSize-i}, true);
        auto temp5 = temp3;
        temp3.assign(temp4);
        temp4.assign(temp5);
    }

    if (_calcV) {
        auto temp2 = V({0,0, 0,curSize}, true);
        for(int i = 0; i < curSize/2; ++i) {
            auto temp3 = temp2({0,0, i,i+1}, true);
            auto temp4 = temp2({0,0, curSize-1-i,curSize-i}, true);
            auto temp5 = temp3;
            temp3.assign(temp4);
            temp4.assign(temp5);
        }
    }
}


//////////////////////////////////////////////////////////////////////////
template<typename T>
void SVD<T>::DivideAndConquer(int col1, int col2, int row1W, int col1W, int shift) {

    // requires rows = cols + 1;
    const int n = col2 - col1 + 1;
    const int k = n/2;
    const T almostZero = DataTypeUtils::min<T>();
    T alphaK;
    T betaK;
    T r0;
    T lambda, phi, c0, s0;
    auto l = NDArrayFactory::create<T>(_u.ordering(), {1, k}, _u.getContext());
    auto f = NDArrayFactory::create<T>(_u.ordering(), {1, n-k-1}, _u.getContext());

    if(n < _switchSize) {

        JacobiSVD<T> jac(_m({col1,col1+n+1, col1,col1+n}, true), _calcU, _calcV, _fullUV);

        if (_calcU) {
            auto temp = _u({col1,col1+n+1, col1,col1+n+1}, true);
            temp.assign(jac._u);
        }
        else {
            auto temp1 = _u({0,1, col1,col1+n+1}, true);
            temp1.assign(jac._u({0,1, 0,0}, true));
            auto temp2 = _u({1,2, col1,col1+n+1}, true);
            temp2.assign(jac._u({n,n+1, 0,0}, true));
        }

        if (_calcV) {
            auto temp = _v({row1W,row1W+n, col1W,col1W+n}, true);
            temp.assign(jac._v);
        }

        auto temp = _m({col1+shift,col1+shift+n+1, col1+shift,col1+shift+n}, true);
        temp.assign(0.);
        auto diag = _m.diagonal('c');
        diag({col1+shift, col1+shift+n, 0,0}, true).assign(jac._s({0,n, 0,0}, true));

        return;
    }

    alphaK = _m.e<T>(col1 + k, col1 + k);
    betaK  = _m.e<T>(col1 + k + 1, col1 + k);

    DivideAndConquer(k + 1 + col1, col2, k + 1 + row1W, k + 1 + col1W, shift);
    DivideAndConquer(col1, k - 1 + col1, row1W, col1W + 1, shift + 1);

    if (_calcU) {
        lambda = _u.e<T>(col1 + k, col1 + k);
        phi    = _u.e<T>(col1 + k + 1, col2 + 1);
    }
    else {
        lambda = _u.e<T>(1, col1 + k);
        phi    = _u.e<T>(0, col2 + 1);
    }

    r0 = math::nd4j_sqrt<T, T>((math::nd4j_abs<T>(alphaK * lambda) * math::nd4j_abs<T>(alphaK * lambda)) + math::nd4j_abs<T>(betaK * phi) * math::nd4j_abs<T>(betaK * phi));

    if(_calcU) {
        l.assign(_u({col1+k,  col1+k+1,  col1,col1+k}, true));
        f.assign(_u({col1+k+1,col1+k+2,  col1+k+1,col1+n}, true));
    }
    else {
        l.assign(_u({1,2, col1, col1+k}, true));
        f.assign(_u({0,1, col1+k+1, col1+n}, true));
    }

    if (_calcV)
        _v.p(row1W+k, col1W, 1.f);

    if (r0 < almostZero){
        c0 = 1.;
        s0 = 0.;
    }
    else {
        c0 = alphaK * lambda / r0;
        s0 = betaK * phi / r0;
    }

    if (_calcU) {

        auto temp = _u({col1,col1+k+1, col1+k,col1+k+1}, true);
        NDArray q1(temp);

        for (int i = col1 + k - 1; i >= col1; --i) {
           auto temp = _u({col1,col1+k+1, i+1,i+2}, true);
            temp.assign(_u({col1, col1+k+1, i, i+1}, true));
        }

        _u({col1,col1+k+1, col1,col1+1}, true).assign(q1 * c0);
        _u({col1,col1+k+1, col2+1,col2+2}, true).assign(q1 * (-s0));
        _u({col1+k+1,col1+n+1, col1, col1+1}, true).assign(static_cast<const NDArray&>(_u({col1+k+1, col1+n+1, col2+1, col2+2}, true)) * s0);
        _u({col1+k+1,col1+n+1, col2+1,col2+2}, true) *= c0;
    }
    else  {

        T q1 = _u.e<T>(0, col1 + k);

        for (int i = col1 + k - 1; i >= col1; --i)
            _u.p(0, i+1, _u.e<T>(0, i));

        _u.p(0, col1, q1 * c0);
        _u.p(0, col2+1, -q1*s0);
        _u.p(1, col1, _u.e<T>(1, col2+1) * s0);
        _u.p(1, col2 + 1,  _u.e<T>(1, col2 + 1) * c0);
        _u({1,2,  col1+1, col1+k+1}, true) = 0.f;
        _u({0,1,  col1+k+1, col1+n}, true) = 0.f;
    }

   _m.p(col1 + shift, col1 + shift, r0);
    auto temp1 = _m({col1+shift+1,col1+shift+k+1, col1+shift,col1+shift+1}, true);
    temp1.assign(l*alphaK);
    auto temp2 = _m({col1+shift+k+1,col1+shift+n, col1+shift,col1+shift+1}, true);
    temp2.assign(f*betaK);

    deflation(col1, col2, k, row1W, col1W, shift);

    NDArray UofSVD, VofSVD, singVals;
    calcBlockSVD(col1 + shift, n, UofSVD, singVals, VofSVD);

    if(_calcU) {
        auto pTemp = _u({col1, col1+n+1, col1,col1+n+1}, true);
        auto temp = pTemp;
        pTemp.assign(mmul(temp, UofSVD));
    }
    else {
        auto pTemp = _u({0,0, col1,col1+n+1}, true);
        auto temp = pTemp;
        pTemp.assign(mmul(temp, UofSVD));
    }

    if (_calcV) {
        auto pTemp = _v({row1W,row1W+n, row1W,row1W+n}, true);
        auto temp = pTemp;
        pTemp.assign(mmul(temp, VofSVD));
    }

    auto blockM = _m({col1+shift,col1+shift+n, col1+shift,col1+shift+n}, true);
    blockM = 0.f;
    auto diag = blockM.diagonal('c');
    diag.assign(singVals);
}

//////////////////////////////////////////////////////////////////////////
template<typename T>
void SVD<T>::exchangeUV(const HHsequence& hhU, const HHsequence& hhV, const NDArray& U, const NDArray& V) {

    if (_calcU) {

        int colsU = _fullUV ? hhU.rows() : _diagSize;
        auto temp1 = NDArrayFactory::create<T>(_u.ordering(), {hhU.rows(), colsU}, _u.getContext());
        temp1.setIdentity();
        _u = temp1;

        auto temp2 = _u({0,_diagSize, 0,_diagSize}, true);
        temp2.assign(V({0,_diagSize, 0,_diagSize}, true));
        const_cast<HHsequence&>(hhU).mulLeft(_u);
    }

    if (_calcV) {

        int colsV = _fullUV ? hhV.rows() : _diagSize;
        auto temp1 = NDArrayFactory::create<T>(_v.ordering(), {hhV.rows(), colsV}, _v.getContext());
        temp1.setIdentity();
        _v = temp1;

        auto temp2 = _v({0,_diagSize, 0,_diagSize}, true);
        temp2.assign(U({0,_diagSize, 0,_diagSize}, true));
        const_cast<HHsequence&>(hhV).mulLeft(_v);
    }
}

//////////////////////////////////////////////////////////////////////////
template <typename T>
void SVD<T>::evalData(const NDArray& matrix) {

    const T almostZero = DataTypeUtils::min<T>();

    if(matrix.sizeAt(1) < _switchSize) {

        JacobiSVD<T> jac(matrix, _calcU, _calcV, _fullUV);

        if(_calcU)
            _u = jac._u;
        if(_calcV)
            _v = jac._v;

        _s.assign(jac._s);

        return;
    }

    T scale = matrix.reduceNumber(reduce::AMax).e<T>(0);

    if(scale == (T)0.)
        scale = 1.;

    NDArray copy;
    if(_transp)
        copy = matrix.transpose();
    else
        copy = matrix / scale;

    BiDiagonalUp biDiag(copy);

    _u = 0.;
    _v = 0.;

    auto temp1 = biDiag._HHbidiag.transpose();
    auto temp2 = _m({0,_diagSize, 0,0}, true);
    temp2.assign(temp1);


    auto temp3 = _m({_m.sizeAt(0)-1,_m.sizeAt(0), 0,0}, true);
    temp3.assign(0.);

    DivideAndConquer(0, _diagSize - 1, 0, 0, 0);

    for (int i = 0; i < _diagSize; ++i) {
        T a = math::nd4j_abs<T>(_m.e<T>(i, i));
        _s.p(i, a * scale);
        if (a < almostZero) {
            auto temp = _s({i+1,_diagSize, 0,0}, true);
            temp.assign(0.);
            break;
        }
        else if (i == _diagSize-1)
            break;
    }

    if(_transp)
        exchangeUV(biDiag.makeHHsequence('v'), biDiag.makeHHsequence('u'), _v, _u);
    else
        exchangeUV(biDiag.makeHHsequence('u'), biDiag.makeHHsequence('v'), _u, _v);
}


BUILD_SINGLE_TEMPLATE(template class ND4J_EXPORT SVD,,FLOAT_TYPES);



}
}
}