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