diff --git a/configure.ac b/configure.ac index 291e9594a..5eb4d3b8c 100644 --- a/configure.ac +++ b/configure.ac @@ -20,6 +20,7 @@ AC_ENABLE_SHARED ############################################################################## # Check for mpiCC immediately after getting C++ compiler... +AC_PROG_CC AC_PROG_CXX AC_LANG([C++]) @@ -569,6 +570,7 @@ AC_CONFIG_FILES([ Makefile meep-pkgconfig src/Makefile + src/support/Makefile tests/Makefile ]) diff --git a/src/Makefile.am b/src/Makefile.am index b5be070d1..ad571b2d2 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -19,6 +19,9 @@ sources.cpp step.cpp step_db.cpp stress.cpp structure.cpp structure_dump.cpp \ susceptibility.cpp time.cpp update_eh.cpp mpb.cpp update_pols.cpp \ vec.cpp step_generic.cpp meepgeom.cpp GDSIIgeom.cpp $(HDRS) $(BUILT_SOURCES) +SUBDIRS = support +libmeep_la_LIBADD = support/libsupport.la + libmeep_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@ PRELUDE = "/* This file was automatically generated --- DO NOT EDIT */" diff --git a/src/initialize.cpp b/src/initialize.cpp index bb1b41d49..276f259e1 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -36,10 +36,10 @@ namespace meep { #define J BesselJ double J(int m, double kr) { -#if defined(HAVE_LIBGSL) - return gsl_sf_bessel_Jn(m, kr); -#elif defined(HAVE_JN) +#if defined(HAVE_JN) return jn(m, kr); // POSIX/BSD jn function +#elif defined(HAVE_LIBGSL) + return gsl_sf_bessel_Jn(m, kr); #else abort("not compiled with GSL, required for Bessel functions"); return 0; diff --git a/src/near2far.cpp b/src/near2far.cpp index 7c89ca019..a6b570a7f 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -157,9 +157,13 @@ void green3d(std::complex *EH, const vec &x, double freq, double eps, do abort("unrecognized source type"); } -#ifdef HAVE_LIBGSL -#include // hankel function J + iY +#if defined(HAVE_JN) +static std::complex hankel(int n, double x) { + return std::complex(jn(n, x), yn(n, x)); +} +#elif defined(HAVE_LIBGSL) +#include static std::complex hankel(int n, double x) { return std::complex(gsl_sf_bessel_Jn(n, x), gsl_sf_bessel_Yn(n, x)); } diff --git a/src/random.cpp b/src/random.cpp index d13dc1eff..757fdf792 100644 --- a/src/random.cpp +++ b/src/random.cpp @@ -18,11 +18,7 @@ #include "meep.hpp" #include "config.h" -#ifdef HAVE_LIBGSL -#include -#else -#include -#endif +#include "support/meep_mt.h" #include using namespace std; @@ -30,16 +26,9 @@ using namespace std; namespace meep { static bool rand_inited = false; -#ifdef HAVE_LIBGSL -static gsl_rng *rng = NULL; -#endif static void init_rand(void) { if (!rand_inited) { -#ifdef HAVE_LIBGSL - if (rng) gsl_rng_free(rng); - rng = gsl_rng_alloc(gsl_rng_mt19937); -#endif rand_inited = true; // no infinite loop since rand_inited == true set_random_seed(time(NULL) * (1 + my_global_rank())); } @@ -47,43 +36,27 @@ static void init_rand(void) { void set_random_seed(unsigned long seed) { init_rand(); - seed = ((unsigned long)broadcast(0, (int)seed)) + my_rank(); -#ifdef HAVE_LIBGSL - gsl_rng_set(rng, seed); -#else - srand(seed); -#endif + meep_mt_init_genrand(seed); } int random_int(int a, int b) { init_rand(); -#ifdef HAVE_LIBGSL - return ((int)gsl_rng_uniform_int(rng, b - a + 1)) + a; -#else - return a + rand() % (b - a + 1); -#endif + return a + meep_mt_genrand_int32() % (b - a + 1); } double uniform_random(double a, double b) { init_rand(); -#ifdef HAVE_LIBGSL - return a + gsl_rng_uniform(rng) * (b - a); -#else - return a + rand() * (b - a) / RAND_MAX; -#endif + return a + meep_mt_genrand_res53() * (b - a); } double gaussian_random(double mean, double stddev) { init_rand(); -#ifdef HAVE_LIBGSL - return mean + gsl_ran_gaussian(rng, stddev); -#else // Box-Muller algorithm to generate Gaussian from uniform // see Knuth vol II algorithm P, sec. 3.4.1 double v1, v2, s; do { - v1 = uniform_random(-1, 1); - v2 = uniform_random(-1, 1); + v1 = 2*meep_mt_genrand_res53() - 1; + v2 = 2*meep_mt_genrand_res53() - 1; s = v1 * v1 + v2 * v2; } while (s >= 1.0); if (s == 0) { @@ -91,7 +64,6 @@ double gaussian_random(double mean, double stddev) { } else { return mean + v1 * sqrt(-2 * log(s) / s) * stddev; } -#endif } } // namespace meep diff --git a/src/support/Makefile.am b/src/support/Makefile.am new file mode 100644 index 000000000..f65ba7ddb --- /dev/null +++ b/src/support/Makefile.am @@ -0,0 +1,3 @@ +noinst_LTLIBRARIES = libsupport.la + +libsupport_la_SOURCES = mt19937ar.c meep_mt.h diff --git a/src/support/meep_mt.h b/src/support/meep_mt.h new file mode 100644 index 000000000..73f280819 --- /dev/null +++ b/src/support/meep_mt.h @@ -0,0 +1,21 @@ +#ifndef MEEP_MT_H +#define MEEP_MT_H 1 + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +void meep_mt_init_genrand(unsigned long s); +void meep_mt_init_by_array(unsigned long init_key[], int key_length); +unsigned long meep_mt_genrand_int32(void); +long meep_mt_genrand_int31(void); +double meep_mt_genrand_real1(void); +double meep_mt_genrand_real2(void); +double meep_mt_genrand_real3(void); +double meep_mt_genrand_res53(void); + +#ifdef __cplusplus +} +#endif /* __cplusplus */ + +#endif /* MEEP_MT_H */ diff --git a/src/support/mt19937ar.c b/src/support/mt19937ar.c new file mode 100644 index 000000000..b1075e379 --- /dev/null +++ b/src/support/mt19937ar.c @@ -0,0 +1,196 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +/* modified by SGJ for Meep: added meep_mt.h header and added meep_mt prefixes + to functions to avoid namespace pollution */ + +#include +#include "meep_mt.h" + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ + +static unsigned long mt[N]; /* the array for the state vector */ +static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ + +/* initializes mt[N] with a seed */ +void meep_mt_init_genrand(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +void meep_mt_init_by_array(unsigned long init_key[], int key_length) +{ + int i, j, k; + meep_mt_init_genrand(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long meep_mt_genrand_int32(void) +{ + unsigned long y; + static unsigned long mag01[2]={0x0UL, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if init_genrand() has not been called, */ + meep_mt_init_genrand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* generates a random number on [0,0x7fffffff]-interval */ +long meep_mt_genrand_int31(void) +{ + return (long)(meep_mt_genrand_int32()>>1); +} + +/* generates a random number on [0,1]-real-interval */ +double meep_mt_genrand_real1(void) +{ + return meep_mt_genrand_int32()*(1.0/4294967295.0); + /* divided by 2^32-1 */ +} + +/* generates a random number on [0,1)-real-interval */ +double meep_mt_genrand_real2(void) +{ + return meep_mt_genrand_int32()*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on (0,1)-real-interval */ +double meep_mt_genrand_real3(void) +{ + return (((double)meep_mt_genrand_int32()) + 0.5)*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on [0,1) with 53-bit resolution*/ +double meep_mt_genrand_res53(void) +{ + unsigned long a=meep_mt_genrand_int32()>>5, b=meep_mt_genrand_int32()>>6; + return(a*67108864.0+b)*(1.0/9007199254740992.0); +} +/* These real versions are due to Isaku Wada, 2002/01/09 added */ + +#if 0 +int main(void) +{ + int i; + unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4; + init_by_array(init, length); + printf("1000 outputs of genrand_int32()\n"); + for (i=0; i<1000; i++) { + printf("%10lu ", genrand_int32()); + if (i%5==4) printf("\n"); + } + printf("\n1000 outputs of genrand_real2()\n"); + for (i=0; i<1000; i++) { + printf("%10.8f ", genrand_real2()); + if (i%5==4) printf("\n"); + } + return 0; +} +#endif