1 /***************************************************************************
2                           gsl_matrix.cpp  -  GSL GDL library function
3                              -------------------
4     begin                : Dec 9 2011
5     copyright            : (C) 2011 by Alain Coulais
6     email                : alaingdl@users.sourceforge.net
7  ***************************************************************************/
8 
9 /***************************************************************************
10  *                                                                         *
11  *   This program is free software; you can redistribute it and/or modify  *
12  *   it under the terms of the GNU General Public License as published by  *
13  *   the Free Software Foundation; either version 2 of the License, or     *
14  *   (at your option) any later version.                                   *
15  *                                                                         *
16  ***************************************************************************/
17 
18 #include "includefirst.hpp"
19 
20 #include <map>
21 #include <cmath>
22 
23 #include "datatypes.hpp"
24 #include "envt.hpp"
25 #include "basic_fun.hpp"
26 #include "gsl_fun.hpp"
27 #include "dinterpreter.hpp"
28 
29 #include <gsl/gsl_sys.h>
30 #include <gsl/gsl_linalg.h>
31 
32 // constant
33 #include <gsl/gsl_math.h>
34 
35 //#define LOG10E 0.434294
36 
37 namespace lib {
38 
39   using namespace std;
40 
41   const int szdbl=sizeof(DDouble);
42   const int szflt=sizeof(DFloat);
43   const int szlng=sizeof(DLong);
44   const int szlng64=sizeof(DLong64);
45 
ludc_pro(EnvT * e)46   void ludc_pro( EnvT* e)
47   {
48     //  cout << szdbl << " " <<szflt << " " << szlng << " " szlng64 << endl;
49 
50     SizeT nParam=e->NParam(1);
51     //     if( nParam == 0)
52     //       e->Throw( "Incorrect number of arguments.");
53 
54     BaseGDL* p0 = e->GetParDefined( 0);
55 
56     SizeT nEl = p0->N_Elements();
57     if( nEl == 0)
58       e->Throw( "Variable is undefined: " + e->GetParString(0));
59 
60     if (p0->Rank() > 2)
61       e->Throw( "Input must be a square matrix:" + e->GetParString(0));
62 
63     if (p0->Rank() > 1) {
64       if (p0->Dim(0) != p0->Dim(1))
65         e->Throw( "Input must be a square matrix:" + e->GetParString(0));
66     }
67 
68     // status
69     // check here, if not done, res would be pending in case of SetPar() throws
70     // SetPar() only throws in AssureGlobalPar()
71     if (nParam == 2) e->AssureGlobalPar( 1);
72 
73     // only one element matrix
74 
75     if (p0->Type() == GDL_COMPLEXDBL || p0->Type() == GDL_COMPLEX){
76       e->Throw( "Input type cannot be COMPLEX, please use LA_LUDC (not ready)");
77     }
78 
79     //DDoubleGDL* p0D = static_cast<DDoubleGDL*>( p0);
80     DDoubleGDL *p0D = e->GetParAs<DDoubleGDL>(0);
81 
82     gsl_matrix *mat = gsl_matrix_alloc(p0->Dim(0), p0->Dim(0));
83     GDLGuard<gsl_matrix> g1(mat,gsl_matrix_free);
84 
85     memcpy(mat->data, &(*p0D)[0], nEl*szdbl);
86 
87     gsl_permutation * p = gsl_permutation_alloc (p0->Dim(0));
88     GDLGuard<gsl_permutation> g2( p, gsl_permutation_free);
89     int s;
90     gsl_linalg_LU_decomp (mat, p, &s);
91 
92     int debug=0;
93     if (debug) {
94       cout << "permutation order: " << s << endl;
95       cout << "permutation vector:"<< endl;
96       gsl_permutation_fprintf (stdout, p, " %u");
97       cout << endl;
98     }
99 
100     // copying over p0 the updated matrix
101     SizeT dims[2] = {p0->Dim(0), p0->Dim(0)};
102     dimension dim0(dims, (SizeT) 2);
103 
104     BaseGDL** p0Do = &e->GetPar( 0);
105     GDLDelete((*p0Do));
106     *p0Do = new DDoubleGDL(dim0, BaseGDL::NOZERO);
107     memcpy(&(*(DDoubleGDL*) *p0Do)[0], mat->data,
108 	   p0->Dim(0)*p0->Dim(0)*szdbl);
109 
110     int double_flag=0;
111     if (p0->Type() == GDL_DOUBLE) double_flag=1;
112     static int doubleIx=e->KeywordIx("DOUBLE");
113     if (e->KeywordSet(doubleIx)) double_flag=1;
114 
115     // this code will always return GDL_DOUBLE because I don't know how to do :(
116     // AC 13-Feb-2012 : this is not working and I don't know how to do :(
117     // if (double_flag == 0)
118     // { p0->Convert2(GDL_FLOAT, BaseGDL::CONVERT); }
119 
120     // copying over p1 the permutation vector
121     SizeT n = p0->Dim(0);
122     dimension dim1(&n, (SizeT) 1);
123     BaseGDL** p1D = &e->GetPar( 1);
124     GDLDelete((*p1D));
125     if (sizeof(size_t) == szlng) {
126       *p1D = new DLongGDL(dim1, BaseGDL::NOZERO);
127       memcpy(&(*(DLongGDL*) *p1D)[0], p->data,
128 	p0->Dim(0)*szlng);
129     } else {
130       *p1D = new DLong64GDL(dim1, BaseGDL::NOZERO);
131       memcpy(&(*(DLong64GDL*) *p1D)[0], p->data,
132 	p0->Dim(0)*szlng64);
133     }
134 
135 //     gsl_matrix_free(mat);
136 //     gsl_permutation_free(p);
137 
138   }
139 
lusol_fun(EnvT * e)140   BaseGDL* lusol_fun( EnvT* e)
141   {
142     SizeT nParam=e->NParam(1);
143 //    int s;
144 
145     //     if( nParam == 0)
146     //       e->Throw( "Incorrect number of arguments.");
147 
148     // managing first input: Square Matrix
149 
150     BaseGDL* p0 = e->GetParDefined( 0);
151 
152     SizeT nEl = p0->N_Elements();
153     if( nEl == 0)
154       e->Throw( "Variable is undefined: " + e->GetParString(0));
155 
156     if (p0->Rank() > 2)
157       e->Throw( "Input must be a square matrix:" + e->GetParString(0));
158 
159     if (p0->Rank() > 1) {
160       if (p0->Dim(0) != p0->Dim(1))
161         e->Throw( "Input must be a square matrix:" + e->GetParString(0));
162     }
163 
164     // status
165     // check here, if not done, res would be pending in case of SetPar() throws
166     // SetPar() only throws in AssureGlobalPar()
167     if (nParam == 2) e->AssureGlobalPar( 1);
168 
169 
170     // managing seconf input: Vector
171     BaseGDL* p1 = e->GetParDefined(1);
172 
173     SizeT nEl1 = p1->N_Elements();
174     if( nEl1 == 0)
175       e->Throw( "Variable is undefined: " + e->GetParString(1));
176 
177     if (p1->Rank() > 2)
178       e->Throw( "Input must be a Vector:" + e->GetParString(1));
179 
180 
181     // managing third input: Vector
182     BaseGDL* p2 = e->GetParDefined(2);
183 
184     SizeT nEl2 = p2->N_Elements();
185     if( nEl2 == 0)
186       e->Throw( "Variable is undefined: " + e->GetParString(2));
187 
188     if (p2->Rank() > 2)
189       e->Throw( "Input must be a Vector:" + e->GetParString(2));
190 
191     if (p0->Type() == GDL_COMPLEXDBL || p0->Type() == GDL_COMPLEX){
192       e->Throw( "Input type cannot be COMPLEX, please use LA_LUDC (not ready)");
193     }
194 
195     DDoubleGDL *p0D = e->GetParAs<DDoubleGDL>(0);
196     gsl_matrix *mat = gsl_matrix_alloc(p0->Dim(0), p0->Dim(0));
197     GDLGuard<gsl_matrix> g1(mat,gsl_matrix_free);
198     memcpy(mat->data, &(*p0D)[0], nEl*szdbl);
199 
200     gsl_permutation *p = gsl_permutation_alloc (nEl1);
201     GDLGuard<gsl_permutation> g2(p,gsl_permutation_free);
202     if (sizeof(size_t) == szlng) {
203       DLongGDL* p1L =e->GetParAs<DLongGDL>(1);
204       memcpy(p->data, &(*p1L)[0], nEl1*szlng);
205     } else {
206       DLong64GDL* p1L =e->GetParAs<DLong64GDL>(1);
207       memcpy(p->data, &(*p1L)[0], nEl1*szlng64);
208     }
209 
210     DDoubleGDL *p2D = e->GetParAs<DDoubleGDL>(2);
211     gsl_vector *b = gsl_vector_alloc(nEl2);
212     GDLGuard<gsl_vector> g3(b,gsl_vector_free); // b was NOT freed before
213     memcpy(b->data, &(*p2D)[0], nEl1*szdbl);
214 
215     gsl_vector *x = gsl_vector_alloc(nEl2);
216     GDLGuard<gsl_vector> g4(x,gsl_vector_free);
217 
218     // computation by GSL
219     gsl_linalg_LU_solve (mat, p, b, x);
220 
221     int debug=0;
222     if (debug) {
223 //      cout << "permutation order: " << s << endl;
224       cout << "permutation vector:";
225       gsl_permutation_fprintf (stdout, p, " %u");
226       cout << endl;
227       cout << "input vector:";
228       gsl_vector_fprintf (stdout, b, " %g");
229       cout << endl;
230       cout << "result vector:";
231       gsl_vector_fprintf (stdout, x, " %g");
232       cout << endl;
233     }
234 
235     DDoubleGDL* res = new DDoubleGDL( p2->Dim(), BaseGDL::NOZERO);
236     memcpy(&(*res)[0], x->data, nEl1*szdbl);
237 
238 //     gsl_matrix_free(mat);Parameter
239 //     gsl_vector_free(x);
240 //     gsl_permutation_free(p);
241 //     b ???
242 
243     int double_flag=0;
244     if (p0->Type() == GDL_DOUBLE || p2->Type() == GDL_DOUBLE) double_flag=1;
245     static int doubleIx=e->KeywordIx("DOUBLE");
246     if (e->KeywordSet(doubleIx)) double_flag=1;
247 
248     if (double_flag)
249       {	return res; }
250     else
251       { return res->Convert2(GDL_FLOAT, BaseGDL::CONVERT); }
252   }
253 
determ_fun(EnvT * e)254   BaseGDL* determ_fun( EnvT* e) {
255 
256     // managing first input: Square Matrix
257 
258     BaseGDL* p0 = e->GetParDefined( 0);
259 
260     SizeT nEl = p0->N_Elements();
261     if( nEl == 0)
262       e->Throw( "Variable is undefined: " + e->GetParString(0));
263 
264     if (p0->Rank() > 2)
265       e->Throw( "Input must be a square matrix:" + e->GetParString(0));
266 
267     if (p0->Rank() > 1) {
268       if (p0->Dim(0) != p0->Dim(1))
269         e->Throw( "Input must be a square matrix:" + e->GetParString(0));
270     }
271 
272     if (p0->Type() == GDL_COMPLEXDBL || p0->Type() == GDL_COMPLEX){
273       e->Throw( "Input type cannot be COMPLEX, please use LA_DETERM (not ready)");
274     }
275 
276     DDoubleGDL *p0D = e->GetParAs<DDoubleGDL>(0);
277     gsl_matrix *mat = gsl_matrix_alloc(p0->Dim(0), p0->Dim(0));
278     GDLGuard<gsl_matrix> g1(mat,gsl_matrix_free);
279     memcpy(mat->data, &(*p0D)[0], nEl*szdbl);
280 
281     gsl_permutation *p = gsl_permutation_alloc(p0->Dim(0));
282     GDLGuard<gsl_permutation> g2(p,gsl_permutation_free);
283 
284     int sign;
285     double determ=0.0;
286 
287     // computation by GSL
288     gsl_linalg_LU_decomp (mat, p, &sign);
289     determ=gsl_linalg_LU_det (mat, sign);
290 
291     int debug=0;
292     if (debug) {
293       cout << "Determ : " << determ << endl;
294     }
295 
296 //     gsl_matrix_free(mat);
297 //     gsl_permutation_free(p);
298 
299     DDoubleGDL* res = new DDoubleGDL(1, BaseGDL::NOZERO);
300     (*res)[0]=determ;
301 
302     int double_flag=0;
303     if (p0->Type() == GDL_DOUBLE) double_flag=1;
304     static int doubleIx=e->KeywordIx("DOUBLE");
305     if (e->KeywordSet(doubleIx)) double_flag=1;
306 
307     if (double_flag)
308       {	return res; }
309     else
310       { return res->Convert2(GDL_FLOAT, BaseGDL::CONVERT); }
311   }
312 
TDMAsolver8(SizeT M,double a[],double b[],double c[],double d[],double x[])313   int TDMAsolver8(SizeT M, double a[], double b[], double c[], double d[], double x[])
314   {
315     /* Tri Diagonal Matrix Algorithm(a.k.a Thomas algorithm) solver
316      TDMA solver, a b c d can be NumPy array type or Python list type.
317      refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
318      */
319     if (b[0] == 0) return 1;
320     double w;
321     DLong i; //not SizeT as unsigned loops fail miserably with decrementing (--i)
322 
323     for (i = 1; i < M; ++i) {
324       if (b[i-1] == 0) return 1;
325       w = a[i] / b[i - 1];
326       b[i] -= w * c[i - 1];
327       d[i] -= w * d[i - 1];
328     }
329     x[M - 1] = d[M - 1] / b[M - 1];
330 
331     for (i = M - 2; i >= 0; --i) {
332       x[i] = (d[i] - c[i] * x[i + 1]) / b[i];
333     }
334     return 0;
335   }
336 
TDMAsolver4(SizeT M,double a[],double b[],double c[],double d[],float x[])337   int TDMAsolver4(SizeT M, double a[], double b[], double c[], double d[], float x[])
338   {
339     /* Tri Diagonal Matrix Algorithm(a.k.a Thomas algorithm) solver
340      TDMA solver, a b c d can be NumPy array type or Python list type.
341      refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
342      */
343     double w;
344     DLong i; //not SizeT as unsigned loops fail miserably with decrementing (--i)
345 
346     for (i = 1; i < M; ++i) {
347       if (b[i-1] == 0) return 1;
348       w = a[i] / b[i - 1];
349       b[i] -= w * c[i - 1];
350       d[i] -= w * d[i - 1];
351     }
352     x[M - 1] = d[M - 1] / b[M - 1];
353 
354     for (i = M - 2; i >= 0; --i) {
355       x[i] = (d[i] - c[i] * x[i + 1]) / b[i];
356     }
357     return 0;
358   }
trisol_fun(EnvT * e)359   BaseGDL* trisol_fun( EnvT* e) {
360 
361     // managing first input: Square Matrix
362 
363     BaseGDL* p0 = e->GetParDefined(0); // sub-diag elements
364     BaseGDL* p1 = e->GetParDefined(1); // diagonal elements
365     BaseGDL* p2 = e->GetParDefined(2); // sup-diag elements
366     BaseGDL* p3 = e->GetParDefined(3); // right-hand side vector
367 
368     SizeT nEl0 = p0->N_Elements();
369     if( nEl0 == 0) e->Throw( "Variable is undefined: " + e->GetParString(0));
370     SizeT nEl1 = p1->N_Elements();
371     if( nEl1 == 0) e->Throw( "Variable is undefined: " + e->GetParString(1));
372     SizeT nEl2 = p2->N_Elements();
373     if( nEl2 == 0) e->Throw( "Variable is undefined: " + e->GetParString(2));
374     SizeT nEl3 = p3->N_Elements();
375     if( nEl3 == 0) e->Throw( "Variable is undefined: " + e->GetParString(3));
376 
377     //    cout << nEl0 << " " << nEl1 << " " << nEl2 << " " << nEl3 << " " << endl;
378 
379     SizeT nEl = nEl0;
380     if (nEl1 != nEl) e->Throw( "Argument: " + e->GetParString(1)+" does not have correct size");
381     if (nEl2 != nEl) e->Throw( "Argument: " + e->GetParString(2)+" does not have correct size");
382     if (nEl3 != nEl) e->Throw( "Argument: " + e->GetParString(3)+" does not have correct size");
383 
384     int complex_flag=0;
385     if (p0->Type() == GDL_COMPLEXDBL || p0->Type() == GDL_COMPLEX) complex_flag=1;
386     if (p1->Type() == GDL_COMPLEXDBL || p1->Type() == GDL_COMPLEX) complex_flag=1;
387     if (p2->Type() == GDL_COMPLEXDBL || p2->Type() == GDL_COMPLEX) complex_flag=1;
388     if (p3->Type() == GDL_COMPLEXDBL || p3->Type() == GDL_COMPLEX) complex_flag=1;
389     if (complex_flag) {
390       e->Throw( "Input type cannot be COMPLEX, please use LA_TRISOL (not ready)");
391     }
392 
393     // computations are done in Double type, conversion at the end
394 
395     DDoubleGDL *p0D = e->GetParAs<DDoubleGDL>(0); //A
396     DDoubleGDL *p1D= e->GetParAs<DDoubleGDL>(1)->Dup();//B modified in call to TDMASolver
397     Guard<BaseGDL> g1(p1D);
398     DDoubleGDL *p2D= e->GetParAs<DDoubleGDL>(2);//C
399     DDoubleGDL *p3D= e->GetParAs<DDoubleGDL>(3)->Dup();//D idem
400     Guard<BaseGDL> g3(p3D);
401     DLong double_flag=0;
402     static int doubleIx=e->KeywordIx("DOUBLE");
403     e->AssureLongScalarKWIfPresent(doubleIx,double_flag);
404     if (double_flag) {
405 
406     DDoubleGDL* res = new DDoubleGDL(nEl, BaseGDL::NOZERO);
407       int err=TDMAsolver8(nEl,(DDouble*)p0D->DataAddr(),(DDouble*)p1D->DataAddr(),(DDouble*)p2D->DataAddr(),(DDouble*)p3D->DataAddr(), (DDouble*)res->DataAddr());
408       if (err > 0) {
409         GDLDelete(res);
410         e->Throw("TRISOL: Error "+i2s(err)+" in tridag");
411       }
412       return res;
413     } else {
414       DFloatGDL* res = new DFloatGDL(nEl, BaseGDL::NOZERO);
415       int err=TDMAsolver4(nEl,(DDouble*)p0D->DataAddr(),(DDouble*)p1D->DataAddr(),(DDouble*)p2D->DataAddr(),(DDouble*)p3D->DataAddr(), (DFloat*)res->DataAddr());
416       if (err > 0) {
417         GDLDelete(res);
418         e->Throw("TRISOL: Error "+i2s(err)+" in tridag");
419       }
420       return res;
421     }
422   }
423 }
424 
425