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