6eb5526fabab54b3c09e2315d146498c20dc08a1
[scilab.git] / scilab / modules / linear_algebra / sci_gateway / cpp / sci_lu.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 "lu.h"
27 #include "doublecomplex.h"
28 }
29 /*--------------------------------------------------------------------------*/
30
31 types::Function::ReturnValue sci_lu(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     types::Double* pDbl     = NULL;
34     types::Double* pDblL    = NULL;
35     types::Double* pDblU    = NULL;
36     types::Double* pDblE    = NULL;
37     double* pdL             = NULL;
38     double* pdU             = NULL;
39     double* pData           = NULL;
40     int iMinRowsCols        = 0;
41
42     if (in.size() != 1)
43     {
44         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "lu", 1);
45         return types::Function::Error;
46     }
47
48     if (_iRetCount < 2 || _iRetCount > 3)
49     {
50         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "lu", 2, 3);
51         return types::Function::Error;
52     }
53
54     if ((in[0]->isDouble() == false))
55     {
56         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lu";
57         return Overload::call(wstFuncName, in, _iRetCount, out);
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             Scierror(999, _("%s: Cannot allocate more memory.\n"), "lu");
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         out.push_back(types::Double::Empty());
79         out.push_back(types::Double::Empty());
80         if (_iRetCount == 3)
81         {
82             out.push_back(types::Double::Empty());
83         }
84         return types::Function::OK;
85     }
86
87     if ((pDbl->getRows() == -1) || (pDbl->getCols() == -1)) // Rhs(1)=k*eye() => Lhs(1)=eye()
88     {
89         // Lhs(2)=k*eye(), Lhs(3)=eye()
90         pDblL = new types::Double(-1, -1, pDbl->isComplex());
91         pDblL->set(0, 1);
92         out.push_back(pDblL);
93         out.push_back(pDbl);
94         if (_iRetCount == 3)
95         {
96             pDblE = new types::Double(-1, -1, pDbl->isComplex());
97             pDblE->set(0, 1);
98             out.push_back(pDblE);
99         }
100         return types::Function::OK;
101     }
102
103     iMinRowsCols = std::min(pDbl->getRows(), pDbl->getCols());
104     pDblL = new types::Double(pDbl->getRows(), iMinRowsCols, pDbl->isComplex());
105     pDblU = new types::Double(iMinRowsCols, pDbl->getCols(), pDbl->isComplex());
106
107     if (pDbl->isComplex())
108     {
109         pdL = (double*) MALLOC(pDblL->getSize() * sizeof(doublecomplex));
110         pdU = (double*) MALLOC(pDblU->getSize() * sizeof(doublecomplex));
111     }
112     else
113     {
114         pdL = pDblL->get();
115         pdU = pDblU->get();
116     }
117
118     if (_iRetCount == 3)
119     {
120         pDblE = new types::Double(pDbl->getRows(), pDbl->getRows());
121     }
122
123     int iRet = iLuM(pData, pDbl->getRows(), pDbl->getCols(), pDbl->isComplex(), pdL, pdU, (pDblE ? pDblE->get() : NULL));
124
125     if (iRet != 0)
126     {
127         Scierror(999, _("%s: LAPACK error n°%d.\n"), "lu", iRet);
128         FREE((doublecomplex*)pdL);
129         FREE((doublecomplex*)pdU);
130         delete pDblL;
131         delete pDblU;
132         delete pDblE;
133         return types::Function::Error;
134     }
135
136     if (pDbl->isComplex())
137     {
138         vGetPointerFromDoubleComplex((doublecomplex*)pdL, pDblL->getSize(), pDblL->getReal(), pDblL->getImg());
139         FREE((doublecomplex*)pdL);
140         vGetPointerFromDoubleComplex((doublecomplex*)pdU, pDblU->getSize(), pDblU->getReal(), pDblU->getImg());
141         FREE((doublecomplex*)pdU);
142     }
143
144     if (pDbl->isComplex())
145     {
146         vFreeDoubleComplexFromPointer((doublecomplex*)pData);
147     }
148
149     out.push_back(pDblL);
150     out.push_back(pDblU);
151     if (_iRetCount == 3)
152     {
153         out.push_back(pDblE);
154     }
155
156     return types::Function::OK;
157 }
158 /*--------------------------------------------------------------------------*/
159