fix bad performance with sparse() function 94/15894/5
Antoine ELIAS [Mon, 2 Feb 2015 09:10:48 +0000 (10:10 +0100)]
and fix umfpack tests

Change-Id: I605475cffeb9defbfdaed6719914a79f6794d032

scilab/modules/api_scilab/src/cpp/api_sparse.cpp
scilab/modules/ast/includes/types/sparse.hxx
scilab/modules/ast/src/cpp/operations/types_divide.cpp
scilab/modules/ast/src/cpp/types/sparse.cpp
scilab/modules/sparse/sci_gateway/cpp/sci_sparse.cpp

index 6e6085d..6c1c9eb 100644 (file)
@@ -1,17 +1,17 @@
 /*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2009 - DIGITEO - Antoine ELIAS
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution.  The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
- *
- * Please note that piece of code will be rewrited for the Scilab 6 family
- * However, the API (profile of the functions in the header files) will be
- * still available and supported in Scilab 6.
- */
+* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+* Copyright (C) 2009 - DIGITEO - Antoine ELIAS
+*
+* This file must be used under the terms of the CeCILL.
+* This source file is licensed as described in the file COPYING, which
+* you should have received as part of this distribution.  The terms
+* are also available at
+* http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+*
+* Please note that piece of code will be rewrited for the Scilab 6 family
+* However, the API (profile of the functions in the header files) will be
+* still available and supported in Scilab 6.
+*/
 #include "sparse.hxx"
 #include "context.hxx"
 #include "gatewaystruct.hxx"
@@ -48,8 +48,8 @@ SciErr getComplexSparseMatrix(void* _pvCtx, int* _piAddress, int* _piRows, int*
 SciErr getCommonSparseMatrix(void* _pvCtx, int* _piAddress, int _iComplex, int* _piRows, int* _piCols, int* _piNbItem, int** _piNbItemRow, int** _piColPos, double** _pdblReal, double** _pdblImg)
 {
     SciErr sciErr = sciErrInit();
-    int iPos   = 0;
-    int iType  = 0;
+    int iPos = 0;
+    int iType = 0;
 
     if (_piAddress == NULL)
     {
@@ -145,7 +145,7 @@ SciErr allocCommonSparseMatrix(void* _pvCtx, int _iVar, int _iComplex, int _iRow
     return sciErr;
 }
 
-SciErr fillCommonSparseMatrix(void* _pvCtx, int *_piAddress, int _iComplex, int _iRows, int _iCols, int _iNbItem, const int* _piNbItemRow, const int* _piColPos, const double* _pdblReal, const double* _pdblImg, int* _piTotalSize)
+SciErr fillCommonSparseMatrix(void* _pvCtx, int **_piAddress, int _iComplex, int _iRows, int _iCols, int _iNbItem, const int* _piNbItemRow, const int* _piColPos, const double* _pdblReal, const double* _pdblImg, int* _piTotalSize)
 {
     SciErr sciErr = sciErrInit();
 
@@ -155,32 +155,54 @@ SciErr fillCommonSparseMatrix(void* _pvCtx, int *_piAddress, int _iComplex, int
         return sciErr;
     }
 
-    Sparse* pSparse = (Sparse*)_piAddress;
+    //convert to ij, val, dims format to call sparse constructor
 
+    //dims
+    Double* dims = new Double(1, 2, false);
+    dims->get()[0] = (double)_iRows;
+    dims->get()[1] = (double)_iCols;
+
+    //ij
+    Double* ij = new Double(_iNbItem, 2);
+    double* pI = ij->get();
+    double* pJ = ij->get() + _iNbItem;
+
+    int idx = 0;
+    for (int i = 0; i < _iRows; i++)
+    {
+        for (int j = 0; j < _piNbItemRow[i]; j++)
+        {
+            pI[idx] = i + 1;
+            pJ[idx] = *_piColPos++;
+            ++idx;
+        }
+    }
+
+    Double* val = new Double(_iNbItem, 1, _iComplex == 1);
+    double* pR = val->get();
     if (_iComplex)
     {
-        for (int i = 0; i < _iRows; i++)
+        double* pI = val->getImg();
+        for (int i = 0; i < _iNbItem; ++i)
         {
-            for (int j = 0; j < _piNbItemRow[i]; j++)
-            {
-                int iIndex = (*_piColPos++ - 1) * _iRows + i;
-                std::complex<double> cplx(*_pdblReal++, *_pdblImg++);
-                pSparse->set(iIndex, cplx);
-            }
+            pR[i] = _pdblReal[i];
+            pI[i] = _pdblImg[i];
         }
     }
     else
     {
-        for (int i = 0; i < _iRows; i++)
+        for (int i = 0; i < _iNbItem; ++i)
         {
-            for (int j = 0; j < _piNbItemRow[i]; j++)
-            {
-                int iIndex = (*_piColPos++ - 1) * _iRows + i;
-                pSparse->set(iIndex, *_pdblReal++);
-            }
+            pR[i] = _pdblReal[i];
         }
     }
 
+    Sparse* pSparse = new Sparse(*val, *ij, *dims);
+    delete dims;
+    delete val;
+    delete ij;
+
+    *_piAddress = (int*)pSparse;
     *_piTotalSize = (int)pSparse->nonZeros();
 
     return sciErr;
@@ -214,7 +236,9 @@ SciErr createCommonSparseMatrix(void* _pvCtx, int _iVar, int _iComplex, int _iRo
     GatewayStruct* pStr = (GatewayStruct*)_pvCtx;
     InternalType** out = pStr->m_pOut;
 
-    types::Sparse* pSparse = new Sparse(_iRows, _iCols, _iComplex == 1);
+    int iTotalSize = 0;
+    types::Sparse* pSparse = NULL;
+    sciErr = fillCommonSparseMatrix(_pvCtx, (int**)&pSparse, _iComplex, _iRows, _iCols, _iNbItem, _piNbItemRow, _piColPos, _pdblReal, _pdblImg, &iTotalSize);
     if (pSparse == NULL)
     {
         addErrorMessage(&sciErr, API_ERROR_CREATE_NAMED_SPARSE, _("%s: Unable to create variable in Scilab memory"), _iComplex ? "createComplexSparseMatrix" : "createSparseMatrix");
@@ -224,8 +248,6 @@ SciErr createCommonSparseMatrix(void* _pvCtx, int _iVar, int _iComplex, int _iRo
     int rhs = _iVar - *getNbInputArgument(_pvCtx);
     out[rhs - 1] = pSparse;
 
-    int iTotalSize = 0;
-    sciErr = fillCommonSparseMatrix(_pvCtx, (int*)pSparse, _iComplex, _iRows, _iCols, _iNbItem, _piNbItemRow, _piColPos, _pdblReal, _pdblImg, &iTotalSize);
     return sciErr;
 }
 
@@ -268,15 +290,15 @@ SciErr createCommonNamedSparseMatrix(void* _pvCtx, const char* _pstName, int _iC
         return sciErr;
     }
 
-    types::Sparse* pSparse = new Sparse(_iRows, _iCols, _iComplex == 1);
+    int iTotalSize = 0;
+    types::Sparse* pSparse = NULL;
+    sciErr = fillCommonSparseMatrix(_pvCtx, (int**)&pSparse, _iComplex, _iRows, _iCols, _iNbItem, _piNbItemRow, _piColPos, _pdblReal, _pdblImg, &iTotalSize);
     if (pSparse == NULL)
     {
         addErrorMessage(&sciErr, API_ERROR_CREATE_NAMED_SPARSE, _("%s: Unable to create %s named \"%s\""), _iComplex ? "createNamedComplexSparseMatrix" : "createNamedSparseMatrix", _("sparse matrix"), _pstName);
         return sciErr;
     }
 
-    int iTotalSize = 0;
-    sciErr = fillCommonSparseMatrix(_pvCtx, (int*)pSparse, _iComplex, _iRows, _iCols, _iNbItem, _piNbItemRow, _piColPos, _pdblReal, _pdblImg, &iTotalSize);
 
     wchar_t* pwstName = to_wide_string(_pstName);
     symbol::Context::getInstance()->put(symbol::Symbol(pwstName), pSparse);
@@ -301,8 +323,8 @@ SciErr readCommonNamedSparseMatrix(void* _pvCtx, const char* _pstName, int _iCom
     int* piColPos = 0;
     int iOne = 1;
 
-    double* pdblReal   = NULL;
-    double* pdblImg            = NULL;
+    double* pdblReal = NULL;
+    double* pdblImg = NULL;
 
     SciErr sciErr = getVarAddressFromName(_pvCtx, _pstName, &piAddr);
     if (sciErr.iErr)
index 538de89..d62c2da 100644 (file)
@@ -520,6 +520,7 @@ private :
      */
     template<typename DestIter>
     void create(int rows, int cols, Double SPARSE_CONST& src, DestIter o, std::size_t n);
+    void create2(int rows, int cols, Double SPARSE_CONST& src, Double SPARSE_CONST& idx);
 
     /** utility function used by insert functions conceptually : sp[destTrav]= src[srcTrav]
         @param src : data source
@@ -724,9 +725,7 @@ struct EXTERN_AST SparseBool : GenericType
     BoolSparse_t* matrixBool;
 
 private:
-    template<typename DestIter>
-    void create(int rows, int cols, Bool SPARSE_CONST& src, DestIter o, std::size_t n);
-
+    void create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx);
 };
 template<typename T>
 struct SparseTraits
index 53e8162..8a6909e 100644 (file)
@@ -106,7 +106,7 @@ InternalType *GenericRDivide(InternalType *_pLeftOperand, InternalType *_pRightO
                     sciprint(_("Warning : Division by zero...\n"));
                 }
                 break;
-                //            default : throw ast::ScilabError(_W("Operator / : Error %d not yet managed.\n"), iResult);
+            //            default : throw ast::ScilabError(_W("Operator / : Error %d not yet managed.\n"), iResult);
             default :
                 sciprint(_("Operator / : Error %d not yet managed.\n"), iResult);
         }
@@ -459,18 +459,18 @@ int RDivideSparseByDouble(types::Sparse* _pSp, types::Double* _pDouble, Internal
 
     size_t iSize = _pSp->nonZeros();
     int* Col = new int[iSize];
-    int* Row = new int[iSize];
+    int* Row = new int[_pSp->getRows()];
     _pSp->getColPos(Col);
     _pSp->getNbItemByRow(Row);
     int* iPositVal = new int[iSize];
 
-    int j = 0;
-    for (int i = 0; i < iSize; j++)
+    int idx = 0;
+    for (int i = 0; i < _pSp->getRows(); i++)
     {
-        for (int k = 0; k < Row[j]; k++)
+        for (int j = 0; j < Row[i]; j++)
         {
-            iPositVal[i] = (Col[i] - 1) * _pSp->getRows() + j;
-            i++;
+            iPositVal[idx] = (Col[idx] - 1) * _pSp->getRows() + i;
+            ++idx;
         }
     }
 
index 8b55f1b..0799e76 100644 (file)
@@ -299,22 +299,32 @@ Sparse::Sparse(int _iRows, int _iCols, bool cplx)
 
 Sparse::Sparse(Double SPARSE_CONST& src)
 {
-    create(src.getRows(), src.getCols(), src, RowWiseFullIterator(src.getRows(), src.getCols()), src.getSize());
+    //compute idx
+    int size = src.getSize();
+    int row = src.getRows();
+    Double* idx = new Double(src.getSize(), 2);
+    double* p = idx->get();
+    for (int i = 0; i < size; ++i)
+    {
+        p[i]        = (double)(i % row) + 1;
+        p[i + size] = (double)(i / row) + 1;
+    }
+    create2(src.getRows(), src.getCols(), src, *idx);
+    delete idx;
 }
 
 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
 {
-    double SPARSE_CONST* const endOfRow(idx.getReal() + idx.getRows());
-    create( static_cast<int>(*std::max_element(idx.getReal(), endOfRow))
-            , static_cast<int>(*std::max_element(endOfRow, endOfRow + idx.getRows()))
-            , src, makeIteratorFromVar(idx), idx.getSize() / 2 );
+    int idxrow = idx.getRows();
+    int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
+    int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
+
+    create2(rows, cols, src, idx);
 }
 
 Sparse::Sparse(Double SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
 {
-    create(static_cast<int>(dims.getReal(0, 0))
-           , static_cast<int>(dims.getReal(0, 1))
-           , src, makeIteratorFromVar(idx), idx.getSize() / 2);
+    create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
 }
 
 Sparse::Sparse(RealSparse_t* realSp, CplxSparse_t* cplxSp):  matrixReal(realSp), matrixCplx(cplxSp)
@@ -337,7 +347,7 @@ Sparse::Sparse(RealSparse_t* realSp, CplxSparse_t* cplxSp):  matrixReal(realSp),
 
 Sparse::Sparse(Double SPARSE_CONST& xadj, Double SPARSE_CONST& adjncy, Double SPARSE_CONST& src, std::size_t r, std::size_t c)
 {
-    Adjacency a(xadj.getReal(), adjncy.getReal());
+    Adjacency a(xadj.get(), adjncy.get());
     create(static_cast<int>(r), static_cast<int>(c), src, makeIteratorFromVar(a), src.getSize());
 }
 
@@ -354,7 +364,7 @@ void Sparse::create(int rows, int cols, Double SPARSE_CONST& src, DestIter o, st
     if (src.isComplex())
     {
         matrixReal = 0;
-        matrixCplx =  new CplxSparse_t( rows, cols);
+        matrixCplx = new CplxSparse_t(rows, cols);
         matrixCplx->reserve((int)n);
         mycopy_n(makeMatrixIterator<std::complex<double> >(src, RowWiseFullIterator(src.getRows(), src.getCols())), n, makeMatrixIterator<std::complex<double> >(*matrixCplx, o));
     }
@@ -362,13 +372,67 @@ void Sparse::create(int rows, int cols, Double SPARSE_CONST& src, DestIter o, st
     {
         matrixReal = new RealSparse_t(rows, cols);
         matrixReal->reserve((int)n);
-        matrixCplx =  0;
-        mycopy_n(makeMatrixIterator<double >(src,  RowWiseFullIterator(src.getRows(), src.getCols())), n
+        matrixCplx = 0;
+        mycopy_n(makeMatrixIterator<double >(src, RowWiseFullIterator(src.getRows(), src.getCols())), n
                  , makeMatrixIterator<double>(*matrixReal, o));
     }
     finalize();
 }
 
+void Sparse::create2(int rows, int cols, Double SPARSE_CONST& src, Double SPARSE_CONST& idx)
+{
+    int nnz = src.getSize();
+    double* i = idx.get();
+    double* j = i + idx.getRows();
+    double* valR = src.get();
+
+    if (src.isComplex())
+    {
+        matrixReal = 0;
+
+        typedef Eigen::Triplet<std::complex<double> > T;
+        std::vector<T> tripletList;
+        tripletList.reserve((int)nnz);
+
+        double* valI = src.getImg();
+
+        for (int k = 0; k < nnz; ++k)
+        {
+            tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, std::complex<double>(valR[k], valI[k])));
+        }
+
+        matrixCplx = new CplxSparse_t(rows, cols);
+        matrixCplx->setFromTriplets(tripletList.begin(), tripletList.end());
+        m_iRows = matrixCplx->rows();
+        m_iCols = matrixCplx->cols();
+    }
+    else
+    {
+        matrixCplx = 0;
+
+        typedef Eigen::Triplet<double> T;
+        std::vector<T> tripletList;
+        tripletList.reserve((int)nnz);
+
+        for (int k = 0; k < nnz; ++k)
+        {
+            tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, valR[k]));
+        }
+
+        matrixReal = new RealSparse_t(rows, cols);
+        matrixReal->setFromTriplets(tripletList.begin(), tripletList.end());
+
+        m_iRows = matrixReal->rows();
+        m_iCols = matrixReal->cols();
+    }
+
+    m_iSize = m_iCols * m_iRows;
+    m_iDims = 2;
+    m_piDims[0] = m_iRows;
+    m_piDims[1] = m_iCols;
+    finalize();
+}
+
 void Sparse::fill(Double& dest, int r, int c) SPARSE_CONST
 {
     Sparse & cthis(const_cast<Sparse&>(*this));
@@ -1562,7 +1626,7 @@ InternalType* Sparse::extract(typed_list* _pArgs)
         }
     }
 
-    finalize();
+    pOut->finalize();
 
     //free pArg content
     cleanIndexesArguments(_pArgs, &pArg);
@@ -1884,7 +1948,7 @@ template<typename S> struct GetReal: std::unary_function<typename S::InnerIterat
     }
 };
 template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
-        : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
+    : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
 {
     double operator()( Sparse::CplxSparse_t::InnerIterator it) const
     {
@@ -2119,17 +2183,27 @@ bool Sparse::reshape(int _iNewRows, int _iNewCols)
 
 SparseBool::SparseBool(Bool SPARSE_CONST& src)
 {
-    create(src.getRows(), src.getCols(), src, RowWiseFullIterator(src.getRows(), src.getCols()), src.getSize());
+    //compute idx
+    int size = src.getSize();
+    int row = src.getRows();
+    Double* idx = new Double(src.getSize(), 2);
+    double* p = idx->get();
+    for (int i = 0; i < size; ++i)
+    {
+        p[i] = (double)(i % row) + 1;
+        p[i + size] = (double)(i / row) + 1;
+    }
+    create2(src.getRows(), src.getCols(), src, *idx);
 }
 /* @param src : Bool matrix to copy into a new sparse matrix
 @param idx : Double matrix to use as indexes to get values from the src
 **/
 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
 {
-    double SPARSE_CONST* const endOfRow(idx.getReal() + idx.getRows());
-    create( static_cast<int>(*std::max_element(idx.getReal(), endOfRow))
-            , static_cast<int>(*std::max_element(endOfRow, endOfRow + idx.getRows()))
-            , src, makeIteratorFromVar(idx), idx.getSize() / 2 );
+    int idxrow = idx.getRows();
+    int rows = static_cast<int>(*std::max_element(idx.get(), idx.get() + idxrow));
+    int cols = static_cast<int>(*std::max_element(idx.get() + idxrow, idx.get() + idxrow * 2));
+    create2(rows, cols, src, idx);
 }
 
 /* @param src : Bool matrix to copy into a new sparse matrix
@@ -2138,7 +2212,7 @@ SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
 **/
 SparseBool::SparseBool(Bool SPARSE_CONST& src, Double SPARSE_CONST& idx, Double SPARSE_CONST& dims)
 {
-    create((int)dims.getReal(0, 0) , (int)dims.getReal(0, 1) , src, makeIteratorFromVar(idx), (int)idx.getSize() / 2);
+    create2(static_cast<int>(dims.get(0)), static_cast<int>(dims.get(1)), src, idx);
 }
 
 SparseBool::SparseBool(int _iRows, int _iCols) : matrixBool(new BoolSparse_t(_iRows, _iCols))
@@ -2171,25 +2245,34 @@ SparseBool::SparseBool(BoolSparse_t* src) : matrixBool(src)
     m_piDims[1] = m_iCols;
 }
 
-template<typename DestIter>
-void SparseBool::create(int rows, int cols, Bool SPARSE_CONST& src, DestIter o, std::size_t n)
+void SparseBool::create2(int rows, int cols, Bool SPARSE_CONST& src, Double SPARSE_CONST& idx)
 {
-    m_iCols = cols;
-    m_iRows = rows;
+    int nnz = src.getSize();
+    double* i = idx.get();
+    double* j = i + idx.getRows();
+    int* val = src.get();
+
+    typedef Eigen::Triplet<bool> T;
+    std::vector<T> tripletList;
+    tripletList.reserve((int)nnz);
+
+    for (int k = 0; k < nnz; ++k)
+    {
+        tripletList.push_back(T(static_cast<int>(i[k]) - 1, static_cast<int>(j[k]) - 1, val[k] == 1));
+    }
+
+    matrixBool = new BoolSparse_t(rows, cols);
+    matrixBool->setFromTriplets(tripletList.begin(), tripletList.end());
+
+    m_iRows = matrixBool->rows();
+    m_iCols = matrixBool->cols();
     m_iSize = cols * rows;
     m_iDims = 2;
     m_piDims[0] = m_iRows;
     m_piDims[1] = m_iCols;
-
-    matrixBool = new BoolSparse_t(rows, cols);
-
-    matrixBool->reserve((int)n);
-    mycopy_n(makeMatrixIterator<int>(src,  RowWiseFullIterator(src.getRows(), src.getCols())), n
-             , makeMatrixIterator<bool>(*matrixBool, o));
     finalize();
 }
 
-
 bool SparseBool::toString(std::wostringstream& ostr) const
 {
     ostr << ::toString(*matrixBool, 0);
index da2c332..8703743 100644 (file)
@@ -23,73 +23,6 @@ extern "C"
 
 using namespace types;
 
-namespace
-{
-template<typename DenseType>
-GenericType* buildSparse(Double SPARSE_CONST& indices, DenseType SPARSE_CONST& vals, Double SPARSE_CONST* pDim)
-{
-    GenericType* pRes = NULL;
-    InternalType* pVals = (InternalType*)&vals;
-    if (indices.isEmpty() && vals.isDouble() && pVals->getAs<Double>()->isEmpty())
-    {
-        //sparse([],[], ...)
-        if (pDim == NULL)
-        {
-            //sparse([], []) -> KO
-            return NULL;
-        }
-
-        //check pDims is a row vector(1,2)
-        if (pDim->getRows() != 1 || pDim->getCols() != 2)
-        {
-            return NULL;
-        }
-
-        //can create an empty sparse matrix
-        return new Sparse((int)pDim->get(0), (int)pDim->get(1));
-    }
-
-    if (indices.getRows() != vals.getSize() || (indices.getCols() != 2 && indices.getCols() != 0))
-    {
-        return NULL;
-    }
-
-    // indices must be >=1
-    if (*std::min_element(indices.getReal(), indices.getReal() + indices.getSize()) < 1)
-    {
-        return NULL;
-    }
-
-    if (pDim == NULL)
-    {
-        pRes = new typename types::SparseTraits<DenseType>::type(vals, indices);
-    }
-    else
-    {
-        // if three arguments are given the first two are the same for the case of two arugments and the third is a '1 x 2' vector for the size of the matrix
-        if (pDim->getRows() != 1 || pDim->getCols() != 2)
-        {
-            return NULL;
-        }
-
-        double* endOfRow(indices.getReal() + indices.getRows());
-        if (*std::max_element(indices.getReal(), endOfRow) > pDim->getReal(0, 0))
-        {
-            return NULL;
-        }
-
-        if (*std::max_element(endOfRow, endOfRow + indices.getRows()) > pDim->getReal(0, 1))
-        {
-            return NULL;
-        }
-
-        pRes = new typename types::SparseTraits<DenseType>::type(vals, indices, *pDim);
-    }
-    return pRes;
-}
-}
-
-
 Function::ReturnValue sci_sparse(typed_list &in, int _piRetCount, typed_list &out)
 {
     bool isValid = true;
@@ -201,22 +134,40 @@ Function::ReturnValue sci_sparse(typed_list &in, int _piRetCount, typed_list &ou
             }
         }
 
-        types::GenericType* pGT1 = in[0]->getAs<types::GenericType>();
+        Double* ij = in[0]->getAs<Double>();
         types::GenericType* pGT2 = in[1]->getAs<types::GenericType>();
 
-        if (pGT2->getSize() != pGT1->getRows())
+        if (pGT2->getSize() != ij->getRows())
         {
-            Scierror(999, _("%s: Wrong size for input argument #%d: A matrix of size %d expected.\n"), "sparse", 2, pGT1->getRows());
+            Scierror(999, _("%s: Wrong size for input argument #%d: A matrix of size %d expected.\n"), "sparse", 2, ij->getRows());
             return Function::Error;
         }
 
+        bool alloc = false;
+        if (pDims == nullptr)
+        {
+            int size = ij->getRows();
+            double* i = ij->get();
+            double* j = i + ij->getRows();
+            pDims = new Double(1, 2, false);
+            pDims->set(0, *std::max_element(i, i + size));
+            pDims->set(1, *std::max_element(j, j + size));
+        }
+
         if (in[1]->isDouble())
         {
-            pRetVal = buildSparse(*in[0]->getAs<Double>(), *in[1]->getAs<Double>(), pDims);
+            Double* dbl = pGT2->getAs<Double>();
+            pRetVal = new Sparse(*dbl, *ij, *pDims);
         }
         else
         {
-            pRetVal = buildSparse(*in[0]->getAs<Double>(), *in[1]->getAs<Bool>(), pDims);
+            Bool* b = pGT2->getAs<Bool>();
+            pRetVal = new SparseBool(*b, *ij, *pDims);
+        }
+
+        if (alloc)
+        {
+            delete pDims;
         }
     }