Least Action

Nontrivializing triviality..and vice versa.

Archive for January 2011

Evaluation of Complex Error Functions erf(z) using GSL/C

with 5 comments

For a small QFT calculation, I needed to numerically evaluate the imaginary error function erfi(x) = erf(i x). But it turns out that GSL (and most other numerical recipe code I could find) can only deal with erf(x), where x is real. Here’s a poor man’s implementation of erf(z) through a standard Taylor expansion,

erf(z) = \frac{2}{\sqrt{\pi}}\sum_{n = 0}^{\infty}\frac{(-1)^n z^{2n+1}}{n!(2n+1)}

The catch here is to deal with propagation of errors in the complex Taylor series, and to also somehow benchmark the results. Well, I haven’t been able to think about this yet, but I was able to confirm that for real arguments, my Taylor series code is about as good as the GSL error function gsl_sf_erf(x).

I am working on a CUDA implementation of this now, because in my project, I need to perform a numerical integration over the error function, which is quite intensive even for Mathematica. For now, I’m just sharing the serial implementation. (Hint: if you can write wrapper for each of the gsl functions inside my Taylor series calculating method, you can embed them within a __global__ kernel call through a struct. Alternatively — and less fun — you can just parallelize the for loop.)

/* Function to compute erf(z) using a Taylor series expansion
/* Author: Vivek Saxena
/* Last updated: January 28, 2011 21:11 hrs */
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
#include <stdio.h>

#define PI 3.1415926543
double cz = 2/sqrt(PI);
const int TERMS = 10; // no of terms to use in the Taylor series

gsl_complex erfTaylor( gsl_complex z, int trunc ){
  gsl_complex res = gsl_complex_rect(0,0),
  num  = gsl_complex_rect(0,0),
  den  = gsl_complex_rect(1,0),
  snum = gsl_complex_rect(1,0),
  temp = gsl_complex_rect(0,0);
  for(int i = 0; i < trunc; i++){
     snum = gsl_complex_rect( cz * gsl_sf_pow_int(-1, i), 0 );
     num  = gsl_complex_mul(snum, gsl_complex_pow_real(z, 2*i+1));
     den  = gsl_complex_rect((2*i + 1)*gsl_sf_fact(i),0);
     temp = gsl_complex_div(num, den);
     res  = gsl_complex_add(res, temp);
  }

 return res;
}

int main ( void ){
   printf( "Real error function\n\n");
   for ( float i = 0; i <= 1; i += 0.01 ){
     float gslerror = gsl_sf_erf(i);
     float taylor   = GSL_REAL( erfTaylor( gsl_complex_rect(i, 0), 10 ) );
     printf("erf(%f): gsl = %f, taylor = %f, mag error = %f\n", i, gslerror, taylor, abs(gslerror-taylor));
   }

   gsl_complex t, arg;
   printf( "\n\nImaginary error function\n\n");
   for (float i = 0; i <= 1; i += 0.01 ){
      /* this would be your generic argument z in erf(z).
       * so if z = x + iy, then
       * arg = gsl_complex_rect(x, y);
      */
      arg = gsl_complex_rect(0, i);
      t   = erfTaylor( arg, TERMS );
      printf("erf(%f + i %f) = %f + i %f\n", GSL_REAL(arg), GSL_IMAG(arg), GSL_REAL(t), GSL_IMAG(t));
   }
   return 0;
}

Written by Vivek

January 28, 2011 at 21:12

CUDA on Ubuntu Maverick Meerkat 10.10

with one comment

(To be edited)

If you’re getting errors like

/usr/bin/ld: cannot find -lGL
/usr/bin/ld: cannot find -lGLU
/usr/bin/ld: cannot find -lX11
/usr/bin/ld: cannot find -lXi
/usr/bin/ld: cannot find -lXmu
/usr/bin/ld: cannot find -lglut

then you need to do the following

sudo apt-get install libxi libxi-dev

Running make again will result in the following errors:

/usr/bin/ld: cannot find -lGL
/usr/bin/ld: cannot find -lGLU
/usr/bin/ld: cannot find -lXmu
/usr/bin/ld: cannot find -lglut

Now, let’s execute:

sudo apt-get install freeglut3 freeglut3-dev

This brings down the errors to

/usr/bin/ld: cannot find -lXmu

So, we just have to do one more apt-get:

sudo apt-get install libxmu6 libxmu-dev

—-
If you get errors suggesting that your libcudart.so.1 is missing, it means your LD_CONFIG_PATH isn’t set right. To set it permanently, use

sudo ldconfig -v /usr/local/cuda/lib64/

on a 64-bit system and

sudo ldconfig -v /usr/local/cuda/lib/

on a 32 bit system.

Written by Vivek

January 25, 2011 at 22:49

CUDA on a Dell XPS 15 in Windows 7 64-bit

leave a comment »

I just figured out how to get NVIDIA CUDA to work on my laptop. You need to replace the generic Dell driver with this one:

266.58_notebook_winvista_win7_64bit_international_whql.exe

You need to do a custom clean installation and make sure the PhysX box is checked.

I had to do some CUDA programming on the Windows partition, so now I have to figure out how to configure all my IDEs to work with CUDA. I will try and post detailed configuration info for Netbeans at least.

Written by Vivek

January 22, 2011 at 11:07