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