scilab/modules/linear_algebra/sci_gateway/cpp/sci_lsq.cpp

   1 /*
   2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
   4 *
   5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
   6  *
   7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
   8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
   9  * This file was originally licensed under the terms of the CeCILL v2.1,
  10  * and continues to be available under such terms.
  11  * For more information, see the COPYING file which you should have received
  12  * along with this program.
  13 *
  14 */
  15 /*--------------------------------------------------------------------------*/
  16
  17 #include "linear_algebra_gw.hxx"
  18 #include "function.hxx"
  19 #include "double.hxx"
  20 #include "overload.hxx"
  21
  22 extern "C"
  23 {
  24 #include "localization.h"
  25 #include "Scierror.h"
  26 #include "lsq.h"
  27 #include "doublecomplex.h"
  28 }
  29 /*--------------------------------------------------------------------------*/
  30
  31 types::Function::ReturnValue sci_lsq(types::typed_list &in, int _iRetCount, types::typed_list &out)
  32 {
  33     types::Double* pDbl[2]      = {NULL, NULL};
  34     types::Double* pDblResult   = NULL;
  35     double* pData[2]            = {NULL, NULL};
  36     double* pResult             = NULL;
  37     double* pdTol               = NULL;
  38     bool bComplexArgs           = false;
  39     int iRank                   = 0;
  40     double dblTol               = 0.0;
  41
  42     if (in.size() < 2 || in.size() > 3)
  43     {
  44         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "lsq", 2, 3);
  45         return types::Function::Error;
  46     }
  47
  48     if (_iRetCount > 2)
  49     {
  50         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "lsq", 1, 2);
  51         return types::Function::Error;
  52     }
  53
  54     if ((in[0]->isDouble() == false))
  55     {
  56         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lsq";
  57         return Overload::call(wstFuncName, in, _iRetCount, out);
  58     }
  59
  60     if (in.size() <=  3)
  61     {
  62         if ((in[1]->isDouble() == false))
  63         {
  64             std::wstring wstFuncName = L"%" + in[1]->getShortTypeStr() + L"_lsq";
  65             return Overload::call(wstFuncName, in, _iRetCount, out);
  66         }
  67     }
  68
  69     if (in.size() == 3)
  70     {
  71         if ((in[2]->isDouble() == false) || (in[2]->getAs<types::Double>()->isComplex()) || (in[2]->getAs<types::Double>()->isScalar() == false))
  72         {
  73             Scierror(256, _("%s: Wrong type for input argument #%d: A Real expected.\n"), "lsq", 3);
  74             return types::Function::Error;
  75         }
  76
  77         dblTol = in[2]->getAs<types::Double>()->get(0);
  78         pdTol = &dblTol;
  79     }
  80
  81     pDbl[0] = in[0]->getAs<types::Double>();
  82     pDbl[1] = in[1]->getAs<types::Double>();
  83
  84     if (pDbl[0]->getRows() != pDbl[1]->getRows())
  85     {
  86         Scierror(265, _("%s: %s and %s must have equal number of rows.\n"), "lsq", "A", "B");
  87         return types::Function::Error;
  88     }
  89
  90     if ((pDbl[0]->getCols() == 0) || (pDbl[1]->getCols() == 0))
  91     {
  92         out.push_back(types::Double::Empty());
  93         if (_iRetCount == 2)
  94         {
  95             out.push_back(types::Double::Empty());
  96         }
  97         return types::Function::OK;
  98     }
  99
 100     if (pDbl[0]->isComplex() || pDbl[1]->isComplex())
 101     {
 102         bComplexArgs = true;
 103     }
 104
 105     pDbl[0] = pDbl[0]->clone();
 106     pDbl[1] = pDbl[1]->clone();
 107
 108     for (int i = 0; i < 2; i++)
 109     {
 110         if (pDbl[i]->getCols() == -1)
 111         {
 112             for (int j=0; j<=i; j++)
 113             {
 114                 pDbl[j]->killMe();
 115             }
 116             Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "lsq", i + 1);
 117             return types::Function::Error;
 118         }
 119
 120         if (bComplexArgs)
 121         {
 122             pData[i] = (double*)oGetDoubleComplexFromPointer(pDbl[i]->getReal(), pDbl[i]->getImg(), pDbl[i]->getSize());
 123             if (!pData[i])
 124             {
 125                 for (int j=0; j<=i; j++)
 126                 {
 127                     pDbl[j]->killMe();
 128                 }
 129                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "lsq");
 130                 return types::Function::Error;
 131             }
 132         }
 133         else
 134         {
 135             pData[i] = pDbl[i]->getReal();
 136         }
 137     }
 138
 139     pDblResult = new types::Double(pDbl[0]->getCols(), pDbl[1]->getCols(), bComplexArgs);
 140
 141     if (bComplexArgs)
 142     {
 143         pResult = (double*)MALLOC(pDbl[0]->getCols() * pDbl[1]->getCols() * sizeof(doublecomplex));
 144     }
 145     else
 146     {
 147         pResult = pDblResult->get();
 148     }
 149
 150     int iRet = iLsqM(pData[0], pDbl[0]->getRows(), pDbl[0]->getCols(), pData[1], pDbl[1]->getCols(), bComplexArgs, pResult, pdTol, ((_iRetCount == 2) ? &iRank : NULL));
 151
 152     pDbl[0]->killMe();
 153     pDbl[1]->killMe();
 154
 155     if (iRet != 0)
 156     {
 157         if (iRet == -1)
 158         {
 159             Scierror(999, _("%s: Allocation failed.\n"),  "lsq");
 160         }
 161         else
 162         {
 163             Scierror(999, _("%s: LAPACK error n°%d.\n"),  "lsq", iRet);
 164         }
 165
 166         if (bComplexArgs)
 167         {
 168             vFreeDoubleComplexFromPointer((doublecomplex*)pResult);
 169             vFreeDoubleComplexFromPointer((doublecomplex*)pData[0]);
 170             vFreeDoubleComplexFromPointer((doublecomplex*)pData[1]);
 171         }
 172         pDblResult->killMe();
 173         return types::Function::Error;
 174     }
 175
 176     if (bComplexArgs)
 177     {
 178         vGetPointerFromDoubleComplex((doublecomplex*)(pResult), pDblResult->getSize(), pDblResult->getReal(), pDblResult->getImg());
 179         vFreeDoubleComplexFromPointer((doublecomplex*)pResult);
 180         vFreeDoubleComplexFromPointer((doublecomplex*)pData[0]);
 181         vFreeDoubleComplexFromPointer((doublecomplex*)pData[1]);
 182     }
 183
 184     out.push_back(pDblResult);
 185     if (_iRetCount == 2)
 186     {
 187         types::Double* pDblRank = new types::Double(1, 1);
 188         pDblRank->set(0, iRank);
 189         out.push_back(pDblRank);
 190     }
 191
 192     return types::Function::OK;
 193 }
 194 /*--------------------------------------------------------------------------*/