From octave-sources-request at bevo dot che dot wisc dot edu Fri Jun 18 15:22:50 1999 Subject: Octave and GSL-RNG interface From: Thomas Walter To: octave-sources at bevo dot che dot wisc dot edu Date: Fri, 18 Jun 1999 22:23:24 +0200 Hello, I wrote an interface between octave and the GSL random number generator part. The source at the end contains two callable functions. The first is a setup and the second calls the RNG from GSL. Note seeding at runtime is currently not supported. As suggested in the help-octave list I used for timing a call like: [tst0,tsu0,tss0]=cputime;r=rand(1000000,1);[tst1,tsu1,tss1]=cputime;tst1-tst0,tsu1-tsu0,tss1-tss0 and [tst0,tsu0,tss0]=cputime;r=rng_rand(1000000,1);[tst1,tsu1,tss1]=cputime;tst1-tst0,tsu1-tsu0,tss1-tss0 Use this large number or you may not see a time difference depending which RNG you use from GSL. Notes about compiling are at the head of the source below. I used 'gsl-0.4.1', 'octave-2.0.14' and 'egcs-1.1b' on Linux. Hints are welcome. Nice weekend Thomas ======== gsl_rng.cc start ================= // Hi emacs! My mode is -*- C++ -*- // // This implements an interface to the GSL random number generators // // Copyright (C) 1999 Thomas Walter // $Id: gsl_rng.cc,v 1.6 1999/06/18 18:06:17 twalter Exp $ // // // ======================= Compile instruction ================ // // If you want to compile the source to the oct file use // something similar to the following in a 'Makefile' // It assume the use of 'g++' currently from 'egcs-1.1b' // // RM = rm -f // LN_S = ln -s // // DBUG_SIMPLE = -W -Wall // DBUG = $(DBUG_SIMPLE) #$(DBUG_FULL) // // #MKOCTCMD = ./mkoctfile // MKOCTCMD = /usr/local/bin/mkoctfile // # request for verbose and strip // MKOCTFLAGS = -v -s // // gsl_rng.oct : gsl_rng.cc // env SH_LDFLAGS='-shared -Wl,-rpath,/usr/local/lib/gsl -L/usr/local/lib/gsl' CXX="g++ -pipe" CPPFLAGS="$(DBUG)" $(MKOCTCMD) $(MKOCTFLAGS) -I/usr/local/include/gsl -lgslrng -lgslerr gsl_rng.cc // # $(LN_S) gsl_rng.oct rng_setup_by_env.oct // # $(LN_S) gsl_rng.oct rng_rand.oct // $(RM) gsl_rng.o // // You may have to adapt the 'SH_LDFLAGS' and 'CXX' and 'MKOCTCMD'. // Enable the 'LN_S' command only on the first time. // // ======================= Compile instruction ================ // #include //#include /* for INT_MAX */ #include // I want 'feval()' #ifdef __cplusplus extern "C" { #endif /* __cplusplus */ #include #include #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ static const char *name_of_current_rng = NULL; static gsl_rng *rng_in_use = NULL; DEFUN_DLD (rng_setup_by_env, /* nargin not used */, /* nargout not used */, "rng_setup_by_env (): Enable the default RNG of GSL. The name is either from the environment variable 'GSL_RNG_TYPE' or the built in type of GLS wich is 'gsl_rng_mt19937'.") { const gsl_rng_type *used_rng_type; // Use either the RNG defined by environment variable 'GSL_RNG_TYPE' // or the default 'mt19937'. // Use either the SEED defined by environment variable 'GSL_RNG_SEED' // or the default seed '0'. // This signals the RNG to use its internal defaults. used_rng_type = gsl_rng_env_setup (); name_of_current_rng = used_rng_type->name; if (rng_in_use != NULL) { gsl_rng_free (rng_in_use); rng_in_use = NULL; } rng_in_use = gsl_rng_alloc (used_rng_type); return (name_of_current_rng); } static bool do_rand (Matrix& cv, const int nx, const int ny) { octave_value retval; if ((nx < 1) || (nx > (INT_MAX - 1))) { // cout << "Matrix/Vector too large" << endl; return (false); } // If I don't have a RNG ... if (rng_in_use == NULL) { feval ("rng_setup_by_env", 1); // ... I get one. } cv.resize (nx, ny); #if 1 // Test direct access to the data double *data = cv.fortran_vec (); unsigned int i; const unsigned int nxny = nx * ny; for (i = 0; i < nxny; i++) { data[i] = gsl_rng_uniform (rng_in_use); } #else // Test access as Matrix int i; for (i = 0; i < nx; i++) { int j; for (j = 0; j < ny; j++) { cv (i, j) = gsl_rng_uniform (rng_in_use); } } #endif return (true); } DEFUN_DLD (rng_rand, args , /* nargout */, "rng_rand (n): Return n random values with uniform distribution using default RNG.\nValid values for arument n: 1 <= n < INT_MAX\nThe returned value(s) are all in the intervall '[0,1['\n\n\ rng_rand -- generate a random value from a uniform distribution rng_rand (N) -- generate N x N matrix rng_rand (size (A)) -- generate matrix the size of A rng_rand (N, M) -- generate N x M matrix ===== The folowing are not implemented ========== rng_rand (\"seed\") -- get current seed rng_rand (\"seed\", N) -- set seed") { const int nargin = args.length (); octave_value retval; // Sanity checks if ((nargin < 0) || (nargin > 2)) { print_usage ("rng_rand"); return (retval); } // Init with default. This is used if there are no arguments int how_many = 1; int how_many_cols = 1; if (nargin > 0) { if (args (0).is_matrix_type ()) // Required to support ' a = rand (size(M)) ' { // This part is ugly, but seems to work const Matrix m = args (0).matrix_value (); const int mr = m.rows (); const int mc = m.columns (); how_many = m (0, 0); // row and column vectors first element how_many_cols = (mc == 1) ? m (1, 0) : m (0, 1); //how_many_cols = (mc == 1) ? mr : mc; #if DEBUG_ME cout << "m = " << m < 1) && (mc > 1)) { #if DEBUG_ME cout << "The values I got are : mr = " << mr << " mc = " << mc << endl; cout << "The values I got are : how_many = " << how_many << " how_many_cols = " << how_many_cols << endl; cout << "The values I got are : m(0,0) = " << m(0,0) << " m(1,1) = " << m(1,1) << endl; cout << "The values I got are : m(1,0) = " << m(1,0) << " m(0,1) = " << m(0,1) << endl; print_usage ("rng_rand"); return (retval); #else // !DEBUG_ME cout << "please use 'rng_rand (size(m))'" << endl; how_many = mr; how_many_cols = mc; #endif // !DEBUG_ME } // ((mr > 1) && (mc > 1)) } else // !is_matrix_type { how_many = (nargin < 1) ? 1 : rint (abs (args (0).double_value ())); // octave-2.0.14 how_many_cols = (nargin < 2) ? how_many : rint (abs (args (1).double_value ())); // octave-2.0.14 // const int how_many = arg.int_value (); // octave-2.1.14 } // !is_matrix_type } // (nargin > 0) // // Now the real work // Matrix cv; bool success = do_rand (cv, how_many, how_many_cols); if (!success) { print_usage ("rng_rand"); return (retval); } // Why have I to cast it in 'octave-2.0.14'? // Without the cast the compiler reports an error ????? return (octave_value (cv)); } ======== gsl_rng.cc end ================= -- Platzangst: Der dauerhafte Zustand eines Luftballons. ---------------------------------------------- Dipl. Phys. Thomas Walter Inst. f. Physiklische Chemie II Egerlandstr. 3 Tel.: ++9131-85 27326 / 27330 91058 Erlangen, Germany email: walter at pctc dot chemie dot uni-erlangen dot de