*Bug #14330 fixed - luget was really slow.
[scilab.git] / scilab / modules / sparse / sci_gateway / cpp / sci_luget.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_luget(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 n1 = 0;
39     int nl = 0;
40     int nu = 0;
41     int nonZerosP = 0;
42     int nonZerosL = 0;
43     int nonZerosU = 0;
44     int nonZerosQ = 0;
45
46     bool cplx = false;
47     const void* pData = NULL;
48
49     double* dblP = NULL;
50     double* dblL = NULL;
51     double* dblU = NULL;
52     double* dblQ = NULL;
53     int* itemsRowP = NULL;
54     int* itemsRowL = NULL;
55     int* itemsRowU = NULL;
56     int* itemsRow = NULL;
57     int* itemsRowQ = NULL;
58     int* fmatindex = NULL;
59
60     //check input parameters
61     if (in.size() != 1)
62     {
63         Scierror(999, _("%s: Wrong number of input argument(s): %d  expected.\n"), "luget", 1);
64         return types::Function::Error;
65     }
66
67     if (_iRetCount != 4)
68     {
69         Scierror(999, _("%s: Wrong number of output argument(s): %d expected.\n"), "luget", 4);
70         return types::Function::Error;
71     }
72
73     if (in[0]->isPointer() == false)
74     {
75         Scierror(999, _("%s: Wrong type for argument %d:  Handle to sparse lu factors expected.\n"), "luget", 1);
76         return types::Function::Error;
77     }
78
79     types::Pointer *pPointerIn = in[0]->getAs<types::Pointer>();
80     n1 = pPointerIn->getCols();
81     cplx = pPointerIn->isComplex();
82     pData = pPointerIn->get();
83     fmatindex = (int*)pData;
84
85     C2F(lusiz1)(fmatindex, &nl, &nu, &ierr);
86     if (ierr > 0)
87     {
88         Scierror(999, _("Wrong value for argument #%d: the lu handle is no more valid.\n"), 1);
89         return types::Function::Error;
90     }
91
92     dblP = new double[n1];
93     dblL = new double[nl];
94     dblU = new double[nu];
95     dblQ = new double[n1];
96
97     types::Sparse *p = new types::Sparse(n1, n1, cplx);
98     types::Sparse *l = new types::Sparse(n1, n1, cplx);
99     types::Sparse *u = new types::Sparse(n1, n1, cplx);
100     types::Sparse *q = new types::Sparse(n1, n1, cplx);
101
102     itemsRowP = new int[n1 + n1];
103     itemsRowL = new int[nl + n1];
104     itemsRowU = new int[nu + n1];
105     itemsRowQ = new int[n1 + n1];
106
107     C2F(luget1)(fmatindex, itemsRowP, dblP, itemsRowL, dblL, itemsRowU, dblU, itemsRowQ, dblQ, &ierr);
108
109     int tpL = n1;
110     int tpU = n1;
111
112     for (int i = 0; i < n1; i++)
113     {
114         p->set(i, itemsRowP[i + n1] - 1, dblP[i], false);
115         q->set(i, itemsRowQ[i + n1] - 1, dblQ[i], false);
116
117         for (int j = 0; j < itemsRowL[i]; j++)
118         {
119             l->set(i, itemsRowL[j + tpL] - 1, dblL[j + tpL - n1], false);
120         }
121         tpL += itemsRowL[i];
122         for (int j = 0; j < itemsRowU[i]; j++)
123         {
124             u->set(i, itemsRowU[j + tpU] - 1, dblU[j + tpU - n1], false);
125         }
126         tpU += itemsRowU[i];
127     }
128
129     p->finalize();
130     q->finalize();
131     u->finalize();
132     l->finalize();
133
134     out.push_back(p);
135     out.push_back(l);
136     out.push_back(u);
137     out.push_back(q);
138
139     delete[] dblP;
140     delete[] dblL;
141     delete[] dblU;
142     delete[] dblQ;
143     delete[] itemsRowP;
144     delete[] itemsRowL;
145     delete[] itemsRowU;
146     delete[] itemsRowQ;
147
148     return types::Function::OK;
149 }