linear_algebra plugged.
[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 }
27 /*--------------------------------------------------------------------------*/
28
29 types::Function::ReturnValue sci_qr(types::typed_list &in, int _iRetCount, types::typed_list &out)
30 {
31     types::Double* pDbl     = NULL;
32     types::Double* pDblQ    = NULL;
33     types::Double* pDblR    = NULL;
34     types::Double* pDblE    = NULL;
35     types::Double* pDblRk   = NULL;
36     double* pdQ             = NULL;
37     double* pdR             = NULL;
38     double* pData           = NULL;
39     double dTol             = -1.;
40     int iRowsToCompute      = 0;
41
42     if(in.size() != 1 && in.size() != 2)
43     {
44         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d to %d expected.\n"), L"qr", 1, 2);
45         return types::Function::Error;
46     }
47
48     if(_iRetCount > 4)
49     {
50         ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d to %d expected.\n"), L"qr", 1, 4);
51         return types::Function::Error;
52     }
53
54     if((in[0]->isDouble() == false))
55     {
56         std::wstring wstFuncName = L"%"  + in[0]->getShortTypeStr() + L"_qr";
57         return Overload::call(wstFuncName, in, _iRetCount, out, new ExecVisitor());
58     }
59
60     pDbl = in[0]->getAs<types::Double>()->clone()->getAs<types::Double>();
61
62     if(pDbl->isComplex())
63     {
64         pData = (double*)oGetDoubleComplexFromPointer(pDbl->getReal(), pDbl->getImg(), pDbl->getSize());
65         if(!pData)
66         {
67             ScierrorW(999,_W("%ls: Cannot allocate more memory.\n"),L"qr");
68             return types::Function::Error;
69         }
70     }
71     else
72     {
73         pData = pDbl->getReal();
74     }
75
76     if((pDbl->getCols() == 0) || (pDbl->getRows() == 0))
77     {
78         if(_iRetCount == 4)
79         {
80             types::Double* zero = new types::Double(1,1);
81             zero->set(0,0);
82             out.push_back(types::Double::Empty());
83             out.push_back(types::Double::Empty());
84             out.push_back(zero);
85             out.push_back(types::Double::Empty());
86         }
87         else
88         {
89             out.push_back(types::Double::Empty());
90             for(int i=1; i < _iRetCount; i++)
91             {
92                 out.push_back(types::Double::Empty());
93             }
94         }
95         return types::Function::OK;
96     }
97
98     if((pDbl->getRows() == -1) || (pDbl->getCols() == -1)) // manage eye case
99     {
100         ScierrorW(271,_W("%ls: Size varying argument a*eye(), (arg %d) not allowed here.\n"), L"qr", 1);
101         return types::Function::Error;
102     }
103
104     iRowsToCompute = pDbl->getRows();
105
106     if(in.size() == 2)
107     {
108        if(in[1]->isString() == true)
109         {/* /!\ original code did not check that string is "e" so any [matrix of] string is accepted as "e" ! */
110             /*
111             types::String* pStr = in[1]->getAs<types::String>();
112             if((wcslen(pStr->get(0)) == 1) && (pStr->get(0)[0] == L'e'))
113             */
114             iRowsToCompute = Min(pDbl->getRows(), pDbl->getCols());
115         }
116         else if(in[1]->isDouble() == true)
117         {/* /!\ original code do not check anything (real && 1x1 matrix)*/
118             dTol = in[1]->getAs<types::Double>()->get(0);
119         }
120         else
121         {
122             ScierrorW(999,_W("%ls: Wrong type for input argument #%d: A real or a string expected.\n"), L"qr", 2);
123             return types::Function::Error;
124         }
125     }
126
127     pDblQ = new types::Double(pDbl->getRows(), iRowsToCompute, pDbl->isComplex());
128     pDblR = new types::Double(iRowsToCompute, pDbl->getCols(), pDbl->isComplex());
129
130     if(pDbl->isComplex())
131     {
132                 pdQ = (double*) MALLOC(pDbl->getRows() * iRowsToCompute * sizeof(doublecomplex));
133                 pdR = (double*) MALLOC(iRowsToCompute * pDbl->getCols() * sizeof(doublecomplex));
134     }
135     else
136     {
137         pdQ = pDblQ->get();
138         pdR = pDblR->get();
139     }
140
141     if(_iRetCount >= 3)/* next alloc for E needed only for _iRetCount>=3 */
142     {
143         pDblE = new types::Double(pDbl->getCols(), pDbl->getCols());
144     }
145
146     if(_iRetCount >= 4)/* next alloc for Rk needed only for _iRetCount>=4 */
147     {
148         pDblRk = new types::Double(1,1);
149     }
150
151     int iRet = iQrM(pData, pDbl->getRows(), pDbl->getCols(), pDbl->isComplex(), iRowsToCompute, dTol, pdQ, pdR, (pDblE ? pDblE->get() : NULL), (pDblRk ? pDblRk->get() : NULL));
152
153     if(iRet != 0)
154     {
155             ScierrorW(999, _W("%ls: LAPACK error n°%d.\n"), L"qr",iRet);
156         return types::Function::Error;
157     }
158
159     if(pDbl->isComplex())
160     {
161         vGetPointerFromDoubleComplex((doublecomplex*)pdQ, pDblQ->getSize(), pDblQ->getReal(), pDblQ->getImg());
162         vFreeDoubleComplexFromPointer((doublecomplex*)pdQ);
163
164         vGetPointerFromDoubleComplex((doublecomplex*)pdR, pDblR->getSize(), pDblR->getReal(), pDblR->getImg());
165         vFreeDoubleComplexFromPointer((doublecomplex*)pdR);
166     }
167
168     if(pDbl->isComplex())
169     {
170         vFreeDoubleComplexFromPointer((doublecomplex*)pData);
171     }
172
173     out.push_back(pDblQ);
174     if(_iRetCount >= 2)
175     {
176         out.push_back(pDblR);
177     }
178
179     if(_iRetCount == 3)
180     {
181         out.push_back(pDblE);
182     }
183     else if(_iRetCount == 4)
184     {
185         out.push_back(pDblRk);
186         out.push_back(pDblE);
187     }
188
189     return types::Function::OK;
190 }
191 /*--------------------------------------------------------------------------*/
192