835b4f7ed3c2c89aff0f17fb3afe6a8d8032f117
[scilab.git] / 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
41     if (in.size() < 2 || in.size() > 3)
42     {
43         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "lsq", 2, 3);
44         return types::Function::Error;
45     }
46
47     if (_iRetCount > 2)
48     {
49         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "lsq", 1, 2);
50         return types::Function::Error;
51     }
52
53     if ((in[0]->isDouble() == false))
54     {
55         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lsq";
56         return Overload::call(wstFuncName, in, _iRetCount, out);
57     }
58
59     pDbl[0] = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>();
60
61     if (in.size() <=  3)
62     {
63         if ((in[1]->isDouble() == false))
64         {
65             std::wstring wstFuncName = L"%" + in[1]->getShortTypeStr() + L"_lsq";
66             return Overload::call(wstFuncName, in, _iRetCount, out);
67         }
68         pDbl[1] = in[1]->getAs<types::Double>()->clone()->getAs<types::Double>();
69     }
70
71     if (in.size() == 3)
72     {
73         if ((in[2]->isDouble() == false) || (in[2]->getAs<types::Double>()->isComplex()) || (in[2]->getAs<types::Double>()->isScalar() == false))
74         {
75             Scierror(256, _("%s: Wrong type for input argument #%d: A Real expected.\n"), "lsq", 3);
76             return types::Function::Error;
77         }
78
79         double dblTol = in[2]->getAs<types::Double>()->get(0);
80         pdTol = &dblTol;
81     }
82
83     if (pDbl[0]->getRows() != pDbl[1]->getRows())
84     {
85         Scierror(265, _("%s: %s and %s must have equal number of rows.\n"), "lsq", "A", "B");
86         return types::Function::Error;
87     }
88
89     if ((pDbl[0]->getCols() == 0) || (pDbl[1]->getCols() == 0))
90     {
91         out.push_back(types::Double::Empty());
92         if (_iRetCount == 2)
93         {
94             out.push_back(types::Double::Empty());
95         }
96         return types::Function::OK;
97     }
98
99     if (pDbl[0]->isComplex() || pDbl[1]->isComplex())
100     {
101         bComplexArgs = true;
102     }
103     for (int i = 0; i < 2; i++)
104     {
105         if (pDbl[i]->getCols() == -1)
106         {
107             Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "lsq", i + 1);
108             return types::Function::Error;
109         }
110
111         if (bComplexArgs)
112         {
113             pData[i] = (double*)oGetDoubleComplexFromPointer(pDbl[i]->getReal(), pDbl[i]->getImg(), pDbl[i]->getSize());
114             if (!pData[i])
115             {
116                 Scierror(999, _("%s: Cannot allocate more memory.\n"), "lsq");
117                 return types::Function::Error;
118             }
119         }
120         else
121         {
122             pData[i] = pDbl[i]->getReal();
123         }
124     }
125
126     pDblResult = new types::Double(pDbl[0]->getCols(), pDbl[1]->getCols(), bComplexArgs);
127
128     if (bComplexArgs)
129     {
130         pResult = (double*)MALLOC(pDbl[0]->getCols() * pDbl[1]->getCols() * sizeof(doublecomplex));
131     }
132     else
133     {
134         pResult = pDblResult->get();
135     }
136
137     int iRet = iLsqM(pData[0], pDbl[0]->getRows(), pDbl[0]->getCols(), pData[1], pDbl[1]->getCols(), bComplexArgs, pResult, pdTol, ((_iRetCount == 2) ? &iRank : NULL));
138
139     if (iRet != 0)
140     {
141         if (iRet == -1)
142         {
143             Scierror(999, _("%s: Allocation failed.\n"),  "lsq");
144         }
145         else
146         {
147             Scierror(999, _("%s: LAPACK error n°%d.\n"),  "lsq", iRet);
148         }
149         return types::Function::Error;
150     }
151
152     if (bComplexArgs)
153     {
154         vGetPointerFromDoubleComplex((doublecomplex*)(pResult), pDblResult->getSize(), pDblResult->getReal(), pDblResult->getImg());
155         vFreeDoubleComplexFromPointer((doublecomplex*)pResult);
156         vFreeDoubleComplexFromPointer((doublecomplex*)pData[0]);
157         vFreeDoubleComplexFromPointer((doublecomplex*)pData[1]);
158     }
159
160     out.push_back(pDblResult);
161     if (_iRetCount == 2)
162     {
163         types::Double* pDblRank = new types::Double(1, 1);
164         pDblRank->set(0, iRank);
165         out.push_back(pDblRank);
166     }
167
168     return types::Function::OK;
169 }
170 /*--------------------------------------------------------------------------*/
171