1c87e19e565e26a20c172bfebf6df55ea0c273e4
[scilab.git] / scilab / modules / sparse / sci_gateway / cpp / sci_lusolve.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2014 - Scilab Enterprises - Anais AUBERT
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 #include "sparse_gw.hxx"
17 #include "function.hxx"
18 #include "types.hxx"
19 #include "double.hxx"
20 #include "sparse.hxx"
21 #include "pointer.hxx"
22 #include "overload.hxx"
23
24 extern "C"
25 {
26 #include "Scierror.h"
27 #include "localization.h"
28 #include "elem_common.h"
29 #include "lu.h"
30 }
31
32 types::Function::ReturnValue sci_lusolve(types::typed_list &in, int _iRetCount, types::typed_list &out)
33 {
34     double abstol   = 0;
35     double reltol   = 0.001;
36     int nrank       = 0;
37     int ierr        = 0;
38     int m1          = 0;
39     int n1          = 0;
40     int m2          = 0;
41     int n2          = 0;
42     int nonZeros    = 0;
43     bool fact       = false;
44
45     const void *pData   = NULL;
46     double* dbl         = NULL;
47     int* colPos         = NULL;
48     int* itemsRow       = NULL;
49     int *fmatindex      = NULL;
50
51     //check input parameters
52     if (in.size() > 2)
53     {
54         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "lusolve", 2);
55         return types::Function::Error;
56     }
57
58     if (_iRetCount > 1)
59     {
60         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "lusolve", 1);
61         return types::Function::Error;
62     }
63
64     if (in[0]->isPointer())
65     {
66         types::Pointer *pPointerIn = in[0]->getAs<types::Pointer>();
67         m1 = pPointerIn->getRows();
68         n1 = pPointerIn->getCols();
69         pData = pPointerIn->get();
70         fmatindex = (int*)pData;
71         fact = false;
72     }
73     else if (in[0]->isSparse())
74     {
75         types::Sparse *pSpIn = in[0]->getAs<types::Sparse>();
76         m1 = pSpIn->getRows();
77         n1 = pSpIn->getCols();
78         if (m1 != n1)
79         {
80             Scierror(77, _("%s: Wrong size for input argument #%d: Square matrix expected.\n"), "lusolve", 1);
81             return types::Function::Error;
82         }
83
84         if (pSpIn->isComplex())
85         {
86             Scierror(77, _("%s: Wrong type for argument %d: Real matrix expected.\n"), "lusolve", 1);
87             return types::Function::Error;
88         }
89
90         nonZeros = (int)pSpIn->nonZeros();
91         dbl = new double[nonZeros];
92         pSpIn->outputValues(dbl, NULL);
93         colPos = new int[nonZeros];
94         itemsRow = new int[m1];
95         pSpIn->getColPos(colPos);
96         pSpIn->getNbItemByRow(itemsRow);
97
98         fmatindex = new int[1];
99         abstol = nc_eps_machine();
100         C2F(lufact1)(dbl, itemsRow, colPos, &m1, &nonZeros, fmatindex, &abstol, &reltol, &nrank, &ierr);
101         fact = true;
102
103         delete dbl;
104         delete colPos;
105         delete itemsRow;
106     }
107     else
108     {
109         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lusolve";
110         return Overload::call(wstFuncName, in, _iRetCount, out);
111     }
112
113     if ((in[1]->isSparse() == false) && (in[1]->isDouble() == false))
114     {
115         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_lusolve";
116         return Overload::call(wstFuncName, in, _iRetCount, out);
117     }
118
119     if (in[1]->isSparse() )
120     {
121         Scierror(999, _("%s not yet implemented for full input parameter.\n"), "lusolve");
122         return types::Function::Error;
123     }
124
125     if (in[1]->isDouble() )
126     {
127         types::Double *pDblIn = in[1]->getAs<types::Double>();
128
129         m2 = pDblIn->getRows();
130         n2 = pDblIn->getCols();
131
132         if (m2 != m1)
133         {
134             Scierror(999, _("%s: Wrong size for input argument #%d: Incompatible dimensions.\n"), "lusolve", 2);
135             return types::Function::Error;
136         }
137
138         double *dbl  = pDblIn->getReal();
139         types::Double *pDblOut = new types::Double(m2, n2, pDblIn->isComplex());
140         double *oReal = pDblOut->get();
141
142         if (pDblIn->isComplex())
143         {
144             double *imag = pDblIn->getImg();
145             double *oImg = pDblOut->getImg();
146             for (int i = 0; i < n2; i++)
147             {
148                 int iPos = i * m2;
149                 C2F(lusolve1)(fmatindex, dbl + iPos, oReal + iPos, &ierr);
150                 if (ierr > 0)
151                 {
152                     Scierror(999, _("Wrong value for argument #%d: the lu handle is no more valid.\n"), 1);
153                     return types::Function::Error;
154                 }
155                 C2F(lusolve1)(fmatindex, imag + iPos, oImg + iPos, &ierr);
156                 if (ierr > 0)
157                 {
158                     Scierror(999, _("Wrong value for argument #%d: the lu handle is no more valid.\n"), 1);
159                     return types::Function::Error;
160                 }
161             }
162         }
163         else
164         {
165             for (int i = 0; i < n2; i++)
166             {
167                 C2F(lusolve1)(fmatindex, &dbl[i * m2], &oReal[i * m2], &ierr);
168                 if (ierr > 0)
169                 {
170                     Scierror(999, _("Wrong value for argument #%d: the lu handle is no more valid.\n"), 1);
171                     return types::Function::Error;
172                 }
173             }
174         }
175
176         if (fact)
177         {
178             C2F(ludel1)(fmatindex, &ierr);
179         }
180
181         out.push_back(pDblOut);
182     }
183
184     return types::Function::OK;
185 }