70f2d038ddbfb97dcd90a3f0e50cf9d56c6946af
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_qr.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2009 - DIGITEO - Bernard HUGUENEY
4 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
5 *
6 * This file must be used under the terms of the CeCILL.
7 * This source file is licensed as described in the file COPYING, which
8 * you should have received as part of this distribution.  The terms
9 * are also available at
10 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 *
12 */
13 /*--------------------------------------------------------------------------*/
14
15 #include "linear_algebra_gw.hxx"
16 #include "function.hxx"
17 #include "double.hxx"
18 #include "overload.hxx"
19 #include "execvisitor.hxx"
20
21 extern "C"
22 {
23 #include "localization.h"
24 #include "Scierror.h"
25 #include "qr.h"
26 #include "doublecomplex.h"
27 }
28 /*--------------------------------------------------------------------------*/
29
30 types::Function::ReturnValue sci_qr(types::typed_list &in, int _iRetCount, types::typed_list &out)
31 {
32     types::Double* pDbl     = NULL;
33     types::Double* pDblQ    = NULL;
34     types::Double* pDblR    = NULL;
35     types::Double* pDblE    = NULL;
36     types::Double* pDblRk   = NULL;
37     double* pdQ             = NULL;
38     double* pdR             = NULL;
39     double* pData           = NULL;
40     double dTol             = -1.;
41     int iRowsToCompute      = 0;
42
43     if (in.size() != 1 && in.size() != 2)
44     {
45         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "qr", 1, 2);
46         return types::Function::Error;
47     }
48
49     if (_iRetCount > 4)
50     {
51         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "qr", 1, 4);
52         return types::Function::Error;
53     }
54
55     if ((in[0]->isDouble() == false))
56     {
57         std::wstring wstFuncName = L"%"  + in[0]->getShortTypeStr() + L"_qr";
58         return Overload::call(wstFuncName, in, _iRetCount, out, new ast::ExecVisitor());
59     }
60
61     pDbl = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>();
62
63     if (pDbl->isComplex())
64     {
65         pData = (double*)oGetDoubleComplexFromPointer(pDbl->getReal(), pDbl->getImg(), pDbl->getSize());
66         if (!pData)
67         {
68             Scierror(999, _("%s: Cannot allocate more memory.\n"), "qr");
69             return types::Function::Error;
70         }
71     }
72     else
73     {
74         pData = pDbl->getReal();
75     }
76
77     if ((pDbl->getCols() == 0) || (pDbl->getRows() == 0))
78     {
79         if (_iRetCount == 4)
80         {
81             types::Double* zero = new types::Double(1, 1);
82             zero->set(0, 0);
83             out.push_back(types::Double::Empty());
84             out.push_back(types::Double::Empty());
85             out.push_back(zero);
86             out.push_back(types::Double::Empty());
87         }
88         else
89         {
90             out.push_back(types::Double::Empty());
91             for (int i = 1; i < _iRetCount; i++)
92             {
93                 out.push_back(types::Double::Empty());
94             }
95         }
96         return types::Function::OK;
97     }
98
99     if ((pDbl->getRows() == -1) || (pDbl->getCols() == -1)) // manage eye case
100     {
101         Scierror(271, _("%s: Size varying argument a*eye(), (arg %d) not allowed here.\n"), "qr", 1);
102         return types::Function::Error;
103     }
104
105     iRowsToCompute = pDbl->getRows();
106
107     if (in.size() == 2)
108     {
109         if (in[1]->isString() == true)
110         {
111             /* /!\ original code did not check that string is "e" so any [matrix of] string is accepted as "e" ! */
112             /*
113             types::String* pStr = in[1]->getAs<types::String>();
114             if((wcslen(pStr->get(0)) == 1) && (pStr->get(0)[0] == L'e'))
115             */
116
117             if (_iRetCount == 4)
118             {
119                 Scierror(999, _("%s: Wrong type for input argument #%d: Real scalar expected.\n"), "qr", 2);
120                 return types::Function::Error;
121             }
122
123             iRowsToCompute = Min(pDbl->getRows(), pDbl->getCols());
124         }
125         else if (in[1]->isDouble() == true)
126         {
127             /* /!\ original code do not check anything (real && 1x1 matrix)*/
128             dTol = in[1]->getAs<types::Double>()->get(0);
129         }
130         else
131         {
132             Scierror(999, _("%s: Wrong type for input argument #%d: A real or a string expected.\n"), "qr", 2);
133             return types::Function::Error;
134         }
135     }
136
137     pDblQ = new types::Double(pDbl->getRows(), iRowsToCompute, pDbl->isComplex());
138     pDblR = new types::Double(iRowsToCompute, pDbl->getCols(), pDbl->isComplex());
139
140     if (pDbl->isComplex())
141     {
142         pdQ = (double*) MALLOC(pDbl->getRows() * iRowsToCompute * sizeof(doublecomplex));
143         pdR = (double*) MALLOC(iRowsToCompute * pDbl->getCols() * sizeof(doublecomplex));
144     }
145     else
146     {
147         pdQ = pDblQ->get();
148         pdR = pDblR->get();
149     }
150
151     if (_iRetCount >= 3) /* next alloc for E needed only for _iRetCount>=3 */
152     {
153         pDblE = new types::Double(pDbl->getCols(), pDbl->getCols());
154     }
155
156     if (_iRetCount >= 4) /* next alloc for Rk needed only for _iRetCount>=4 */
157     {
158         pDblRk = new types::Double(1, 1);
159     }
160
161     int iRet = iQrM(pData, pDbl->getRows(), pDbl->getCols(), pDbl->isComplex(), iRowsToCompute, dTol, pdQ, pdR, (pDblE ? pDblE->get() : NULL), (pDblRk ? pDblRk->get() : NULL));
162
163     if (iRet != 0)
164     {
165         Scierror(999, _("%s: LAPACK error n°%d.\n"), "qr", iRet);
166         return types::Function::Error;
167     }
168
169     if (pDbl->isComplex())
170     {
171         vGetPointerFromDoubleComplex((doublecomplex*)pdQ, pDblQ->getSize(), pDblQ->getReal(), pDblQ->getImg());
172         vFreeDoubleComplexFromPointer((doublecomplex*)pdQ);
173
174         vGetPointerFromDoubleComplex((doublecomplex*)pdR, pDblR->getSize(), pDblR->getReal(), pDblR->getImg());
175         vFreeDoubleComplexFromPointer((doublecomplex*)pdR);
176     }
177
178     if (pDbl->isComplex())
179     {
180         vFreeDoubleComplexFromPointer((doublecomplex*)pData);
181     }
182
183     out.push_back(pDblQ);
184     if (_iRetCount >= 2)
185     {
186         out.push_back(pDblR);
187     }
188
189     if (_iRetCount == 3)
190     {
191         out.push_back(pDblE);
192     }
193     else if (_iRetCount == 4)
194     {
195         out.push_back(pDblRk);
196         out.push_back(pDblE);
197     }
198
199     return types::Function::OK;
200 }
201 /*--------------------------------------------------------------------------*/
202