* Bug 15321: lu() was leaking memory
[scilab.git] / scilab / modules / linear_algebra / src / c / lu.c
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2009-2009 - DIGITEO - Bernard HUGUENEY
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 #include "machine.h"
16 #include "core_math.h"
17 #include "sci_malloc.h"
18 #include "lu.h"
19 #include "doublecomplex.h"
20
21 extern void C2F(zgetrf)(int*, int*, double*, int*, int*, int*);
22 extern void C2F(dgetrf)(int*, int*, double*, int*, int*, int*);
23 extern void C2F(dlaset)(char const uplo[1] /* "U"pper, "L"ower, or full*/, int const * piRows, int const * piCols
24                         , double const * pAlpha, double const * pBeta, double* pData, int const * pLDData);
25 extern void C2F(dlaswp)(int const* n, double* a, int const* ldA, int const* k1, int const* k2, int const* iPiv, int const* incX);
26
27
28 int iLuM(double* pData, int iRows, int iCols, int complexArg, double* pdblLData, double* pdblUData, double* pdblEData )
29 {
30     int ret;
31     int* piPivot;
32     int* piWork = NULL;
33     int iMinRowsCols;
34     double* pdblWork = NULL;
35     iMinRowsCols = Min( iRows, iCols);
36     piPivot = (int*) MALLOC(iMinRowsCols * sizeof(int));
37     if (pdblEData == NULL) /* must allocate more tmp memory when not computing E */
38     {
39         piWork = (int*) MALLOC(iRows * sizeof(int));
40         pdblWork = (double*) MALLOC(iRows * iMinRowsCols * ( complexArg ? sizeof(doublecomplex) : sizeof(double)));
41     }
42     ret = iLu(pData, iRows, iCols, complexArg, pdblLData, pdblUData, pdblEData, piPivot, piWork, pdblWork);
43     if (pdblEData == NULL)
44     {
45         FREE(piWork);
46         FREE(pdblWork);
47     }
48     FREE(piPivot);
49     return ret;
50 }
51 int iLu(double* pData, int iRows, int iCols, int complexArg, double* pdblLData, double* pdblUData, double* pdblEData
52         , int* piPivot, int* piWork, double* pdblWork)
53 {
54     int ret = 0 ;
55     int iInfo;
56     if (complexArg)
57     {
58         C2F(zgetrf)(&iRows, &iCols, pData, &iRows, piPivot, &iInfo);
59     }
60     else
61     {
62         C2F(dgetrf)(&iRows, &iCols, pData, &iRows, piPivot, &iInfo);
63     }
64     if (iInfo < 0)
65     {
66         ret = iInfo;
67     }
68     else
69     {
70         int iMinRowsCols;
71         int i, j;
72         double* pdCurrentL;
73         double* pdCurrentU;
74         double* pdCurrentArg;
75         iMinRowsCols = Min(iRows, iCols);
76         /* fill L matrix with 1.[+0i] on the diagonal */
77         /*
78         /!\ if E was not requested, we will have to swap rows before returning the result, so we do not
79         fill the real L matrix, but pdblWork instead. We will then copy and swap rows from pdblWork to the result matrix.
80               */
81         pdCurrentL = (pdblWork == NULL) ? pdblLData : pdblWork ;
82         pdCurrentArg = pData;
83         for ( j = 0; j != iMinRowsCols; ++j)
84         {
85             for ( i = 0; i != iRows; ++i, ++pdCurrentL, ++pdCurrentArg)
86             {
87                 if (i > j)
88                 {
89                     *pdCurrentL = *pdCurrentArg;
90                     if (complexArg)
91                     {
92                         /* copy imaginary part */
93                         *++pdCurrentL = *++pdCurrentArg;
94                     }
95                 }
96                 else
97                 {
98                     *pdCurrentL = (i == j) ? 1. : 0.;
99                     if (complexArg)
100                     {
101                         /* set imaginary part to 0. and advance pdCurrentArg */
102                         *++pdCurrentL = 0.;
103                         ++pdCurrentArg;
104                     }
105                 }
106             }
107         }
108
109
110         /* init U */
111         pdCurrentU = pdblUData;
112         pdCurrentArg = pData;
113         for ( j = 0; j != iCols; ++j)
114         {
115             pdCurrentArg = pData + j * iRows * (complexArg ? 2 : 1);
116             for ( i = 0; i != iMinRowsCols; ++i, ++pdCurrentU, ++pdCurrentArg)
117             {
118                 if ( i <= j)
119                 {
120                     *pdCurrentU = *pdCurrentArg;
121                     if (complexArg)
122                     {
123                         /* copy imaginary part */
124                         *++pdCurrentU = *++pdCurrentArg;
125                     }
126                 }
127                 else
128                 {
129                     *pdCurrentU = 0.;
130                     if (complexArg)
131                     {
132                         *++pdCurrentU = 0.;
133                         ++pdCurrentArg;
134                     }
135                 }
136             }
137         }
138         if ( pdblEData == NULL)
139         {
140             for (i = 0; i != iRows; ++i)
141             {
142                 piWork[i] = i;
143             }
144             for (i = 0; i != iMinRowsCols; ++i)
145             {
146                 int ip;
147                 ip = piPivot[i] - 1; /* -1 because piPivot contains Fortran index*/
148                 if ( ip != i)
149                 {
150                     int swapTmp;
151                     swapTmp = piWork[i];
152                     piWork[i] = piWork[ip];
153                     piWork[ip] = swapTmp;
154                 }
155             }
156             for (i = 0; i != iRows; ++i)
157             {
158                 for (j = 0; j != iMinRowsCols; ++j)
159                 {
160                     if (complexArg)
161                     {
162                         *(((doublecomplex*)pdblLData) + piWork[i] + j * iRows) = *(((doublecomplex*)pdblWork) + i + j * iRows);
163                     }
164                     else
165                     {
166                         *(pdblLData + piWork[i] + j * iRows) = *(pdblWork + i + j * iRows);
167                     }
168                 }
169             }
170         }/* end if E was not requested */
171         else
172         {
173             double const dblZero = 0., dblOne = 1.;
174             int const iOne = 1;
175             C2F(dlaset)("F", &iRows, &iRows, &dblZero, &dblOne, pdblEData, &iRows);
176             C2F(dlaswp)(&iRows, pdblEData, &iRows, &iOne, &iMinRowsCols, piPivot, &iOne);
177         }
178     }/* end if iInfo reported an error in [z|d]getrf */
179     return ret;
180 }