GURLS_EXPORT void svd(const gMat2D<float>& A, gMat2D<float>& U, gVec<float>& W, gMat2D<float>& Vt) {

    char jobu = 'S', jobvt = 'S';
    int m = A.rows();
    int n = A.cols();
    int k = std::min<int>(m, n);

    if ((int)W.getSize() < k) {
        throw gException("The length of vector W must be at least equal to the minimum dimension of the input matrix A");
    }
    if ((int)U.rows() < m || (int)U.cols() < k) {
        throw gException("Please check the dimensions of the matrix U where to store the singular vectors");
    }
    if ((int)Vt.rows() < k || (int)Vt.cols() < n) {
        throw gException("Please check the dimensions of the matrix Vt where to store the rigth singular vectors");
    }

    int lda = A.cols();
    int ldu = U.cols();
    int ldvt = Vt.cols();
    int info, lwork = std::max<int>(3*k+std::max<int>(m,n), 5*k);
    float* work = new float[lwork];
    float* copy = new float[m*n];
    A.asarray(copy, m*n);
    sgesvd_(&jobu, &jobvt, &n, &m, copy, &lda, W.getData(), Vt.getData(), &ldvt, U.getData(), &ldu, work, &lwork, &info);
    delete[] work;
    delete[] copy;
}

GURLS_EXPORT void eig(const gMat2D<float>& A, gVec<float>& Wr, gVec<float>& Wi)
{
    if (A.cols() != A.rows())
        throw gException("The input matrix A must be squared");

    float* Atmp = new float[A.getSize()];
    copy(Atmp, A.getData(), A.getSize());

    char jobvl = 'N', jobvr = 'N';
    int n = A.cols(), lda = A.cols(), ldvl = 1, ldvr = 1;
    int info, lwork = 4*n;
    float* work = new float[lwork];

    sgeev_(&jobvl, &jobvr, &n, Atmp, &lda, Wr.getData(), Wi.getData(), NULL, &ldvl, NULL, &ldvr, work, &lwork, &info);

    delete[] Atmp;
    delete[] work;

    if(info != 0)
    {
        std::stringstream str;
        str << "Eigenvalues/eigenVectors computation failed, error code " << info << ";" << std::endl;
        throw gException(str.str());
    }
}

GURLS_EXPORT void svd(const gMat2D<float>& A, gMat2D<float>& U, gVec<float>& W, gMat2D<float>& Vt)
{
    float* Ubuf;
    float* Sbuf;
    float* Vtbuf;

    int Urows, Ucols;
    int Slen;
    int Vtrows, Vtcols;

    gurls::svd(A.getData(), Ubuf, Sbuf, Vtbuf,
             A.rows(), A.cols(),
             Urows, Ucols, Slen, Vtrows, Vtcols);


    U.resize(Urows, Ucols);
    copy(U.getData(), Ubuf, U.getSize());

    W.resize(Slen);
    copy(W.getData(), Sbuf, Slen);

    Vt.resize(Vtrows, Vtcols);
    copy(Vt.getData(), Vtbuf, Vt.getSize());

    delete [] Ubuf;
    delete [] Sbuf;
    delete [] Vtbuf;
}

void BigArray<T>::setRow(unsigned long row, const gVec<T>& value)
{
    if(row >= this->numrows)
        throw gException(Exception_Index_Out_of_Bound);

    if(value.getSize() != this->numcols)
        throw gException(Exception_Inconsistent_Size);

    setMatrix(row, 0, value.getData(), 1, this->numcols);
}

GURLS_EXPORT float dot(const gVec<float>& x, const gVec<float>& y)
{
    if ( x.getSize() != y.getSize() )
        throw gException(gurls::Exception_Inconsistent_Size);

    int n = x.getSize();
    int incr = 1;

    return sdot_(&n, const_cast<float*>(x.getData()), &incr, const_cast<float*>(y.getData()), &incr);
}

GURLS_EXPORT double dot(const gVec<double>& x, const gVec<double>& y) {

    if ( x.getSize() != y.getSize() ) {
        throw gException(gurls::Exception_Inconsistent_Size);
    }

    int n = x.getSize();
    int incr = 1;

    return ddot_(&n, const_cast<double*>(x.getData()), &incr, const_cast<double*>(y.getData()), &incr);

}

GURLS_EXPORT void lu(gMat2D<float>& A, gVec<int>& pv) {

    unsigned int k = std::min<unsigned int>(A.cols(), A.rows());
    if (pv.getSize() != k) {
        throw gException("The lenghth of pv must be equal to the minimun dimension of A");
    }
    int info;
    int m = A.rows();
    int n = A.cols();
    int lda = A.cols();
    sgetrf_(&m, &n, A.getData(), &lda, pv.getData(), &info);

}

GURLS_EXPORT void eig(const gMat2D<float>& A, gMat2D<float>& V, gVec<float>& Wr, gVec<float>& Wi) {

    if (A.cols() != A.rows()) {
        throw gException("The input matrix A must be squared");
    }

    char jobvl = 'N', jobvr = 'V';
    int n = A.cols(), lda = A.cols(), ldvl = 1, ldvr = A.cols();
    int info, lwork = 4*n;
    float* work = new float[lwork];
    gMat2D<float> Atmp = A;
    gMat2D<float> Vtmp = V;
    sgeev_(&jobvl, &jobvr, &n, Atmp.getData(), &lda, Wr.getData(), Wi.getData(), NULL, &ldvl, Vtmp.getData(), &ldvr, work, &lwork, &info);
    Vtmp.transpose(V);
    delete[] work;
}

T GurlsWrapper<T>::eval(const gVec<T> &X, unsigned long *index)
{
    if(!trainedModel())
        throw gException("Error, Train Model First");

    gMat2D<T>X_mat(1, X.getSize());
    copy(X_mat.getData(), X.getData(), X.getSize());

    gMat2D<T>* pred_mat = eval(X_mat);

    const T* pred = pred_mat->getData();
    const unsigned long size = pred_mat->getSize();

    const T* max = std::max_element(pred, pred+size);
    T ret = *max;
    if(index != NULL)
        *index = max-pred;

    delete pred_mat;
    return ret;
}

void RecursiveRLSWrapper<T>::update(const gVec<T> &X, const gVec<T> &y)
{
    if(!this->trainedModel())
        throw gException("Error, Train Model First");

    RLSPrimalRecUpdate<T> optimizer;

    const unsigned long d = X.getSize();
    const unsigned long t = y.getSize();

    gMat2D<T>X_mat(1, d);
    copy(X_mat.getData(), X.getData(), d);
    gMat2D<T>y_mat(1, t);
    copy(y_mat.getData(), y.getData(), t);

    GurlsOptionsList* ret = optimizer.execute(X_mat, y_mat, *(this->opt));
    this->opt->removeOpt("optimizer");
    this->opt->addOpt("optimizer", ret);

    ++nTot;

    gMat2D<T>* xtx = new gMat2D<T>(d,d);
    gMat2D<T>* xty = new gMat2D<T>(d,t);

    dot(X.getData(), X.getData(), xtx->getData(), 1, d, 1, d, d, d, CblasTrans, CblasNoTrans, CblasColMajor);
    dot(X.getData(), y.getData(), xty->getData(), 1, d, 1, t, d, t, CblasTrans, CblasNoTrans, CblasColMajor);

    GurlsOptionsList* kernel = this->opt->template getOptAs<GurlsOptionsList>("kernel");

    const gMat2D<T>& XtX = kernel->getOptValue<OptMatrix<gMat2D<T> > >("XtX");
    const gMat2D<T>& Xty = kernel->getOptValue<OptMatrix<gMat2D<T> > >("Xty");

    axpy(d*d, (T)1.0, XtX.getData(), 1, xtx->getData(), 1);
    axpy(d*t, (T)1.0, Xty.getData(), 1, xty->getData(), 1);

    kernel->removeOpt("XtX");
    kernel->addOpt("XtX", new OptMatrix<gMat2D<T> >(*xtx));

    kernel->removeOpt("Xty");
    kernel->addOpt("Xty", new OptMatrix<gMat2D<T> >(*xty));


    unsigned long proportion = static_cast<unsigned long>(gurls::round(1.0/this->opt->getOptAsNumber("hoproportion")));

    if(nTot % proportion == 0)
    {
        const gMat2D<T>& Xva = kernel->getOptValue<OptMatrix<gMat2D<T> > >("Xva");
        const gMat2D<T>& yva = kernel->getOptValue<OptMatrix<gMat2D<T> > >("yva");

        const unsigned long nva = Xva.rows();
        const unsigned long nva_new = nva+1;


        gMat2D<T>* Xva_new = new gMat2D<T>(nva_new, d);

        const T* old_it = Xva.getData();
        T* new_it = Xva_new->getData();
        for(const T* end = new_it+(nva_new*d); new_it< end; old_it+=nva, new_it +=nva_new)
            copy(new_it, old_it, nva);

        copy(Xva_new->getData()+nva, X.getData(), d, nva_new, 1);

        kernel->removeOpt("Xva");
        kernel->addOpt("Xva", new OptMatrix<gMat2D<T> >(*Xva_new));


        gMat2D<T>* yva_new = new gMat2D<T>(nva_new, t);

        old_it = yva.getData();
        new_it = yva_new->getData();
        for(const T* end = new_it+(nva_new*t); new_it< end; old_it+=nva, new_it +=nva_new)
            copy(new_it, old_it, nva);

        copy(yva_new->getData()+nva, y.getData(), t, nva_new, 1);

        kernel->removeOpt("yva");
        kernel->addOpt("yva", new OptMatrix<gMat2D<T> >(*yva_new));

    }
}