linear_algebra plugged.
[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 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution.  The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 /*--------------------------------------------------------------------------*/
13
14 #include "linear_algebra_gw.hxx"
15 #include "function.hxx"
16 #include "double.hxx"
17 #include "overload.hxx"
18 #include "execvisitor.hxx"
19
20 extern "C"
21 {
22 #include "localization.h"
23 #include "Scierror.h"
24 #include "lsq.h"
25 }
26 /*--------------------------------------------------------------------------*/
27
28 types::Function::ReturnValue sci_lsq(types::typed_list &in, int _iRetCount, types::typed_list &out)
29 {
30     types::Double* pDbl[2]      = {NULL, NULL};
31     types::Double* pDblResult   = NULL;
32     double* pData[2]            = {NULL, NULL};
33     double* pResult             = NULL;
34     double* pdTol               = NULL;
35     bool bComplexArgs           = false;
36     int iRank                   = 0;
37
38     if(in.size() < 2 || in.size() > 3)
39     {
40         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d to %d expected.\n"), L"lsq", 2, 3);
41         return types::Function::Error;
42     }
43
44     if(_iRetCount > 2)
45     {
46         ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d to %d expected.\n"), L"lsq", 1, 2);
47         return types::Function::Error;
48     }
49
50     if((in[0]->isDouble() == false))
51     {
52         std::wstring wstFuncName = L"%"  + in[0]->getShortTypeStr() + L"_lsq";
53         return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
54     }
55
56     if(in.size() == 2)
57     {
58         if((in[1]->isDouble() == false))
59         {
60             std::wstring wstFuncName = L"%"  + in[1]->getShortTypeStr() + L"_lsq";
61             return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
62         }
63         pDbl[1] = in[1]->getAs<types::Double>()->clone()->getAs<types::Double>();
64     }
65
66     if(in.size() == 3)
67     {
68         if((in[2]->isDouble() == false) || (in[2]->getAs<types::Double>()->isComplex()) || (in[2]->getAs<types::Double>()->isScalar() == false))
69         {
70             ScierrorW(256, _W("%ls: Wrong type for input argument #%d: A Real expected.\n"), L"lsq", 3);
71             return types::Function::Error;     
72         }
73         *pdTol = in[2]->getAs<types::Double>()->get(0);
74     }
75
76     pDbl[0] = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>();
77
78     if(pDbl[0]->getRows() != pDbl[1]->getRows())
79     {
80         ScierrorW(265, _W("%ls: %ls and %ls must have equal number of rows.\n"), L"lsq", L"A", L"B");
81         return types::Function::Error;
82     }
83
84     if((pDbl[0]->getCols() == 0) || (pDbl[1]->getCols() == 0))
85     {
86         out.push_back(types::Double::Empty());
87         if(_iRetCount == 2)
88         {
89             out.push_back(types::Double::Empty());
90         }
91         return types::Function::OK;
92     }
93
94     if(pDbl[0]->isComplex() || pDbl[1]->isComplex())
95     {
96         bComplexArgs = true;
97     }
98     for(int i=0; i<2; i++)
99     {
100         if(pDbl[i]->getCols() == -1)
101         {
102             ScierrorW(271, _W("%ls: Size varying argument a*eye(), (arg %d) not allowed here.\n"), L"lsq", i+1);
103             return types::Function::Error;
104         }
105
106         if(bComplexArgs)
107         {
108             pData[i] = (double*)oGetDoubleComplexFromPointer(pDbl[i]->getReal(), pDbl[i]->getImg(), pDbl[i]->getSize());
109             if(!pData[i])
110             {
111                 ScierrorW(999,_W("%ls: Cannot allocate more memory.\n"),L"lsq");
112                 return types::Function::Error;
113             }
114         }
115         else
116         {
117             pData[i] = pDbl[i]->getReal();
118         }
119     }
120
121     pDblResult = new types::Double(pDbl[0]->getCols(), pDbl[1]->getCols(), bComplexArgs);
122
123     if(bComplexArgs)
124     {
125         pResult = (double*)MALLOC(pDbl[0]->getCols() * pDbl[1]->getCols() * sizeof(doublecomplex));
126     }
127     else
128     {
129         pResult = pDblResult->get();
130     }
131
132     int iRet = iLsqM(pData[0], pDbl[0]->getRows(), pDbl[0]->getCols(), pData[1], pDbl[1]->getCols(), bComplexArgs, pResult, pdTol, ((_iRetCount == 2) ? &iRank : NULL));
133
134     if(iRet != 0)
135     {
136         if(iRet == -1)
137         {
138                     ScierrorW(999, _W("%ls: Allocation failed.\n"), L"lsq");
139         }
140         else
141         {
142                     ScierrorW(999, _W("%ls: LAPACK error n°%d.\n"), L"lsq",iRet);
143         }
144         return types::Function::Error;
145     }
146
147     if(bComplexArgs)
148     {
149         vGetPointerFromDoubleComplex((doublecomplex*)(pResult), pDblResult->getSize(), pDblResult->getReal(), pDblResult->getImg());
150         vFreeDoubleComplexFromPointer((doublecomplex*)pResult);
151                 vFreeDoubleComplexFromPointer((doublecomplex*)pData[0]);
152         vFreeDoubleComplexFromPointer((doublecomplex*)pData[1]);
153     }
154
155     out.push_back(pDblResult);
156     if(_iRetCount == 2)
157     {
158         types::Double* pDblRank = new types::Double(1,1);
159         pDblRank->set(0, iRank);
160         out.push_back(pDblRank);
161     }
162
163     return types::Function::OK;
164 }
165 /*--------------------------------------------------------------------------*/
166