The Return of Agent Zlerich

michellethompson.org

2012-01-08T04:05:00.002Z

Two years after our friend Nancy snagged the domain, we finally put up some content at http://michellethompson.org/...

Specializations of general banded matrix-vector multiplication

2011-12-24T05:08:00.002Z

I keep needing mixed real/complex level 2 BLAS operations for real-valued, general banded matrices. It's simple to modify Netlib's dgbmv for mixed real/complex data but the result always feels unsatisfactory. Partially because I'd like to avoid a Fortran dependency in my project. Partially because I'd like to build kernels for the special case where I know my matrix's upper and lower bandwidth (as optimizers can rip there according to fellow student Myoungkyu Lee).

I finally got off (err... on?) my duff and reworked dgbmv. The result is a snippet which can be #included from a C99 source file where the proper defines turn it into a real, complex, or mixed real/complex routine. I've uglified it sufficiently that both GNU and Intel's compiler's inform me at -O3 that they've vectorized the loops handling both A and A**T cases for contiguous vectors.

// #include-time parameters available followed by a sample usage.
// Each of these are #undef-ed at the end of this template file
//
// #define GBMV_STATIC    /*empty*/         /* Use for static function    */
// #define GBMV_FUNCTION  dgbmv             /* Function name              */
// #define GBMV_COMPONENT double            /* Type of matrices           */
// #define GBMV_SCALAR    double            /* Type of coeffs and vectors */
// #define GBMV_KL        const int kl,     /* Matrix lower bandwidth     */
// #define GBMV_KU        const int ku,     /* Matrix upper bandwidth     */

#ifndef __GBMV_INTERNAL_ONCE
#define __GBMV_INTERNAL_ONCE
static inline int ndxmin(int a, int b) { return a < b ? a : b; }
static inline int ndxmax(int a, int b) { return a > b ? a : b; }
#endif /* __GBMV_INTERNAL_ONCE */

GBMV_STATIC int GBMV_FUNCTION(
    char trans, const int m, const int n, GBMV_KL GBMV_KU
    const GBMV_SCALAR alpha, const GBMV_COMPONENT *       restrict a, int lda,
                             const    GBMV_SCALAR * const restrict x, const int incx,
    const GBMV_SCALAR beta,           GBMV_SCALAR * const restrict y, const int incy)
{
    // Logic from http://www.netlib.org/blas/dgbmv.f revised
    // for zero-indexed, vectorization-ready code on GCC and Intel.

    trans = toupper(trans);  // Simplifies case-insensitive checks

    // Test the input parameters
    if (trans != 'N' && trans != 'T' && trans != 'C') {
        return 1;
    } else if (m < 0) {
        return 2;
    } else if (n < 0) {
        return 3;
    } else if (kl < 0) {
        return 4;
    } else if (ku < 0) {
        return 5;
    } else if (lda <= kl + ku) {
        return 8;
    } else if (incx == 0) {
        return 10;
    } else if (incy == 0) {
        return 13;
    }

    // Quick return if possible
    if (m == 0 || n == 0 || (alpha == 0 && beta == 1)) {
        return 0;
    }

    // Set the lengths of the vectors x and y.
    const int lenx = trans == 'N' ? n : m;
    const int leny = trans == 'N' ? m : n;

    // Start the operations. In this version the elements of A are accessed
    // sequentially with one pass through the band part of A.

    // First form y := beta*y (ignoring the irrelevant incy sign)...
    if (beta != 1) {
        const int abs_incy = abs(incy);
        if (abs_incy == 1) {  // ...for contiguous y
            if (beta == 0) {
                memset(y, 0, leny*sizeof(y[0]));
            } else {
                for (int i = 0; i < leny; ++i) {
                    y[i] *= beta;
                }
            }
        } else {              // ...for strided y
            if (beta == 0) {
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = 0; i < leny; ++i) {
                    y[i*abs_incy] = 0;
                }
            } else {
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = 0; i < leny; ++i) {
                    y[i*abs_incy] *= beta;
                }
            }
        }
    }

    // Quick return when the matrix used is irrelevant
    if (alpha == 0) {
        return 0;
    }

    // Set up the start points in x and y.
    int kx = incx > 0 ? 0 : incx*(1 - lenx);
    int ky = incy > 0 ? 0 : incy*(1 - leny);

    // Banded matrix dereference is always of form a[ku + i + j*(lda - 1)].
    // Incorporate the ku offset and decrement lda to speed indexing in loops.
    a += ku;
    --lda;

    // Perform the banded matrix-vector accumulation
    const int klp1 = kl + 1;

    if (trans == 'N') {   // Form y := alpha*A*x + y...

        int jx = kx;
        if (incy == 1) {  // ...for contiguous y
            for (int j = 0; j < n; ++j) {
                const int il = ndxmax(0, j - ku);
                const int iu = ndxmin(m, j + klp1);
                const GBMV_SCALAR temp = alpha*x[jx];
                for (int i = il; i < iu; ++i) {
                    y[i] += temp*a[i];
                }
                jx += incx;
                a += lda;
            }
        } else {          // ...for strided y
            for (int j = 0; j < n; ++j) {
                const int il = ndxmax(0, j - ku);
                const int iu = ndxmin(m, j + klp1);
                const GBMV_SCALAR temp = alpha*x[jx];
                int iy = ky;
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = il; i < iu; ++i) {
                    y[iy] += temp*a[i];
                    iy += incy;
                }
                jx += incx;
                if (j >= ku) {
                    ky += incy;
                }
                a += lda;
            }
        }

    } else {              // Form y := alpha*A**T*x + y...

        int jy = ky;
        if (incx == 1) {  // ...for contiguous x
            for (int j = 0; j < n; ++j) {
                const int il = ndxmax(0, j - ku);
                const int iu = ndxmin(m, j + klp1);
                GBMV_SCALAR temp = 0;
                for (int i = il; i < iu; ++i) {
                    temp += a[i]*x[i];
                }
                y[jy] += alpha * temp;
                jy += incy;
                a += lda;
            }
        } else {          // ...for strided x
            for (int j = 0; j < n; ++j) {
                const int il = ndxmax(0, j - ku);
                const int iu = ndxmin(m, j + klp1);
                GBMV_SCALAR temp = 0;
                int ix = kx;
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = il; i < iu; ++i) {
                    temp += a[i]*x[ix];
                    ix += incx;
                }
                y[jy] += alpha*temp;
                jy += incy;
                if (j >= ku) {
                    kx += incx;
                }
                a += lda;
            }
        }
    }

    return 0;
}

#undef GBMV_STATIC
#undef GBMV_FUNCTION
#undef GBMV_COMPONENT
#undef GBMV_SCALAR
#undef GBMV_KL
#undef GBMV_KU

The same tricks transfer over to dsbmv as well (though the studious will notice indexing differences relative to the Netlib version):

// #include-time parameters available followed by a sample usage.
// Each of these are #undef-ed at the end of this template file
//
// #define SBMV_STATIC    /*empty*/         /* Use for static function    */
// #define SBMV_FUNCTION  dsbmv             /* Function name              */
// #define SBMV_COMPONENT double            /* Type of matrices           */
// #define SBMV_SCALAR    double            /* Type of coeffs and vectors */
// #define SBMV_K         const int k,      /* Matrix bandwidth           */

#ifndef __SBMV_INTERNAL_ONCE
#define __SBMV_INTERNAL_ONCE
static inline int ndxmin(int a, int b) { return a < b ? a : b; }
static inline int ndxmax(int a, int b) { return a > b ? a : b; }
#endif /* __SBMV_INTERNAL_ONCE */

SBMV_STATIC int SBMV_FUNCTION(
    char uplo, const int n, SBMV_K
    const SBMV_SCALAR alpha, const SBMV_COMPONENT *       restrict a, int lda,
                             const    SBMV_SCALAR * const restrict x, const int incx,
    const SBMV_SCALAR beta,           SBMV_SCALAR * const restrict y, const int incy)
{
    // Logic from http://www.netlib.org/blas/dsbmv.f revised
    // for zero-indexed, vectorization-ready code on GCC and Intel.

    uplo = toupper(uplo);  // Simplifies case-insensitive checks

    // Test the input parameters
    if (uplo != 'U' && uplo != 'L') {
        return 1;
    } else if (n < 0) {
        return 2;
    } else if (k < 0) {
        return 3;
    } else if (lda <= k) {
        return 6;
    } else if (incx == 0) {
        return 8;
    } else if (incy == 0) {
        return 11;
    }

    // Quick return if possible
    if (n == 0 || (alpha == 0 && beta == 1)) {
        return 0;
    }

    // Start the operations. In this version the elements of A are accessed
    // sequentially with one pass through the band part of A.

    // First form y := beta*y (ignoring the irrelevant incy sign)...
    if (beta != 1) {
        const int abs_incy = abs(incy);
        if (abs_incy == 1) {  // ...for contiguous y
            if (beta == 0) {
                memset(y, 0, n * sizeof(y[0]));
            } else {
                for (int i = 0; i < n; ++i) {
                    y[i] *= beta;
                }
            }
        } else {              // ...for strided y
            if (beta == 0) {
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = 0; i < n; ++i) {
                    y[i*abs_incy] = 0;
                }
            } else {
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = 0; i < n; ++i) {
                    y[i*abs_incy] *= beta;
                }
            }
        }
    }

    // Quick return when the matrix used is irrelevant
    if (alpha == 0) {
        return 0;
    }

    // Set up the start points in x and y.
    int kx = incx > 0 ? 0 : incx*(1 - n);
    int ky = incy > 0 ? 0 : incy*(1 - n);

    // Perform the banded matrix-vector accumulation

    if (uplo == 'U') {  // Form y := alpha*A*x + y for A upper storage...

        // Banded matrix dereference is always of form a[ku + i + j*(lda - 1)].
        // Incorporate the ku offset and decrement lda to speed indexing.
        a += k;
        --lda;

        if (incx == 1 && incy == 1) {  // ...for contiguous x and y
            for (int j = 0; j < n; ++j) {
                const int il = ndxmax(0, j - k);
                const SBMV_SCALAR temp1 = alpha*x[j];
                SBMV_SCALAR       temp2 = 0;
                for (int i = il; i < j; ++i) {
                    y[i]  += temp1*a[i];
                    temp2 += a[i]*x[i];
                }
                y[j] += temp1*a[j] + alpha*temp2;
                a += lda;
            }
        } else {                       // ...for strided x and/or strided y
            int jx = kx;
            int jy = ky;
            for (int j = 0; j < n; ++j) {
                const int il = ndxmax(0, j - k);
                const SBMV_SCALAR temp1 = alpha*x[jx];
                SBMV_SCALAR temp2       = 0;
                int ix = kx;
                int iy = ky;
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = il; i < j; ++i) {
                    y[iy] += temp1*a[i];
                    temp2 += a[i]*x[ix];
                    ix += incx;
                    iy += incy;
                }
                y[jy] += temp1*a[j] + alpha*temp2;
                jx += incx;
                jy += incy;
                if (j >= k) {
                    kx += incx;
                    ky += incy;
                }
                a += lda;
            }
        }

    } else {            // Form y := alpha*A*x + y for A lower storage...

        // Banded matrix dereference is always of form a[ku + i + j*(lda - 1)]
        // where here ku == 0.  Decrement lda to speed indexing.
        --lda;

        const int kp1 = k + 1;
        if (incx == 1 && incy == 1) {  // ...for contiguous x and y
            for (int j = 0; j < n; ++j) {
                const int iu = ndxmin(n, j + kp1);
                const SBMV_SCALAR temp1 = alpha*x[j];
                SBMV_SCALAR       temp2 = 0;
                y[j] += temp1*a[j];
                for (int i = j + 1; i < iu; ++i) {
                    y[i]  += temp1*a[i];
                    temp2 += a[i]*x[i];
                }
                y[j] += alpha*temp2;
                a += lda;
            }
        } else {                       // ...for strided x and/or strided y
            int jx = kx;
            int jy = ky;
            for (int j = 0; j < n; ++j) {
                const int iu = ndxmin(n, j + kp1);
                const SBMV_SCALAR temp1 = alpha*x[jx];
                SBMV_SCALAR temp2       = 0;
                y[jy] += temp1*a[j];
                int ix = jx;
                int iy = jy;
#ifdef __INTEL_COMPILER
#pragma unroll
#endif
                for (int i = j + 1; i < iu; ++i) {
                    ix += incx;
                    iy += incy;
                    y[iy] += temp1*a[i];
                    temp2 += a[i]*x[ix];
                }
                y[jy] += alpha*temp2;
                jx += incx;
                jy += incy;
                a += lda;
            }
        }

    }

    return 0;
}

#undef SBMV_STATIC
#undef SBMV_FUNCTION
#undef SBMV_COMPONENT
#undef SBMV_SCALAR
#undef SBMV_K

Hopefully these will be useful to someone else.

Update 11 January 2012:Of course, now that I've found the reference XBLAS, pretty much everything here could (should?) have been done within XBLAS' macro framework. XBLAS definitely a cleaner mixed precision approach.

GNU parallel for parameter sweeps

2011-09-08T16:14:00.000Z

Following up on previous failures, I've chosen to test a portion of my compressible channel code Suzerain across a small parameter sweep. This problem is ill-suited for a traditional HPC resource as many of the jobs will fail and none of them have to run for very long.

The problem is perfect for GNU parallel:

#!/bin/bash

# Build the fixed parameters for channel initialization
fixed="--Ma=1.5 --Re=3000 --Pr=0.7 --gamma=1.4 --beta=0.7 --k=8"

# Build the parameter sweep affix for GNU parallel, placeholders, and a case description
sweep=""                       ; case="case{#}"                      ; args="$case"
sweep="$sweep ::: 96 120"      ; Nx="--Nx={1}"                       ; args="$args $Nx"
sweep="$sweep ::: 72  96"      ; Ny="--Ny={2}"                       ; args="$args $Ny"
sweep="$sweep :::  2   3"      ; htdelta="--htdelta={3}"             ; args="$args $htdelta"
sweep="$sweep ::: 60  96"      ; Nz="--Nz={4}"                       ; args="$args $Nz"
sweep="$sweep ::: 125 150"     ; fluct_percent="--fluct_percent={5}" ; args="$args $fluct_percent"
sweep="$sweep ::: 12345 67890" ; fluct_seed="--fluct_seed={6}"       ; args="$args $fluct_seed"
sweep="$sweep :::  2   5"      ; alpha="--alpha={7}"                 ; args="$args $alpha"
sweep="$sweep ::: .4  .5"      ; npower="--npower={8}"               ; args="$args $npower"
sweep="$sweep ::: 0:1 .5:.75"  ; fluct_kxfrac="--fluct_kxfrac={9}"   ; args="$args $fluct_kxfrac"
sweep="$sweep ::: 0:1 .5:.75"  ; fluct_kzfrac="--fluct_kzfrac={10}"  ; args="$args $fluct_kzfrac"

# Remove any stale data
rm -rfv case* bad.case*

# Generate a master table of test cases for fun and profit
parallel -k -j 1 echo $args $sweep | column -t > master

# Generate inputs for a single test case
parallel "mkdir $case && cd $case && echo $args > args"                                                             \
     " && ../bin/channel_init initial.h5 $fixed $Nx $Ny $htdelta $Nz $alpha $npower"                                \
     " && ../bin/channel_explicit initial.h5 --advance_nt=0 $fluct_percent $fluct_seed $fluct_kxfrac $fluct_kzfrac" \
     $sweep

# Run each case for 1000 time steps to shake out non-robust scenarios
parallel --eta -j 3 -u                                                              \
     "cd {} && ../bin/channel_explicit restart0.h5 --advance_nt=1000 --status_nt=5" \
     ::: case*
 
# Flag test cases that glaringly failed with a ".bad" suffix
parallel "mv {//} bad.{//}" ::: $(grep -l "TimeController halted" case*/log.dat)

The first two "blocks" define the parameter sweep using parallel's feature where "::: a b ::: 1 2" turns into the outer product "a 1", "a 2", "b 1", "b 2". The next bit builds a table of these outer products so I can refer to it later. The remaining three bits perform some IO-heavy initialization, some compute-intensive tasks, and then a post-processing pass on each of the entries in said outer product. This logic, and the associated runtime process juggling to get nice batch throughput, would be hell without GNU parallel. Thanks Ole.

ça marche (mal?)

2011-08-11T15:21:00.001Z

Suzerain, my pseudospectral, compressible DNS code, is coming along. After some manufactured solution love I'm onto channel problems. With some added solenoidal velocity perturbations my laminar solutions transitioned to turbulent ones quite nicely.

However, I'm currently puzzling over why my Coleman-like case is showing stationary mean state behavior...

...but momentum fluctuations that grow in time:

Currently suspect things are:

Inadequate resolution ('tis quite lousy)
Pressure growing in time (a probable symptom I need to check but not a root cause)
Funkiness associated with the discrete mass conservation fixes (my numerics are not conservative and holding mass constant requires some integral constraint futzing)
Inadequate dealiasing (unlikely, currently using enough to handle quadratic nonlinearities in my Fourier directions but my equations are worse-than-cubic).

Thank you to Nick Malaya for poring over plots with me this morning.

Welcome to a dumber planet.

2011-08-09T02:52:00.002Z

Consider IBM's stock price over the last several years...

...and then today's news that IBM Terminated The Blue Waters Contract because of financial, rather than technical, reasons:

When HPCwire spoke with Herb Schultz, marketing manager for IBM's Deep Computing unit, last year, he outlined a new business model that would apply a lot more scrutiny to how the company positioned its high-end supercomputers. "There is really no appetite in IBM anymore -- with some of the leadership changes over the last few years -- for revenue that has no profit with it," he told us back in November 2010.... From NCSA's perspective, the system met all of its technical requirements. In particular, they appeared confident the machine, based on Power 755 servers, would indeed be able to deliver a sustained petaflop from its 10-petaflop peak performance.... The most curious aspect of the IBM pull-out is that they had already delivered three racks of Power 755 servers to NCSA.

It is a sad day when a screamingly profitable U.S. company reneges on a contract so fundamental to the next decade of basic U.S. science. Especially a contract where even the worst imaginable outcome for IBM has no potential to noticeably impact shareholder returns for even a single quarter.

FFTW 3.3 MPI Cheat Sheet

2011-07-09T23:11:00.012Z

Because I believe they provide a nice way to perform P3DFFT- and 2DECOMP&FFT-like MPI-parallel pencil decompositions, I have have spent some time staring at FFTW's 3.3 alpha and beta releases. In particular, poring over their Distributed-memory FFTW with MPI capabilities.

Using 3.3-alpha1 I wrote a small library (called underling) which mimicked P3DFFT 2.3's data movement capabilities. I wholly isolated the data movement from computing the FFTs. This, unlike P3DFFT and 2DECOMP, provides well-defined memory layout at any stage during the pencil transposes. Functionally my first approach was quite sound. However, it performs suboptimal on-node memory reshuffling (a design mistake on my part) and can stand to be improved.

I am revisiting the assumption of separating the parallel FFTs from the parallel data movement. Doing so allows using the higher-level 2D r2c/c2r planning APIs freshly documented for 3.3-beta. It should require less needless memory reshuffling when combined with appropriate FFTW_MPI_TRANPOSED_IN and FFTW_MPI_TRANSPOSED_OUT flag choices.

Because, though wonderfully written, the relevant sections of the FFTW MPI documentation do not provide a cheat sheet, I have cooked my own for the various piece parts I may use in a second attempt. Perhaps someone will find it useful. Be sure to check the FFTW MPI reference for full details, especially the local_size calls necessary to obtain data distribution information.

Please tell me if you catch any mistakes. Many thanks to Steven G. Johnson for correcting my earlier misunderstanding of the real-to-complex and complex-to-real 2D DFT semantics.

							transposed in					transposed out					transposed in\|out
transpose	in	n0/P	×	n1	×	nc	n1	×	n0/P	×	nc	n0/P	×	n1	×	nc	n1	×	n0/P	×	nc
transpose	out	n1/P	×	n0	×	nc	n1/P	×	n0	×	nc	n0	×	n1/P	×	nc	n0	×	n1/P	×	nc
c2c 2D	in	ñ0/P	×	ñ1	×	ñc	ñ1/P	×	ñ0	×	ñc	ñ0/P	×	ñ1	×	ñc	ñ1/P	×	ñ0	×	ñc
c2c 2D	out	ñ0/P	×	ñ1	×	ñc	ñ0/P	×	ñ1	×	ñc	ñ1/P	×	ñ0	×	ñc	ñ1/P	×	ñ0	×	ñc
r2c 2D	in	n0/P	×	2(n1/2+1)	×	nc	2(n1/2+1)/P	×	n0	×	nc	n0/P	×	2(n1/2+1)	×	nc	2(n1/2+1)/P	×	n0	×	nc
r2c 2D	out	ñ0/P	×	(ñ1/2+1)	×	ñc	ñ0/P	×	(ñ1/2+1)	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc
c2r 2D	in	ñ0/P	×	ñ1/2+1	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc	ñ0/P	×	ñ1/2+1	×	ñc	(ñ1/2+1)/P	×	ñ0	×	ñc
c2r 2D	out	n0/P	×	2(n1/2+1)	×	nc	n0/P	×	2(n1/2+1)	×	nc	2(n1/2+1)/P	×	n0	×	nc	2(n1/2+1)/P	×	n0	×	nc

Notation: Real-valued directions are denoted n0 and n1 while complex-valued directions are ñ0 and ñ1. Half-complex storage is denoted ñ/2+1 and its padded, real-valued counterpart is 2(n/2+1). Directions decomposed along a communicator with P processes are denoted n/P. nc stands for "number of components" and corresponds to the advanced planning API's howmany arguments.

My divergence divergence.

2011-07-05T23:22:00.000Z

I've spent three full work days chasing down a verification failure within my spectral, compressible Navier—Stokes code. I can't remember the last time I was so happy to average 1/3rd of a source line per day.

Using Boost Spirit 2.1+ to evaluate constant arithmetic expressions

2011-06-26T15:56:00.005Z

I used Boost Program_options for my research code's input parsing requirements. As expected, it works well. For floating point quantities, I've found having to specify 4.18879020478639 for 4π/3 is a bit cumbersome. When specifying arguments on the command line, I could use any one of a number of hacks to perform floating point arithmetic in bash. But that didn't help me when taking input from response files. I wanted an arithmetic expression utility that could process a C-like grammar. Long double accuracy was a requirement. I looked around a bit. GNU libmatheval and muParser appeared to be good fits but looked to be double precision-only. Though straightforward, I decided against combining the Shunting-yard algorithm with some in-place RPN logic because I wasn't enthusiastic about either finding or writing an appropriate (and admittedly toy simple) lexer.

Already having a Boost pre-requisite in my code, I decided to try out Boost Spirit since it could provide a header-only, precision-agnostic solution. Spirit has gone through a lot of significant revisions and many articles seemed to be written for the older Spirit Classic. The latest official reference manual gave several excellent examples that almost fit. I adapted Spirit's calc3.cpp and roman.cpp sample grammars to accomplish what I wanted. I added ** as an operator for exponentiation and pi as a built-in constant. I also added a smattering of special functions (e.g. exp). The task required about 150 lines of Spirit-based code. Further extensions should be fairly straightforward.
My sample code does the following:

Declares X-macros for the unary and binary functions built into the grammar.
Declares Boost Phoenix-ready functors for the unary and binary functions.
Declares expression::grammar consisting of the rules expression, term, factor, and primary.
Declares a helper called parse to hide needing boost::spirit::ascii::space when performing whitespace-agnostic parsing.
Declares a small test macro called EXPRTEST to reduce test case boilerplate.
Runs a collection of tests to ensure the parser is behaving as expected.

Compared to the samples on which I built, the biggest changes to the grammar were adding primary and reworking factor. I could have shortened up some details with a few using boost::spirit::qi; declarations, but I wanted to keep track of exactly which parts of Spirit I used. Feedback and suggestions are most welcome.

A less hideous version is presented further below.

///////////////////////////////////////////////
// Ensure Boost Spirit version is high enough.
// Compilation errors are astounding otherwise.
///////////////////////////////////////////////

#include <boost/version.hpp>

#if BOOST_VERSION < 104100
#error "Require at least Spirit 2.1 from at least Boost 1.41"
#endif

//////////////////////////////////////////////////////////////////
// Implementation uses X macros to reduce boilerplate.
// See http://drdobbs.com/blogs/cpp/228700289 for X macro details.
//////////////////////////////////////////////////////////////////

#include <cmath>
#include <limits>

#include <boost/math/constants/constants.hpp>
#include <boost/spirit/include/phoenix.hpp>
#include <boost/spirit/include/qi.hpp>

namespace expression {

#define EXPRESSION_FORALL_UNARY_FUNCTIONS(apply)  \
    apply(abs)                                    \
    apply(acos)                                   \
    apply(asin)                                   \
    apply(atan)                                   \
    apply(ceil)                                   \
    apply(cos)                                    \
    apply(cosh)                                   \
    apply(exp)                                    \
    apply(floor)                                  \
    apply(log)                                    \
    apply(log10)                                  \
    apply(sin)                                    \
    apply(sinh)                                   \
    apply(sqrt)                                   \
    apply(tan)                                    \
    apply(tanh)

#define EXPRESSION_FORALL_BINARY_FUNCTIONS(apply) \
    apply(pow)                                    \
    apply(atan2)


namespace detail {

#define DECLARE_UNARY_FUNCTOR(name)                          \
struct name ## _                                             \
{                                                            \
    template <typename X> struct result { typedef X type; }; \
                                                             \
    template <typename X> X operator()(X x) const            \
    {                                                        \
        using namespace std;                                 \
        return name(x);                                      \
    }                                                        \
};
EXPRESSION_FORALL_UNARY_FUNCTIONS(DECLARE_UNARY_FUNCTOR)
#undef DECLARE_UNARY_FUNCTOR

#define DECLARE_BINARY_FUNCTOR(name)                                     \
struct name ## _                                                         \
{                                                                        \
    template <typename X, typename Y> struct result { typedef X type; }; \
                                                                         \
    template <typename X, typename Y> X operator()(X x, Y y) const       \
    {                                                                    \
        using namespace std;                                             \
        return name(x, y);                                               \
    }                                                                    \
};
EXPRESSION_FORALL_BINARY_FUNCTIONS(DECLARE_BINARY_FUNCTOR)
#undef DECLARE_BINARY_FUNCTOR

} // end namespace detail

template <typename FPT, typename Iterator>
struct grammar
    : boost::spirit::qi::grammar<
            Iterator, FPT(), boost::spirit::ascii::space_type
        >
{
    boost::spirit::qi::rule<
            Iterator, FPT(), boost::spirit::ascii::space_type
        > expression, term, factor, primary;

    grammar() : grammar::base_type(expression)
    {
        namespace constants = boost::math::constants;
#define DECLARE_PHOENIX_FUNCTION(name)                    \
        boost::phoenix::function<detail::name ## _> name;
EXPRESSION_FORALL_UNARY_FUNCTIONS(DECLARE_PHOENIX_FUNCTION)
EXPRESSION_FORALL_BINARY_FUNCTIONS(DECLARE_PHOENIX_FUNCTION)
#undef DECLARE_PHOENIX_FUNCTION

        using boost::spirit::qi::real_parser;
        using boost::spirit::qi::real_policies;
        real_parser<FPT,real_policies<FPT> > real;

        using boost::spirit::qi::_1;
        using boost::spirit::qi::_2;
        using boost::spirit::qi::no_case;
        using boost::spirit::qi::_val;

        expression =
            term                   [_val =  _1]
            >> *(  ('+' >> term    [_val += _1])
                |  ('-' >> term    [_val -= _1])
                )
            ;

        term =
            factor                 [_val =  _1]
            >> *(  ('*' >> factor  [_val *= _1])
                |  ('/' >> factor  [_val /= _1])
                )
            ;

        factor =
            primary                [_val =  _1          ]
            >> *(  ("**" >> factor [_val = pow(_val, _1)])
                )
            ;

#define STRINGIFY(s) #s
#define UNARY_FUNCTION_RULE(name)                                             \
    | no_case[STRINGIFY(name)] >> '(' >> expression [_val =  name(_1)] >> ')'
#define BINARY_FUNCTION_RULE(name)                                            \
    | ( no_case[STRINGIFY(name)]>> '('                                        \
        >> expression >> ',' >> expression >> ')' )[_val =  name(_1,_2)]

        primary =
            real                   [_val =  _1]
            |   '(' >> expression  [_val =  _1] >> ')'
            |   ('-' >> primary    [_val = -_1])
            |   ('+' >> primary    [_val =  _1])
            |   no_case["pi"]      [_val =  constants::pi<FPT>()]
                EXPRESSION_FORALL_UNARY_FUNCTIONS(UNARY_FUNCTION_RULE)
                EXPRESSION_FORALL_BINARY_FUNCTIONS(BINARY_FUNCTION_RULE)
            ;

#undef UNARY_RULE
#undef BINARY_RULE
#undef STRINGIFY

    }
};

#undef EXPRESSION_FORALL_UNARY_FUNCTIONS
#undef EXPRESSION_FORALL_BINARY_FUNCTIONS

template <typename FPT, typename Iterator>
bool parse(Iterator &iter,
           Iterator end,
           const grammar<FPT,Iterator> &g,
           FPT &result)
{
    return boost::spirit::qi::phrase_parse(
                iter, end, g, boost::spirit::ascii::space, result);
}

} // end namespace expression


/////////////
// Test cases
/////////////

#include <string>
#include <boost/math/special_functions/fpclassify.hpp>

#define BOOST_TEST_MODULE expression
#include <boost/test/included/unit_test.hpp>

const expression::grammar<double,std::string::const_iterator> eg;
const double close_enough = std::numeric_limits<double>::epsilon();

#define EXPRTEST(casename, expr, expected)               \
BOOST_AUTO_TEST_CASE( casename )                         \
{                                                        \
    const std::string s = expr;                          \
    std::string::const_iterator iter = s.begin();        \
    std::string::const_iterator end  = s.end();          \
    double result;                                       \
    BOOST_REQUIRE(parse(iter, end, eg, result));         \
    BOOST_CHECK(iter == end);                            \
    BOOST_CHECK_CLOSE(result, (expected), close_enough); \
}

EXPRTEST(literal1, "1.234", 1.234)
EXPRTEST(literal2, "4.2e2", 420)
EXPRTEST(literal3, "5e-01", 0.5)
EXPRTEST(literal4, "-3",    -3)
EXPRTEST(literal5, "pi",    boost::math::constants::pi<double>())

EXPRTEST(basicop1, " 2 +\t3\n",   5)       // Whitespace ignored
EXPRTEST(basicop2, " 2 -\t3\n",  -1)
EXPRTEST(basicop3, " 2 *\t3\n",   6)
EXPRTEST(basicop4, " 2 /\t3\n",   2./3.)   // Double division
EXPRTEST(basicop5, " 2 ** 3\n",  8)

EXPRTEST(pemdas1,  "2*3+4*5",   2*3+4*5)
EXPRTEST(pemdas2,  "2*(3+4)*5", 2*(3+4)*5)
EXPRTEST(pemdas3,  "2**3+4",    std::pow(2,3)+4)
EXPRTEST(pemdas4,  "2**2**-3",  std::pow(2.,std::pow(2.,-3.)))

EXPRTEST(unary1,   "-(2)",      -2)
EXPRTEST(unary2,   "-(-2)",      2)
EXPRTEST(unary3,   "+(-2)",     -2)
EXPRTEST(unary4,   "+(+2)",      2)

EXPRTEST(func_abs,   "abs (-1.0)", std::abs(-1.0))
EXPRTEST(func_acos,  "acos( 1.0)", std::acos(1.0))
EXPRTEST(func_asin,  "asin( 1.0)", std::asin(1.0))
EXPRTEST(func_atan,  "atan( 1.0)", std::atan(1.0))
EXPRTEST(func_ceil,  "ceil( 0.5)", std::ceil(0.5))
EXPRTEST(func_cos,   "cos ( 1.0)", std::cos(1.0))
EXPRTEST(func_cosh,  "cosh( 1.0)", std::cosh(1.0))
EXPRTEST(func_exp,   "exp ( 1.0)", std::exp(1.0))
EXPRTEST(func_floor, "floor(0.5)", std::floor(0.5))
EXPRTEST(func_log,   "log ( 1.0)", std::log(1.0))
EXPRTEST(func_log10, "log10(0.5)", std::log10(0.5))
EXPRTEST(func_sin,   "sin ( 1.0)", std::sin(1.0))
EXPRTEST(func_sinh,  "sinh( 1.0)", std::sinh(1.0))
EXPRTEST(func_sqrt,  "sqrt( 1.0)", std::sqrt(1.0))
EXPRTEST(func_tan,   "tan ( 1.0)", std::tan(1.0))
EXPRTEST(func_tanh,  "tanh( 1.0)", std::tanh(1.0))

EXPRTEST(func_pow,   "pow  (2.0, 3.0)", std::pow(2.0,3.0))
EXPRTEST(func_atan2, "atan2(2.0, 3.0)", std::atan2(2.0,3.0))

BOOST_AUTO_TEST_CASE( parse_inf )
{
    const std::string s = "inf";
    std::string::const_iterator iter = s.begin();
    std::string::const_iterator end = s.end();
    double result;
    BOOST_REQUIRE(parse(iter, end, eg, result));
    BOOST_CHECK(iter == end);
    BOOST_CHECK((boost::math::isinf)(result));
}

BOOST_AUTO_TEST_CASE( parse_infinity )
{
    const std::string s = "infinity";
    std::string::const_iterator iter = s.begin();
    std::string::const_iterator end = s.end();
    double result;
    BOOST_REQUIRE(parse(iter, end, eg, result));
    BOOST_CHECK(iter == end);
    BOOST_CHECK((boost::math::isinf)(result));
}

BOOST_AUTO_TEST_CASE( parse_nan )
{
    const std::string s = "nan";
    std::string::const_iterator iter = s.begin();
    std::string::const_iterator end = s.end();
    double result;
    BOOST_REQUIRE(parse(iter, end, eg, result));
    BOOST_CHECK(iter == end);
    BOOST_CHECK((boost::math::isnan)(result));
}

Updated 9 Jan 2012: During a second pass motivated by reducing compilation times, I removed the X-macro nastiness and replaced it with some qi::symbols-based logic. The code in the anonymous namespace could likely be removed if I could figure out the right boost::phoenix::bind hoodoo to use to replace lazy_ufunc_ and lazy_bfunc_ (any suggestions?). I have not yet measured the difference in either memory footprint or runtime performance. Compilation time and memory is greatly reduced for both GNU and Intel compilers. The same unit tests from before continue to pass against the following version:

#include <cmath>
#include <limits>
#include <iterator>

#include <boost/spirit/version.hpp>
#if !defined(SPIRIT_VERSION) || SPIRIT_VERSION < 0x2010
#error "At least Spirit version 2.1 required"
#endif
#include <boost/math/constants/constants.hpp>
#include <boost/spirit/include/phoenix.hpp>
#include <boost/spirit/include/qi.hpp>

namespace expression {

namespace { // anonymous

struct lazy_pow_
{
    template <typename X, typename Y>
    struct result { typedef X type; };

    template <typename X, typename Y>
    X operator()(X x, Y y) const
    {
        return std::pow(x, y);
    }
};

struct lazy_ufunc_
{
    template <typename F, typename A1>
    struct result { typedef A1 type; };

    template <typename F, typename A1>
    A1 operator()(F f, A1 a1) const
    {
        return f(a1);
    }
};

struct lazy_bfunc_
{
    template <typename F, typename A1, typename A2>
    struct result { typedef A1 type; };

    template <typename F, typename A1, typename A2>
    A1 operator()(F f, A1 a1, A2 a2) const
    {
        return f(a1, a2);
    }
};

} // end namespace anonymous

template <typename FPT, typename Iterator>
struct grammar
    : boost::spirit::qi::grammar<
            Iterator, FPT(), boost::spirit::ascii::space_type
        >
{

    // symbol table for constants like "pi"
    struct constant_
        : boost::spirit::qi::symbols<
                typename std::iterator_traits<Iterator>::value_type,
                FPT
            >
    {
        constant_()
        {
            this->add
                ("pi", boost::math::constants::pi<FPT>()  )
            ;
        }
    } constant;

    // symbol table for unary functions like "abs"
    struct ufunc_
        : boost::spirit::qi::symbols<
                typename std::iterator_traits<Iterator>::value_type,
                FPT (*)(FPT)
            >
    {
        ufunc_()
        {
            this->add
                ("abs"   , (FPT (*)(FPT)) std::abs  )
                ("acos"  , (FPT (*)(FPT)) std::acos )
                ("asin"  , (FPT (*)(FPT)) std::asin )
                ("atan"  , (FPT (*)(FPT)) std::atan )
                ("ceil"  , (FPT (*)(FPT)) std::ceil )
                ("cos"   , (FPT (*)(FPT)) std::cos  )
                ("cosh"  , (FPT (*)(FPT)) std::cosh )
                ("exp"   , (FPT (*)(FPT)) std::exp  )
                ("floor" , (FPT (*)(FPT)) std::floor)
                ("log"   , (FPT (*)(FPT)) std::log  )
                ("log10" , (FPT (*)(FPT)) std::log10)
                ("sin"   , (FPT (*)(FPT)) std::sin  )
                ("sinh"  , (FPT (*)(FPT)) std::sinh )
                ("sqrt"  , (FPT (*)(FPT)) std::sqrt )
                ("tan"   , (FPT (*)(FPT)) std::tan  )
                ("tanh"  , (FPT (*)(FPT)) std::tanh )
            ;
        }
    } ufunc;

    // symbol table for binary functions like "pow"
    struct bfunc_
        : boost::spirit::qi::symbols<
                typename std::iterator_traits<Iterator>::value_type,
                FPT (*)(FPT, FPT)
            >
    {
        bfunc_()
        {
            this->add
                ("pow"  , (FPT (*)(FPT, FPT)) std::pow  )
                ("atan2", (FPT (*)(FPT, FPT)) std::atan2)
            ;
        }
    } bfunc;

    boost::spirit::qi::rule<
            Iterator, FPT(), boost::spirit::ascii::space_type
        > expression, term, factor, primary;

    grammar() : grammar::base_type(expression)
    {
        using boost::spirit::qi::real_parser;
        using boost::spirit::qi::real_policies;
        real_parser<FPT,real_policies<FPT> > real;

        using boost::spirit::qi::_1;
        using boost::spirit::qi::_2;
        using boost::spirit::qi::_3;
        using boost::spirit::qi::no_case;
        using boost::spirit::qi::_val;

        boost::phoenix::function<lazy_pow_>   lazy_pow;
        boost::phoenix::function<lazy_ufunc_> lazy_ufunc;
        boost::phoenix::function<lazy_bfunc_> lazy_bfunc;

        expression =
            term                   [_val =  _1]
            >> *(  ('+' >> term    [_val += _1])
                |  ('-' >> term    [_val -= _1])
                )
            ;

        term =
            factor                 [_val =  _1]
            >> *(  ('*' >> factor  [_val *= _1])
                |  ('/' >> factor  [_val /= _1])
                )
            ;

        factor =
            primary                [_val =  _1]
            >> *(  ("**" >> factor [_val = lazy_pow(_val, _1)])
                )
            ;

        primary =
            real                   [_val =  _1]
            |   '(' >> expression  [_val =  _1] >> ')'
            |   ('-' >> primary    [_val = -_1])
            |   ('+' >> primary    [_val =  _1])
            |   no_case[constant]  [_val =  _1]
            |   (no_case[ufunc] >> '(' >> expression >> ')')
                                   [_val = lazy_ufunc(_1, _2)]
            |   (no_case[bfunc] >> '(' >> expression >> ','
                                       >> expression >> ')')
                                   [_val = lazy_bfunc(_1, _2, _3)]
            ;

    }
};

template <typename FPT, typename Iterator>
bool parse(Iterator &iter,
           Iterator end,
           const grammar<FPT,Iterator> &g,
           FPT &result)
{
    return boost::spirit::qi::phrase_parse(
                iter, end, g, boost::spirit::ascii::space, result);
}

} // end namespace expression

Adding enthought.traits.api.HasTraits as an ancestor using Python metaclasses

2011-05-24T21:49:00.001Z

Today, while banging my head against a wall trying to get a SWIG-generated Python module to play nice with Python Traits, I cooked a way to use Python metaclasses to modify a class to add enthought.traits.has_traits.HasTraits as a base class:

from enthought.traits.api import HasTraits
from enthought.traits.has_traits import MetaHasTraits

# See http://onlamp.com/lpt/a/3388 for an intro to Python metaclasses and
# http://code.activestate.com/recipes/204197-solving-the-metaclass-conflict/
# for what can go wrong here
class AddHasTraitsAncestor(MetaHasTraits):
    def __new__(cls, name, bases, dct):
        return MetaHasTraits.__new__(cls, name, cls.modify(bases), dct)

    def __init__(cls, name, bases, dct):
        super(AddHasTraitsAncestor, cls).__init__(name, cls.modify(bases), dct)

    @staticmethod
    def modify(bases):
        if (bases == (object,)):
            return (HasTraits,)
        else:
            new_bases = list(bases)
            new_bases.insert(0, HasTraits)
            return tuple(new_bases)

Adding the line __metaclass__ = AddHasTraitsAncestor to a class definition is enough, in some cases, for the class to now have HasTraits as a base class. Most of the frustrating bits were before I realized HasTraits has its own metaclass MetaHasTraits which AddHasTraitsAncestor had to subclass before the approach worked.

"Uhhh...?" you say. That puzzled look asking "Why didn't you just add the base class yourself?". I used this approach because SWIG was generating the Python class definitions and I was abusing SWIG's pythoncode feature to tweak the SWIG wrappers to be Traits-ready a la

%extend someCPlusPlusClass {

    %pythoncode %{
        __metaclass__ = AddHasTraitsAncestor
    %}

}

Turns out it works and my SWIG wrappers are now Traits-friendly. I'm still working on getting Traits to see the SWIG-based members, however. Suggestions much appreciated.

Using gnuplot from the command line

2011-04-16T02:40:00.011Z

gnuplot is excellent for both quick 'n' dirty and lengthy plotting projects. Still, for truly quick work of just plottin' some column against some other column, I've always felt like gnuplot demands too many keystrokes:

[1000 rhys@tiresias 16]$ gnuplot

 G N U P L O T
 Version 4.2 patchlevel 6
 last modified Sep 2009
 System: Linux 2.6.32-30-generic

 Copyright (C) 1986 - 1993, 1998, 2004, 2007 - 2009
 Thomas Williams, Colin Kelley and many others

 Type `help` to access the on-line reference manual.
 The gnuplot FAQ is available from http://www.gnuplot.info/faq/

 Send bug reports and suggestions to 


Terminal type set to 'x11'
gnuplot> plot 'restart0.mean' using 2:17
gnuplot> exit

Of course, one can turn this process into a one-liner at the cost of quoting:

gnuplot -e "plot 'restart0.mean' using 2:17"

Call me picky, but the signal-to-noise ratio in that one-liner still feels too high.

After some quoting experimentation using printf's %q and %b format specifiers, the following little function definition

gplot () {
   local fileexpr=$(printf "'%q'" $1)
   shift
   local plotexpr=$(printf '%b' "plot ${fileexpr} $*")
   gnuplot -persist -raise -e "$plotexpr"
}

makes

gplot restart0.mean using 2:17

behave like the previous two examples. This function happily hides the file name quoting requirements and remains tab-completion-friendly.

Providing computed column values like gnuplot's using 1:$(cos($2)) on the command line still requires thought, but the straightforward tasks work thoughtlessly. Suggestions for improvements definitely welcome in the comments. One glaring bug is that passing a pipe expression as the file bombs.

Update 20110504: I've decided that I much prefer creating a gplot script containing

#!/bin/bash
file=$1
shift
gnuplot -persist -raise <(cat <<-GNUPLOTSCRIPT
plot '$file' $*
GNUPLOTSCRIPT
)

as it is much simpler to work through the quoting and far more extensible.

Update 20111031: Globbing, labelling, and a handful of other enhancements:

#!/bin/bash

# Fail on first error
set -e

# Create temporary files to hold gnuplot script
tmp1=`mktemp`
tmp2=`mktemp`
trap "rm -f $tmp1 $tmp2" EXIT

# Process (and then remove) command line options
# Build up any post-terminal processing in tmp2 for later use
autotitle=
title="`date +%r`: gplot $*"
terminal="set terminal x11 title '$title' enhanced persist"
cmd=plot
forexpr=
raise=
showhelp=
while getopts "3cf:ghils:re:t:x:y:z:F:SX:Y:Z:" opt
do
  case $opt in
    3) cmd=splot
       ;;
    c) autotitle=true
       ;;
    e) terminal='set term postscript eps enhanced color dashed rounded "Arial,14"'
       pause='' # Disable interactive on EPS output
       echo set output \"$OPTARG\" >> $tmp2
       ;;
    f) forexpr=" for [$OPTARG] "
       ;;
    g) echo "set grid" >> $tmp2
       ;;
    h) showhelp=0
       ;;
    i) pause='pause -1 "Plot interactive until Enter pressed.\nHit h in plot window for help.\n"'
       ;;
    l) echo "set logscale y" >> $tmp2
       ;;
    r) raise="-raise"
       ;;
    s) savescript="$OPTARG"
       ;;
    t) echo set title \"$OPTARG\" >> $tmp2
       ;;
    x) echo set xlabel \"$OPTARG\" >> $tmp2
       ;;
    y) echo set ylabel \"$OPTARG\" >> $tmp2
       ;;
    z) echo set zlabel \"$OPTARG\" >> $tmp2
       ;;
    F) pause="pause $OPTARG; replot; reread"
       ;;
    S) stdin=true
       ;;
    X) echo set xrange \[$OPTARG\] >> $tmp2
       ;;
    Y) echo set yrange \[$OPTARG\] >> $tmp2
       ;;
    Z) echo set zrange \[$OPTARG\] >> $tmp2
       ;;
  esac
done
shift $((OPTIND-1))

if [ x$showhelp != x ]; then
    cat <<-HERE
Usage: gplot [OPTION]... (FILE|EXTGLOB) GNUPLOTCMD...
Use gnuplot to plot one or more files directly from the command line.

  -3             Perform 3D plotting using gnuplot's splot command.
  -c             Populate the key using autotitling.
  -e FILE        Save an Encapsulated Postscript (eps) called FILE.
  -f FOREXPR     Prepend a 'for [FOREXPR]' to the plotting command.
  -g             Show grid lines.
  -h             Show this help message.
  -i             Interactive plotting mode.  Hit 'h' on plot for help.
  -l             Use logarithmic scale for y axis.
  -s FILE        Save the generated gnuplot as a script called FILE.
  -r             Have the window manager raise the plot window.
  -t TITLE       Set TITLE as the plot's title.
  -x XLABEL      Specify XLABEL as the x axis label.
  -y YLABEL      Specify YLABEL as the y axis label.
  -z ZLABEL      Specify ZLABEL as the z axis label.
  -F FREQUENCY   Replot the inputs every FREQUENCY seconds.
  -S             Prior to plotting, read auxililary gunplot from stdin.
  -X XLOW:XHIGH  Specify an explicit x axis range instead of autoscaling.
  -Y YLOW:YHIGH  Specify an explicit y axis range instead of autoscaling.
  -Z ZLOW:ZHIGH  Specify an explicit z axis range instead of autoscaling.

Examples (see gnuplot documentation for complete GNUPLOTCMD details):

  gplot -i foo.dat using 1:2 with linespoints
  gplot -s foo.gp -X 0:1 -Y 0:2 foo.dat using 1:2 with linespoints
  gplot -e foo.eps foo.dat using 1:2 with linespoints
  gplot -3 restart\*.dat using 1:2:3

On error, the failing gnuplot script is shown.
HERE
    exit $showhelp
fi

# Set terminal
echo "$terminal" >> $tmp1

# Slurp any settings built up during getops processing
cat $tmp2 >> $tmp1

# Obtain file(s) to plot from first argument using extended globbing
# Deliberately allow globbing to occur in this assignment
shopt -s extglob
declare -a files=($1)
shift

# Tweak autotitle based on options and incoming argument details
if [ "$autotitle" ]; then
    # Use columnheader option iff only one file provided
    if [ ${#files[@]} -eq 1 -a $cmd != splot ]; then
        echo 'set key autotitle columnheader' >> $tmp1
    else
        echo 'set key autotitle' >> $tmp1
    fi
else
    echo 'set key noautotitle' >> $tmp1
fi

# Possibly slurp standard input for further options
# FIXME Not working for 'echo foo | gplot'
if [ "$stdin" ]; then
    cat <&1 >> $tmp1
fi

# Build gnuplot script to plot possibly multiple files
declare -i count=0
for file in "${files[@]}"
do
    count+=1
    if [ $count -eq 1 ]; then
        echo "$cmd" "$forexpr" \'$file\' $* \\ >> $tmp1
    else
        echo "   "  , "$forexpr" \'$file\' $* \\ >> $tmp1
    fi
done
echo '' >> $tmp1

# If requested, save the plotting commands as an executable script
if [ "$savescript" ]; then
    echo '#!/usr/bin/env gnuplot' > "$savescript"
    cat "$tmp1" >> "$savescript"
    chmod a+x "$savescript"
fi

# If requested, add a pause command to the end of the script. Notice this was
# not added to the $savescript as it would be annoying in that context.
if [ "$pause" ]; then
    echo "$pause" >> $tmp1
fi

# Generate the plot
# Display the script when an error occurs to aid debugging
gnuplot $raise $tmp1 || cat $tmp1

Noise.

2011-04-01T02:04:00.000Z

Following up on some quasi-1D success I'm ready to calculate some 3D fields. I seed them by taking my 1D solution and add noise. This is purportedly a terrible way to kick off turbulence as the fluctuations have to be massive to ultimately take and transition into a correct field. I've been warned...

Needing a way to monitor the fluctuation strengths and being in coefficient-space for a mixed Fourier/B-Spline discretization, I'm trying out the L^2 norm of the fluctuating signal. By Parseval's theorem along with some minor linear algebra, this turns into a sum over the complex-valued coefficient magnitudes: Computation of the L^2 norm in a mixed Fourier/B-spline collocation discretization

Now, for the first case of the L^2 norm over all three directions I need the matrix M which computes the L^2 inner product of two functions given their B-spline coefficient expansions. The code to do this atop some old contributions to the GSL isn't bad, but testing that my banded derivative operator coefficients are correct is awful.

Happily, Mathematica 7 includes BSplineBasis which I can use to cook analytic results to test my implementation. Unhappily, BSplineBasis has issues computing derivatives and an awful amount of mucking around is required to get my analytical test cases: Finding Galerkin L^2-based Operators for B-spline discretizations

In the end, this is all The Right Thing ™ to do but man it eats time. Recognizing this particular rabbit hole and how often such things ensnare my little mind, I sent the following email to my wonderful wife today:

In case you ever wonder what I do with my days....

I stare at things like this for far too long to work out the math.

Then I sit and stare at the math for far too long to figure out the computer implementation.

Then I stare at the computer implementation to test it and clean it up.

Then I use the implementation at a single line within my thesis research code.

Then I go talk to my advisor on Mondays and report that a single line now works this week that didn't work last week.

— Rhys

Steppin'

2011-03-10T23:54:00.001Z

The code I've been working on for, quite literally, years is finally coming together. It's called Suzerain. Suzerain is a mixed Fourier, B-spline spectral DNS code for solving the compressible Navier–Stokes equations for a perfect gas. Nearly every technical thing I've written on this blog has been about it in some way.

My first sanity check has been advancing a quasi-1D, laminar channel flow profile to a steady state condition. Some plots follow with units nondimensionalized by the bulk velocity, channel half width, and wall temperature and pressure. The initial conditions were constant density and temperature along with a parabolic velocity profile. Happily, the simulation relaxes the way I expect.

Onward.

Converting to and from human-readable byte counts

2011-01-25T18:40:00.000Z

Converting some number of bytes (say 1024) into a human-readable byte count (say "1K") seems to be a common problem (simply search or ask StackOverflow). I also needed to go the other direction and turn things like "1.5 MB" into 1,572,864. Here's what I cooked up in C99:

// Adapted from http://stackoverflow.com/questions/3758606/
// how-to-convert-byte-size-into-human-readable-format-in-java
void to_human_readable_byte_count(long bytes,
                                  int si,
                                  double *coeff,
                                  const char **units)
{
    // Static lookup table of byte-based SI units
    static const char *suffix[][2] = { { "B",  "B"   },
                                       { "kB", "KiB" },
                                       { "MB", "MiB" },
                                       { "GB", "GiB" },
                                       { "TB", "TiB" },
                                       { "EB", "EiB" },
                                       { "ZB", "ZiB" },
                                       { "YB", "YiB" } };
    int unit = si ? 1000 : 1024;
    int exp = 0;
    if (bytes > 0) {
        exp = min( (int) (log(bytes) / log(unit)),
                   (int) sizeof(suffix) / sizeof(suffix[0]) - 1);
    }
    *coeff = bytes / pow(unit, exp);
    *units  = suffix[exp][!!si];
}

// Convert strings like the following into byte counts
//    5MB, 5 MB, 5M, 3.7GB, 123b, 456kB
// with some amount of forgiveness baked into the parsing.
long from_human_readable_byte_count(const char *str)
{
    // Parse leading numeric factor
    char *endptr;
    errno = 0;
    const double coeff = strtod(str, &endptr);
    if (errno) return -1;

    // Skip any intermediate white space
    while (isspace(*endptr)) ++endptr;

    // Read off first character which should be an SI prefix
    int exp  = 0;
    int unit = 1024;
    switch (toupper(*endptr)) {
        case 'B':  exp =  0; break;
        case 'K':  exp =  3; break;
        case 'M':  exp =  6; break;
        case 'G':  exp =  9; break;
        case 'T':  exp = 12; break;
        case 'E':  exp = 15; break;
        case 'Z':  exp = 18; break;
        case 'Y':  exp = 21; break;

        case ' ':
        case '\t':
        case '\0': exp =  0; goto done;

        default:   return -1;
    }
    ++endptr;

    // If an 'i' or 'I' is present use SI factor-of-1000 units
    if (toupper(*endptr) == 'I') {
        ++endptr;
        unit = 1000;
    }

    // Next character must be one of B/empty/whitespace
    switch (toupper(*endptr)) {
        case 'B':
        case ' ':
        case '\t': ++endptr;  break;

        case '\0': goto done;

        default:   return -1;
    }

    // Skip any remaining white space
    while (isspace(*endptr)) ++endptr;

    // Parse error on anything but a null terminator
    if (*endptr) return -1;

done:
    return exp ? coeff * pow(unit, exp / 3) : coeff;
}

Using Argp with MPI-based applications

2011-01-15T23:29:00.000Z

Argp is a great, great parser for command line options. When using it in MPI-based applications, there's a catch in that you want only one MPI rank to print --help information, usage warnings, error messages, etc. Otherwise, you get a whole mess of repeated, jumbled output as each MPI rank squawks about the same problem.

Here's a wrapper for Argp for use in MPI-based applications that solves just this nuisance:

#include <stdio.h>
#include <unistd.h>
#include "argp.h"

/**
 * Call <a href="http://www.gnu.org/s/libc/manual/html_node/Argp.html"
 * >Argp</a>'s \c argp_parse in an MPI-friendly way.  Processes
 * with nonzero rank will have their \c stdout and \c stderr redirected
 * to <tt>/dev/null</tt> during \c argp_parse.
 *
 * @param rank MPI rank of this process.  Output from \c argp_parse
 *             will only be observable from rank zero.
 * @param argp      Per \c argp_parse semantics.
 * @param argc      Per \c argp_parse semantics.
 * @param argv      Per \c argp_parse semantics.
 * @param flags     Per \c argp_parse semantics.
 * @param arg_index Per \c argp_parse semantics.
 * @param input     Per \c argp_parse semantics.
 *
 * @return Per \c argp_parse semantics.
 */
error_t mpi_argp_parse(const int rank,
                       const struct argp *argp,
                       int argc,
                       char **argv,
                       unsigned flags,
                       int *arg_index,
                       void *input);

error_t mpi_argp_parse(const int rank,
                       const struct argp *argp,
                       int argc,
                       char **argv,
                       unsigned flags,
                       int *arg_index,
                       void *input)
{
    // Flush stdout, stderr
    if (fflush(stdout))
        perror("mpi_argp_parse error flushing stdout prior to redirect");
    if (fflush(stderr))
        perror("mpi_argp_parse error flushing stderr prior to redirect");

    // Save stdout, stderr so we may restore them later
    int stdout_copy, stderr_copy;
    if ((stdout_copy = dup(fileno(stdout))) < 0)
        perror("mpi_argp_parse error duplicating stdout");
    if ((stderr_copy = dup(fileno(stderr))) < 0)
        perror("mpi_argp_parse error duplicating stderr");

    // On non-root processes redirect stdout, stderr to /dev/null
    if (rank) {
        if (!freopen("/dev/null", "a", stdout))
            perror("mpi_argp_parse error redirecting stdout");
        if (!freopen("/dev/null", "a", stderr))
            perror("mpi_argp_parse error redirecting stderr");
    }

    // Invoke argp per http://www.gnu.org/s/libc/manual/html_node/Argp.html
    error_t retval = argp_parse(argp, argc, argv, flags, arg_index, input);

    // Flush stdout, stderr again
    if (fflush(stdout))
        perror("mpi_argp_parse error flushing stdout after redirect");
    if (fflush(stderr))
        perror("mpi_argp_parse error flushing stderr after redirect");

    // Restore stdout, stderr
    if (dup2(stdout_copy, fileno(stdout)) < 0)
        perror("mpi_argp_parse error reopening stdout");
    if (dup2(stderr_copy, fileno(stderr)) < 0)
        perror("mpi_argp_parse error reopening stderr");

    // Close saved versions of stdout, stderr
    if (close(stdout_copy))
        perror("mpi_argp_parse error closing stdout_copy");
    if (close(stderr_copy))
        perror("mpi_argp_parse error closing stderr_copy");

    // Clear any errors that may have occurred on stdout, stderr
    clearerr(stdout);
    clearerr(stderr);

    // Return what argp_parse returned
    return retval;
}

C header-only unit testing with FCTX

2011-01-04T18:48:00.000Z

Just a quick hat tip to FCTX, a library I've found invaluable these past few months. FCTX provides header-only unit testing for C. Sure, if you're in C++ land there's a ton of xUnit-like frameworks available (with Boost.Test being my favorite), but for vanilla C projects FCTX wins hands down.

As an example, here's something I put together for a Stack Overflow response:

#include <ctype.h>

void trim(char *a)
{
    char *b = a;
    while (isspace(*b))   ++b;
    while (*b)            *a++ = *b++;
    *a = '\0';
    while (isspace(*--a)) *a = '\0';
}

/* See http://fctx.wildbearsoftware.com/ */
#include "fct.h"

FCT_BGN()
{
    FCT_QTEST_BGN(trim)
    {
        { char s[] = "";      trim(s); fct_chk_eq_str("",    s); } // Trivial
        { char s[] = "   ";   trim(s); fct_chk_eq_str("",    s); } // Trivial
        { char s[] = "\t";    trim(s); fct_chk_eq_str("",    s); } // Trivial
        { char s[] = "a";     trim(s); fct_chk_eq_str("a",   s); } // NOP
        { char s[] = "abc";   trim(s); fct_chk_eq_str("abc", s); } // NOP
        { char s[] = "  a";   trim(s); fct_chk_eq_str("a",   s); } // Leading
        { char s[] = "  a c"; trim(s); fct_chk_eq_str("a c", s); } // Leading
        { char s[] = "a  ";   trim(s); fct_chk_eq_str("a",   s); } // Trailing
        { char s[] = "a c  "; trim(s); fct_chk_eq_str("a c", s); } // Trailing
        { char s[] = " a ";   trim(s); fct_chk_eq_str("a",   s); } // Both
        { char s[] = " a c "; trim(s); fct_chk_eq_str("a c", s); } // Both
    }
    FCT_QTEST_END();
}
FCT_END();

The logic isn't rocket science, of course. But testing it in C without resorting to external libraries and complicated makefiles shouldn't be rocket science either. Provided that fct.h is in the same directory, this source will compile and run.

getopt.h:195: error: redefinition of 'struct option'

2010-12-23T14:00:00.002Z

While working with Gnulib recently I kept running into problems with its getopt-gnu module:

In file included from ../../esio/gnulib/argp.h:24:0,
                 from ../../esio/apps/esio_bench.c:7:
../gnulib/getopt.h:195:8: error: redefinition of ‘struct option’
/usr/include/getopt.h:106:8: note: originally defined here
../gnulib/getopt.h:242:12: error: conflicting types for ‘getopt_long’
/usr/include/getopt.h:175:12: note: previous declaration of ‘getopt_long’ was here
../gnulib/getopt.h:246:12: error: conflicting types for ‘getopt_long_only’
/usr/include/getopt.h:179:12: note: previous declaration of ‘getopt_long_only’ was here
make[2]: *** [esio_bench.o] Error 1

which looked suspiciously like someone else's problem but couldn't be resolved by a simple make distclean.

Turns out that I failed to include the required

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

at the top of my source file. Oops. This detail is obvious in hindsight but slightly difficult to dig out of the relevant section of the gnulib manual.

snprintf + realloc = snprintf_realloc

2010-12-16T21:22:00.002Z

Recently, in C99, I needed to format some arbitrary-length string content using snprintf. I could re-use the same buffer because I only needed one result at a time. I wanted to (hopefully) minimize the number of malloc calls necessary. Finally, I wanted to (mostly) maintain snprintf's interface including its error semantics.

My use case looked like this...

    char *buf = NULL;
    size_t len = 0;
    while (/* more work */) {
        if (0 > snprintf_realloc(&buf, &len, /*format specifier*/, ...)) {
            /* error occurred, report it, break */
        }
        /* use buf to perform work */
    }
    if (buf) free(buf);

Here is the snprintf_realloc method I created including the relevant Doxygen and FCTX-based unit tests. The source below should build and run provided that fct.h is also in the same directory:

#include <stdarg.h>
#include <stdlib.h>

#include "fct.h"

/**
 * Call snprintf(*str, *size, format, ...) and reallocate the buffer
 * pointed to by *str as appropriate to contain the entire result.  On
 * exit, *str and *size will contain a pointer to the
 * realloced buffer and its maximum usable size, respectively.
 *
 * The reallocation scheme attempts to reduce the reallocation calls when the
 * same str and size arguments are used repeatedly.  It is
 * valid to pass *str == NULL and *size == 0 and then have
 * the buffer allocated to perfectly fit the result.
 *
 * @param[in,out] str    Pointer to the buffer in which to write the result.
 * @param[in,out] size   Pointer to the initial buffer size.
 * @param[in]     format Format specifier to use in sprintf call.
 * @param[in]     ...    Variable number of arguments corresponding
 *                       to \c format.
 *
 * @return On success, the number of characters (not including the trailing
 *         '\0') written to *str.  On error, a negative value
 *         is returned, *str is freed and *size
 *         is set to zero.
 */
int snprintf_realloc(char **str, size_t *size, const char *format, ...);

int
snprintf_realloc(char **str, size_t *size, const char *format, ...)
{
    int retval, needed;
    va_list ap;
    va_start(ap, format);
    while (   0     <= (retval = vsnprintf(*str, *size, format, ap)) // Error?
           && *size <  (needed = retval + 1)) {                      // Space?
        va_end(ap);
        *size *= 2;                                                  // Space?
        if (*size < needed) *size = needed;                          // Space!
        char *p = realloc(*str, *size);                              // Alloc
        if (p) {
            *str = p;
        } else {
            free(*str);
            *str  = NULL;
            *size = 0;
            return -1;
        }
        va_start(ap, format);
    }
    va_end(ap);
    return retval;
}

FCT_BGN()
{
    FCT_FIXTURE_SUITE_BGN(snprintf_realloc)
    {
        const char s[] = "1234";
        char *ptr;
        size_t size;

        FCT_SETUP_BGN()
        {
            ptr  = NULL;
            size = 0;
        }
        FCT_SETUP_END();

        FCT_TEARDOWN_BGN()
        {
            if (ptr) free(ptr);
        }
        FCT_TEARDOWN_END();

        FCT_TEST_BGN(snprintf_realloc)
        {
            // Initial should cause malloc-like behavior
            fct_chk_eq_int(4, snprintf_realloc(&ptr, &size, "%s", s));
            fct_chk(ptr);
            fct_chk_eq_int(size, 5);
            fct_chk_eq_str(ptr, s);

            // Repeated size should not cause any new buffer allocation
            {
                char *last_ptr = ptr;
                fct_chk_eq_int(4, snprintf_realloc(&ptr, &size, "%s", s));
                fct_chk(last_ptr == ptr);
                fct_chk_eq_int(size, 5);
                fct_chk_eq_str(ptr, s);
            }

            // Request requiring more than twice the space should
            // realloc memory to fit exactly.
            fct_chk_eq_int(12, snprintf_realloc(&ptr, &size, "%s%s%s",
                                                s, s, s));
            fct_chk_eq_int(size, 13);
            fct_chk_eq_str(ptr, "123412341234");

            // Request requiring less than twice the space should
            // cause a doubling of the buffer size.
            fct_chk_eq_int(16, snprintf_realloc(&ptr, &size, "%s%s%s%s",
                                                s, s, s, s));
            fct_chk_eq_int(size, 26);
            fct_chk_eq_str(ptr, "1234123412341234");
        }
        FCT_TEST_END();
    }
    FCT_FIXTURE_SUITE_END();
}
FCT_END()

Valgrind, after suppressing some unrelated FCTX warnings, gives this test a clean bill of health. If anyone's interested, aside from passing NULL pointers in or relying upon undefined snprintf behavior, I'd love to hear of ways to break this particular snippet.

APS DFD 2010 Talk Slides

2010-12-04T01:58:00.000Z

As promised, sort of, here are the slides from my APS DFD 2010 talk:

A WENO-Based Code for Investigating RANS Model Closures for Multicomponent Hydrodynamic Instabilities

Go listen to The Mikie Show

2010-12-04T01:54:00.000Z

An old acquaintance has been running an amazing radio show called The Mikie Show for some time now. If you get the chance, take a listen. Think Mister Rogers on the radio if Mister Rogers happened to be a sound designer/musician with a groanful sense of humor.

Recently my friend Billy, an arborist, was in episode 25. Billy's fun. When asked "Does hugging a tree really help it?" Billy didn't miss a beat and replied "It depends what you're going to do next."

Imminent (but probably dizzy) bundle o' joy

2010-11-18T04:16:00.001Z

My wife amazes me most days. Today especially.

The bump's our son. He should be arriving in the next three to seven weeks.

GNU Parallel

2010-11-12T05:16:00.000Z

I discovered GNU Parallel today. Phenomenal. Especially for simple post processing on multicore systems. Watch the video-- it is quite good. Here are three sample commands taken from a bash script where $1 is set to be a particular prefix...

parallel --eta -j +0 ./gnuplot.sh -g -F {} $1.r???? ::: `seq 1 8`

Runs a custom gnuplot-based shell script eight times against a file glob. Each time provide a different argument in the set {1,2,...,7,8}. The results are stored as $1.f1.gif, $1.f2.gif, etc.

ls -1 $1.r???? | parallel --eta -j +0 'montage -geometry 1280x960 -tile 4x2 {}.f?.gif {}.miff'

Create a tiled montage of the eight output files mentioned above and store the result in $1.miff.

ls -1 $1.r???? | parallel --eta -j +0 'convert {}.miff {}.gif'

Convert the resulting MIFFs into GIFs.

I win. On my system, 8 - 0.

Fortran generic interface mixed with iso_c_binding

2010-10-15T04:58:00.001Z

I'm going crazy trying to figure out some Fortran generic interface funkiness. The compiler's complaining that a function call doesn't match a generic interface despite the fact that only one module procedure statement is in the interface block and the same damned code compiles when I write the concrete procedure name instead of the generic name.

Just so I've got it later, this happens to work:

module testing

  use, intrinsic :: iso_c_binding, only: c_float, c_double
  private :: c_float, c_double

  interface test
    module procedure test1
    module procedure test2
  end interface

contains

  real(c_float) function test1 (x)
    real(c_float), intent(in) :: x(*)
    test1 = 2 * x(1)
  end function test1

  real(c_double) function test2 (x)
    real(c_double), intent(in) :: x(*)
    test2 = 2 * x(1)
  end function test2

end module testing

program main
  use testing
  implicit none
  real(4) :: s(1)
  real(8) :: d(1)
  s = 2.0
  d = 2d0
  write (*,*) test(s)
  write (*,*) test(d)
end program main

This snippet, however, is not the thing that's breaking...

APS DFD 2010 abstract

2010-09-17T18:52:00.000Z

From the APS DFD 2010 schedule:

Session AW: Instability: Richtmyer-Meshkov/Rayleigh-Taylor I

8:00 AM–10:10 AM, Sunday, November 21, 2010
Hyatt Regency Long Beach Room: Regency C

Chair: Arindam Banerjee, Missouri University of Science & Technology

Abstract: AW.00010 : Comparisons of a Reynolds-Averaged Navier--Stokes Model with Self-Similar Solutions for Large Atwood Number Rayleigh--Taylor Mixing

9:57 AM–10:10 AM

Preview Abstract

Authors:

Rhys Ulerich
(University of Texas at Austin)

Oleg Schilling
(Lawrence Livermore National Laboratory

A new high-order, multicomponent, weighted essentially nonoscillatory (WENO) implementation of a three- and four- equation Reynolds-averaged Navier-Stokes (RANS) model incorporating both mechanical and scalar turbulence is used to simulate intermediate-to-large Atwood number Rayleigh-Taylor turbulent mixing. The predicted RANS mixing layer evolution is compared with the analytical self-similar solutions of the transport equations. The terms in the transport equation budgets are compared in detail to their self-similar profiles across the mixing layer. Additionally, the sensitivity of the RANS solutions to variations in the initial conditions and in the model coefficients is explored. The implications of these results for advanced modeling of Rayleigh-Taylor turbulent mixing are discussed.

internal error: assertion failed: lower_expr: bad kind

2010-09-14T18:44:00.001Z

I'm back in Texas.

In case it helps anyone else, sometimes I see errors like the following

 ../../sz/suzerain/diffwave.hpp(65): internal error: assertion failed: lower_expr: bad kind (shared/edgcpfe/lower_il.c, line 14342)

coming from Intel's icpc 10.1 20081024. This usually happens when I use reinterpret_cast while calling a function but screw up the const-correctness of the type to which I'm casting.

I've never been able to get a clean recreate to provide to Intel Support.

Poster from this summer's internship

2010-08-12T02:55:00.002Z

As part of my summer internship at Lawrence Livermore National Laboratory I'll be presenting the following poster tomorrow:

Towards a WENO-based code for investigating RANS model closures for hydrodynamic instabilities

di·ver·tisse·ment

2010-08-09T06:00:00.000Z

View Mines Rd to Big Basin in a larger map

Friday Musings

2010-07-30T16:08:00.001Z

First, Mumford & Son's "Little Lion Man" is amazing. Listen. Purchase. Bonus points for accidentally starting two or more copies and noticing it sounds excellent in a round. Found it on thesixtyone.com, which is a great source for day-to-day working music and occasionally emits gems like this one.

Second, the Richtmyer-Meshkov setup I've been mucking with for two days is working. In it, a shock wave interacts with a perturbed, constant pressure interface between two different densities (details). Courtesy of gnuplot and gifsicle, here's an animated gif showing density evolution for a problem that's periodic in y and reflective at both x boundaries (click to view animation, 3.2M):

Just a roll up?

The reflective condition at the left boundary causes a re-reshock, which isn't usually what folks are interested in for these problems. I'm not entirely sure how to enforce a strict inflow condition for compressible simulations like these-- presumably such an inflow condition would eat the outgoing characteristics and stop the re-reshock seen here.

Glam Shot

2010-07-23T04:04:00.009Z

Today at LLNL I was involved in some photos for a story on the summer students at ISCR:

Daniel from Ghana, Hilary from Dallas, and me

The screen images are from the 2D WENO-based code I've been revising to investigate RANS modeling of Rayleigh-Taylor instabilities with Dr. Oleg Schilling. As a bit of nerd sniping, I made sure Trac, VIM's NERD tree, and some Mathematica reproducing classic results from Roe's JCP 1981 paper were clearly visible.

Unrelated, today I also had a couple of beers with Ondřej Čertík of sympy and theoretical physics fame. Unbelievably swell guy— easily one of the smartest one or two people with whom I've ever shared a brew.

Update: Entertainingly, in the related story the lab released, the caption shows me as a "Livermore scientist".

Burst traffic smoothing for SIP processing elements

2010-07-20T01:43:00.001Z

I received another patent grant from my IBM days. This time its United States Patent 7,742,417:

Mechanisms for burst traffic smoothing for Session Initiation Protocol (SIP) processing elements are provided. A dispatch queue management engine determines whether a received packet is a TCP or UDP packet. If the packet is a TCP packet, the packet is automatically added to the dispatch queue. If the packet is a UDP packet, a value for a drop function f is generated and a random or pseudo-random number r is generated. If r has a predetermined relationship to f, then the UDP packet is added to the dispatch queue, otherwise the UDP packet is discarded. The value for f is based on the current dispatch queue load, the network quality, the retransmission rate, and the allowable drop rate. Thus, the determination as to whether to drop UDP packets or not is configurable by an administrator and also adaptable to the current network and dispatch queue conditions.

The idea for this patent came from dealing with overflowing UDP buffers during a SIP-enabled Java application server stress test. We observed that garbage collection events would pause the application server and cause a backlog of SIP retransmissions over UDP to grow. The invention is a way to selectively toss out UDP packets during such a retransmission burst in a way that would allow the application server to catch up again while maintaining a probabilistically high QoS.

Wireless tunes on the bike

2010-07-15T05:45:00.000Z

I'm quite stoked about combining three new toys:

BlueAnt's F4 Interphone Bluetooth helmet communication system
Sony's TMR-BT8IP Bluetooth iPod transmitter
HJC's IS-16 helmet

I've tried out the first two (to Kids of the K Hole no less, an excellent escapist riding anthem) and am supremely impressed. It'll be nice to not be tethered to the tankbag by a huge, flaky PS/2 cable a la my klunkotronic Chatterbox GMRS X1.

I can't speak to mounting the Interphone's helmet speakers since the new HJC hasn't arrived yet, but the mounting hardware BlueAnt provided looks solid. Also no verdict yet on the intercom functionality, but the full duplex sound should be a welcome improvement compared to the GMRS X1 (unless Pauly abuses said duplex connection by singing along to decidedly non-escapist anthems).

How to parse Fortran command line arguments using NAMELIST

2010-06-25T06:34:00.002Z

While searching around for a pure Fortran 90+ way to parse command line arguments and/or input files, I stumbled across Exploring NAMELIST by John S. Urban. John's article shows a way to use Fortran's NAMELIST feature to quickly parse arbitrary command line arguments while providing sane default behavior.

Like usual, I decided to overdevelop a simple idea. The complete sample is up on github. Assuming you've got an interface declaration in scope for get_command_arguments, here's what you'd see as an end consumer:

PROGRAM namelist_cli

  IMPLICIT NONE

  CHARACTER(LEN=255) :: string, message
  INTEGER :: status

  ! Declare and initialize a NAMELIST for values of interest
  INTEGER  :: i=1,   j=2
  REAL     :: s=3.3, t=4.4
  NAMELIST /cmd/ i, j, s, t

  ! Get command line arguments as a string ready for NAMELIST
  CALL get_command_arguments(string, "&cmd", "/")
  READ (string, NML=cmd, IOSTAT=status, IOMSG=message)
  IF (status /= 0) THEN
    WRITE (*,*) 'Error:   ', status
    WRITE (*,*) 'Message: ', message
    CALL EXIT (status)
  END IF

  ! Echo the results
  WRITE (*, NML=cmd)

END PROGRAM namelist_cli

With a little mucking around, the same NAMELIST could be used to read an input file from disk, allow the user to override individual options on the command line (say, for parametric studies), and also to dump the full setup into an output file for tracking purposes.

Summer internship

2010-05-31T20:11:00.010Z

Following a 2,500 mile motorcycle trip from Austin to the Bay Area with my friend Paul, I'll be starting a summer internship at Lawrence Livermore National Labs tomorrow morning.

Instead of my usual my thesis work, I'll be looking at RANS turbulence models/numerics for Rayleigh-Taylor problems. As preparation, over the past couple of months, I've been spending some time staring at conservation law numerics classics like LeVeque's book as well as lecture notes from Caltech's AE 232 course. On the odd chance it is useful to anyone, I've put my toy code up at http://github.com/RhysU/codefour.

The Fortran 90 code solves inviscid or viscid 1D scalar conservation laws of the form u_t + f(u)_x = 0 or u_t + f(u)_x = \nu u_{xx} using WENO reconstruction and a Lax-Friedrichs flux.

It is based heavily on one provided through Caltech's AE 232 course. Compared to the course's code, the following changes have been made:

Adding 3rd-order WENO reconstruction
Adding 2nd and 4th order viscous capabilities
Adding HDF5-based result output
Adding Doxygen-based documentation
Adding Mathematica-based error analysis against Hopf-Cole-based analytic solutions of the viscid Burger's equation

Made GSL Committer

2010-05-06T14:52:00.001Z

I was just added as a committer for the GNU Scientific Library (GSL) over on Savannah.

New Toy

2010-04-24T21:41:00.000Z

Courtesy of the employer, my new Lenovo W510 came in on Thursday.

make -j8 on a laptop is fun:

Gram's Gob Recipe

2010-03-04T05:11:00.001Z

My grandmother recently mailed me her gob recipe. Hands down, my Gram's are better than Yost's, your neighbor's, your mother's, or even your grandmother's. As a public service, I'm sharing this jewel.

Gram's Gobs

Cookies

Cream together

2 cups sugar,
½ cup crisco,
2 eggs, and
1 cup sour cream.

Then add 1 cup boiling water. Sift and add to the above

4 cups flour,
2 teaspoons baking soda,
½ teaspoon baking powder,
½ teaspoon salt, and
½ cup cocoa.

Drop by tablespoons onto an ungreased baking sheet. Bake at 450° Fahrenheit for exactly 5 minutes. Remove from the sheet at once.

Filling

First, mix together

5 tablespoons flour and
1 cup milk.

Cook the mixture until very thick and then let it cool. Separately, use a mixer to mix until light and fluffy

1 cup crisco or oleo,
1 cup powdered sugar,
¼ teaspoon salt, and
1 teaspoon vanilla.

Add the first mixture to the second mixture and mix well.

Assembly

Spread the filling between two cookies wrapping each one separately. Gram adds...

MPI Bandwidth Benchmarks

2010-02-17T23:04:00.000Z

I need to see if the wall times reported by mpiP for the MPI calls underneath FFTW 3.3alpha1's parallel transposes make sense. The TACC folks were excellent, as usual, and pointed me to the MVAPICH benchmarks.

One quick realization of the bandwidth test on lonestar and ranger shows the following:

Unfortunately, this information does not explain my bottleneck...

Ahhh... Corporate Governance

2010-01-28T03:38:00.000Z

Really?

We have agreed that if our stockholders fail to adopt the Amendment at the Special Meeting or any postponement or adjournment thereof, we will continue to seek to obtain the requisite approval at least as frequently as every six months until such stockholder adoption has been obtained.

Really. Thankfully, they caught and removed the draft paragraph's final "you stupid asshats".

Admittedly, I have not yet fully reviewed all the particulars. The choice paragraph appears at the bottom of page 6.

Using FFTW for in-place matrix transposition

2010-01-20T03:23:00.008Z

The Wikipedia article on in-place matrix transposition tantalizingly hints that FFTW can compute rectangular in-place matrix transposes efficiently (or as efficiently as one can). FFTW's manual, however, is silent. At least as far as section headings, function names, and the index are concerned.

Frigo and Johnson's paper, on page 6, spells out how to tap into the in-place transpose routines: create a rank-0 transform of vector rank 2 using transposed vector strides. I had to bug Matteo Frigo twice (nice guy!) to figure out how to accomplish this using the guru interface.

A C99 solution capable of handling the double-valued either row major or column major case looks like:

#include <assert.h>
#include <ctype.h>
#include <fftw3.h>

fftw_plan plan_transpose(char storage_type,
                         int rows,
                         int cols,
                         double *in,
                         double *out)
{
    const unsigned flags = FFTW_ESTIMATE; /* Do not destroy input */

    fftw_iodim howmany_dims[2];
    switch (toupper(storage_type)) {
        case 'R':
            howmany_dims[0].n  = rows;
            howmany_dims[0].is = cols;
            howmany_dims[0].os = 1;
            howmany_dims[1].n  = cols;
            howmany_dims[1].is = 1;
            howmany_dims[1].os = rows;
            break;
        case 'C':
            howmany_dims[0].n  = rows;
            howmany_dims[0].is = 1;
            howmany_dims[0].os = cols;
            howmany_dims[1].n  = cols;
            howmany_dims[1].is = rows;
            howmany_dims[1].os = 1;
            break;
        default:
            return NULL;
    }
    const int howmany_rank = sizeof(howmany_dims)/sizeof(howmany_dims[0]);

    return fftw_plan_guru_r2r(/*rank*/0, /*dims*/NULL,
                              howmany_rank, howmany_dims,
                              in, out, /*kind*/NULL, flags);
}

You can get a better performing plan if you're willing to use flags = FFTW_MEASURE or its big brothers. I chose to use FFTW_ESTIMATE in the sample above because I didn't want the planning process to destroy the input array.

You use the solution like:

void somelogic()
{
    char storage_type = /* One of 'R', 'r', 'C', 'c' */
    int rows          = /* Number of rows */
    int cols          = /* Number of columns */
    double *in        = /* ... */
    double *out       = /* ... */

    /* Plan the transpose once; transpose is in-place iff in == out */
    fftw_plan transpose = plan_transpose(storage_type, rows, cols, in, out);
    assert(transpose);

    /* Execute the plan potentially many times */
    fftw_execute(transpose);

    /* FFTW New-array Execute functions should be callable, too */
    /* Beware of mixing in-place and out-of-place planning and usage */
    double *another_in  = /* ... */
    double *another_out = /* ... */
    fftw_execute_r2r(transpose, another_in, another_out);

    /* Destroy the plan when completely done */
    fftw_destroy_plan(transpose);
}

There are many restrictions around the the New-array execute interface, so be careful if you go that route.

One very cool thing is that you don't need to know a priori that you're computing an in-place transpose. The planner works for out-of-place transposes as well. Another cool thing is that the plans are useful for rectangular transposes and square transposes alike with more efficient algorithms being used for the latter case.

I wrote myself a toy to test it out. Transposing a Row major, 3-by-2 matrix In-place gives:

[499 rhys@mentes misc]$ ./fftw_transpose I R 3 2

viewing extents [ 3 2 ] with strides [ 2 1 ] {
   0    1 
 100  101 
 200  201 
}
stored contiguously as {
   0    1  100  101  200  201 
}

Transposing in-place with 
(rdft-transpose-cut-3x2
  (rdft-rank0-ip-sq-tiledbuf/1-x2-x2)
  (rdft-rank0-iter-ci/1-x2))

viewing extents [ 2 3 ] with strides [ 3 1 ] {
   0  100  200 
   1  101  201 
}
stored contiguously as {
   0  100  200    1  101  201 
}

Transposing a Column major, 2-by-3 matrix Out-of-place gives:

[500 rhys@mentes misc]$ ./fftw_transpose O C 2 3

viewing extents [ 2 3 ] with strides [ 1 2 ] {
   0    1    2 
 100  101  102 
}
stored contiguously as {
   0  100    1  101    2  102 
}

Transposing out-of-place with 
(rdft-rank0-tiledbuf/1-x3-x2)

viewing extents [ 3 2 ] with strides [ 1 3 ] {
   0  100 
   1  101 
   2  102 
}
stored contiguously as {
   0    1    2  100  101  102 
}

From burning lots of time playing with FFTW's transpose capabilities, it seems that the planner returns NULL (i.e. bombs) whenever you attempt to transpose a non-contiguous matrix. That rules out specifying leading dimensions and/or arbitrary strides as far as I can tell. Someone please correct me if I'm wrong. I'd love to be wrong so I can use FFTW to emulate Intel MKL's mkl_trans.h functionality for systems where a newer version of MKL is not available.

Mathematica's FourierParameters necessary to emulate FFTW_FORWARD and FFTW_BACKWARD

2010-01-13T19:45:00.003Z

These are the FourierParameters necessary to have Mathematica's Fourier and InverseFourier discrete transforms behave like FFTW's forward and backward FFTs:

FFTWForward[expr_] := Fourier[expr, FourierParameters -> {1, -1}];
FFTWBackward[expr_] := InverseFourier[expr, FourierParameters -> {-1, -1}];

In particular, these functions/parameter choices do not scale the result when transforming in either direction.

Similar parameters can be used when performing analytic transforms using FourierTransform and InverseFourierTransform:

FFTWForwardAnalytic[expr_, t_, w_] := FourierTransform[expr, t, w, FourierParameters -> {1, -1}];
FFTWBackwardAnalytic[expr_, w_, t_] := InverseFourierTransform[expr, w, t, FourierParameters -> {-1, -1}];

Providing fill and for_each algorithms for Boost.MultiArray

2010-01-13T06:10:00.003Z

I sometimes want to fill a Boost.MultiArray with a particular value. For concrete implementations like boost::multi_array or boost::multi_array_ref, the task isn't too hard and std::fill works just fine. But getting an arbitrary-dimensional version working that only uses the MultiArray concept takes a little thought. This is important when working with views. I additionally want the solution to take advantage of the contiguous storage guarantees made by multi_array and multi_array_ref.

Because fill looks like for_each with an assignment-like UnaryFunction, I decided to write a MultiArray-specific for_each implementation. First comes a helper functor recursively templated on the dimensionality:

namespace { // anonymous

template<std::size_t D>
struct for_each_functor {

    BOOST_STATIC_ASSERT(D != 0); // Nonsensical behavior for zero dimensions
    BOOST_STATIC_ASSERT(D > 1);  // Ensure not instantiated for specialized values

    // See http://groups.google.com/group/boost-list/browse_thread/thread/e16f32c4411dea08
    // for details about why MultiArray::iterator::reference is used below.
    template<class MultiArray,
             class UnaryFunction >
    void operator()(MultiArray &x, UnaryFunction &f) const {
        for_each_functor<D-1> functor;
        for (typename MultiArray::iterator i = x.begin(); i != x.end(); ++i) {
            typename MultiArray::iterator::reference ri = *i;
            functor(ri,f);
        }
    }
};

template<> struct for_each_functor<1> {
    template<class MultiArray,
             class UnaryFunction >
    void operator()(MultiArray &x, UnaryFunction &f) const {
        std::for_each(x.begin(), x.end(), f);
    }
};

} // namespace anonymous

The functor uses std::for_each for the one dimensional base case. The BOOST_STATIC_ASSERTs are there just to keep me from shooting myself in the foot. A smarter version might specialize for the two- and three-dimensional cases, check the underlying strides, and then be certain to iterate in the faster directions first.

Now that we've got for_each_functor that works on the general MultiArray concept, it's straightforward to implement for_each:

/**
* Invoke \c f for each element in MultiArray \c x.  The
* order in which the invocations take place is undefined.
*
* @param x MultiArray over which to iterate.
* @param f UnaryFunction to invoke on each element of \c x.
*
* @return \c f.
*
* @see SGI's <a href="http://www.sgi.com/tech/stl/UnaryFunction.html">
*      UnaryFunction</a> concept for more information.
*/
template<class MultiArray,
         class UnaryFunction>
inline
UnaryFunction for_each(MultiArray &x,
                       UnaryFunction f) {
    for_each_functor<MultiArray::dimensionality>()(x,f);
    return f;
}

Please note that my for_each type signature is not identical with the STL algorithm's signature.

Since we know that multi_array and multi_array_ref have contiguous storage, we can specialize our version of for_each to linearly walk their underlying memory using the usual algorithms.

/**
* Invoke \c f for each element in <tt>boost::multi_array</tt> \c x.  The order
* in which the invocations take place is undefined.  This specialization
* takes advantage of <tt>boost::multi_array</tt>'s contiguous storage.
*
* @param x <tt>boost::multi_array</tt> over which to iterate.
* @param f UnaryFunction to invoke on each element of \c x.
*
* @return \c f.
*
* @see SGI's <a href="http://www.sgi.com/tech/stl/UnaryFunction.html">
*      UnaryFunction</a> concept for more information.
*/
template<class ValueType, std::size_t NumDims, class Allocator,
         class UnaryFunction>
inline
UnaryFunction for_each(boost::multi_array<ValueType,NumDims,Allocator> &x,
                       UnaryFunction f) {
    return std::for_each(x.data(), x.data() + x.num_elements(), f);
}

/**
* Invoke \c f for each element in <tt>boost::multi_array_ref</tt> \c x.  The
* order in which the invocations take place is undefined.  This specialization
* takes advantage of <tt>boost::multi_array_ref</tt>'s contiguous storage.
*
* @param x <tt>boost::multi_array_ref</tt> over which to iterate.
* @param f UnaryFunction to invoke on each element of \c x.
*
* @return \c f.
*
* @see SGI's <a href="http://www.sgi.com/tech/stl/UnaryFunction.html">
*      UnaryFunction</a> concept for more information.
*/
template<class ValueType, std::size_t NumDims,
         class UnaryFunction>
inline
UnaryFunction for_each(boost::multi_array_ref<ValueType,NumDims> &x,
                       UnaryFunction f)
{
    return std::for_each(x.data(), x.data() + x.num_elements(), f);
}

With for_each complete, we can move on to implementing fill by mapping an assignment-like functor over over a MultiArray.

In my particular application I have MultiArrays with complex-valued elements from a couple of different numeric libraries (e.g. std::complex<T>, T[2]) and I'd like to be able to fill a MultiArray with a value for which a member operator= may not be available. To accommodate this need, I created a templated assignment functor which can be specialized using schmancy boost::enable_if techniques based on both the target and source of the assignment:

/**
* A functor that performs assignment to type \c Target from type \c Source.
* \c Target must be assignable from \c Source.
*
* The \c Enable template parameter allows using <tt>boost::enable_if</tt> to
* specialize or partially specialize the functor per <tt>enable_if</tt>'s
* documentation section <a
* href="http://www.boost.org/doc/libs/1_41_0/libs/utility/enable_if.html">
* 3.1: Enabling template class specializations</a>.
*/
template<class Target, class Source, class Enable = void>
struct assign {

    /**
     * Create an instance which assigns \c s when applied.
     *
     * @param s source of assignment operation occurring via
     *          <tt>operator()</tt>.
     */
    assign(const Source &s) : s_(s) {}

    /**
     * Assign the value provided at construction to \c t.
     *
     * @param t to be assigned.
     */
    void operator()(Target& t) const { t = s_; }

private:
    const Source &s_; /**< Source for assignment operations */
};

The default implementation above makes the extra parameters seem overkill. I'll give an example where it is useful later.

Using the assignment functor and for_each, here's our one-line fill:

/**
* Fill MultiArray \c x with the value \c v.  MultiArray <tt>x</tt>'s elements
* must be assignable from \c v.  The underlying assignment uses
* ::assign to allow specializations of that functor to be found.
*
* @param x MultiArray to fill.
* @param v Value with which to fill \c x.
*/
template<class MultiArray, class V>
void fill(MultiArray &x, const V &v) {
    for_each(x, assign<typename MultiArray::element,V>(v));
}

Finally, to show the spurious-looking assign functor template parameters are useful, here's an assign functor for the complex-valued use case I mentioned above:

/**
* A specialization of the assign functor to properly handle the case where \c
* Target is a recognized complex type according to
* ::suzerain::complex::traits::is_complex.  It uses
* ::suzerain::complex::assign_complex to perform the assignment, and therefore
* supports all types that \c assign_complex does.
*/
template<class Target, class Source>
struct assign<
    Target,
    Source,
    typename boost::enable_if<
        ::suzerain::complex::traits::is_complex<Target>
    >::type >
{
    /**
     * Create an instance which assigns \c s when applied.
     *
     * @param s source of assignment operation occurring via
     *          <tt>operator()</tt>.
     */
    assign(const Source &s) : s_(s) {};

    /**
     * Assign the value provided at construction to \c t.
     *
     * @param t to be assigned.
     */
    void operator()(Target& t) const {
        ::suzerain::complex::assign_complex(t, s_);
    }

private:
    const Source &s_; /**< Source for assignment operations */
};

That very last code snippet is missing details and won't compile. Note that the suzerain::complex::assign_complex invocation may itself be a template and it receives appropriate type information. This turned out to be important when I wanted to fill a complex-valued MultiArray with a real-valued scalar. I needed two versions of assign_complex: one that took as a Source a recognized complex type and one that took a scalar type. boost::enable_if and boost::disable_if made it easy to provide an assign_complex template that did the right thing.

Given the right #includes in the right places, the other snippets should all compile. Please let me know if they do not.

This post topic came from a question I asked on boost-users a couple of days back. Hopefully someone else finds this post and can use this fill routine. Thank you to Ronald Garcia for his response and for maintaining Boost.MultiArray.

Small gotcha when combining Boost.Array, Assign, and Test

2009-12-31T19:40:00.001Z

When attempting to use Boost.Assign's list_of in expressions like BOOST_CHECK_EQUAL(somearray, list_of(x)(y)(z)), where somearray has a type like boost::array<T,N>, I kept running into funky errors that there was no operator<< defined with the appropriate types. So I'd define one at global scope, but the errors kept arising from the depths of BOOST_CHECK_EQUAL. It turns out that I needed to define such an operator<< in the boost namespace:

namespace boost {

template< typename charT, typename traits, typename T, ::std::size_t N >
::std::basic_ostream& operator<<(
        ::std::basic_ostream &os, const ::boost::array &array)
{
    os << '[' << N << "]{ ";
    ::std::copy(array.begin(),
                array.end(),
                ::std::ostream_iterator(os, " "));
    os << '}';
    return os;
}

}

Lit search

2009-12-24T06:45:00.000Z

I finally got around to doing some lit search towards both my dissertation proposal and a poster that I need to present in February. CiteULike's been a nice way to organize the information, especially combined with SyncUThink.

The following is a snapshot of the author's that I'll be wading through:

Orlu Abe Abedin Adams Adrian Adumitroaie Akbary Alazard Aldama Alfonsi Alfredsson Alving Ammer Andersson Andrieu Antonia Arnal Arnette Ashrafian Ashworth Baer Bailly Bake Bakken Balakrishnan Basili Bazilevs Beckie Beddini Behnia Benay Berggren Berland Bertolotti Berzins Betyaev Bhaganagar Bhattacharya Biagioni Bientinesi Bijl Biringen Blaisdell Blottner Bochev Boehm Bogey Boninger Bonnet Borodai Bosques Botella Bottaro Bourassa Bowersox Boyce Boyd Bradshaw Brandt Brasseur Brassuer Brethouwer Browne Brun Buckley Buter Cabrit Cal Calo Cambon Candler Cartwright Castillo Ceccio Ch Chae Chan Chanetz Chang Chauhan Cheban Chen Cheskidov Childs Chiyaka Choi Chong Chow Christon Christophorou Chung Clemens Coleman Collis Comte Constantin Constantinescu Cooper Corbett Cortissoz Creighton Cremonesi Crouch Cutbirth Da Prato Dai Das Davidson De Chant De Graaff Debussche Delery Demko Dexun Dilipkumar Dmitruk Dolling Donatelli Dong Donovan Donzis Doom Dotter Doucet Dowling Duan Dubos Dufresne Dumas Durbin Dussauge Eaton Ecer Egolf Elghobashi Elliott Elsinga Emerson Enrique Erturk Eswaran Ewing Fabignon Favre Fedorchenko Feijóo Feireisl Fernando Fernholz Ferrante Ferziger Fiala Fife Finleyt Fisher Fitzgerald Flandoli Floryan Fox Foysi Fransson Friedlander Friedrich Fu Fueyo Fujii Gad El Hak Gallagher Ganapathisubramani Garcia Garg Garnier Gatski George Gerolymos Geurts Ghosal Ghuloum Glimm Gokcol Goldstein Goyne Grant Grasso Gravemeier Greenside Grosse Guarini Gubinelli Guermond Gullbrand Guocan Gwilliam Haak Haberkorn Hadjadj Hale Hamba Haminh Hamman Haselbacher Hattori Hauke Heinrichs Heinz Helton Hemsch Henkes Henningson Hilburger Hillier Hochstein Holloway Honein Hou Hoyas Huang Hudea Hughes Hwang Imad Ishihara Ivanov Iwamoto Jackson Jacobs Jayaram Jiménez Jiménez Johansson Johnson Jordan Kachanov Kafoussias Kaneda Kasagi Kattelans Katz Kawaguchi Kawai Kawamura Kedia Kenny Khujadze Kiehn Kim Kimmel Kirby Klebanoff Kleiser Klewicki Kline Knight Kopriva Korimilli Kravchenko Krimmelbein Krishnamurty Krishnan Kritsuk Krogstad Kuerten Kunstmann Kwok Lan Langford Lechner Lee Lehoucq Lele Leonard Leonardo Levin Li Liepmann Lin Linan Liu Lo Luchini Luckring Luker Lumley Lumsdaine Lund Luo Lyons Ma Maeder Maekawa Maghami Mahesh Malaya Malik Mallinson Manca Manhart Mansour Martin Mashayek Maslowe Matsubara Matsuo Matthaeus Mattingly Mazzei McComb Mellet Meneguzzi Metzger Miao Miller Mizuno Moin Monkewitz Montgomery-Smith Morinishi Moser Mossi Moussaoui Mukund Muller Munz Murthy Müller Na Nagib Najjar Najm Nakabayashi Narasimha Neemann Newhall Nickels Nicoud Norman Nosetips Oberkampf Oberlack Oden Orszag Orti Osterlund Padoan Padula Paicu Pancaldo Park Parkman Parneix Passoni Patel Paull Pavlović Pego Pekurovsky Periaux Perlin Peterson Petrovan Piomelli Piro Pirozzoli Plesniak Pletcher Poggie Politano Pope Pouquet Prabhu Prudhomme Psihoyios Purtell Quadrio Quincy Röckner Radespiel Reed Reichert Reisenthel Rempfer Reshotko Reynolds Ribert Ringuette Ristorcelli Rizzetta Roach Robinson Rodi Rogallo Rogers Romito Roshko Roy Rubinstein Ruzmaikina Sadowski Sagaut Samimy Sandberg Sandham Sarkar Satofuka Scarano Schiffner Schilling Schlatter Schmitt Schneider Schochet Schwarz Scientist Scovazzi Seckel Sesterhenn Settles Shadid Shariff Sharp Sheikh Shen Shewchuk Shikazono Shvydkoy Shyy Sidi Simens Simos Sinha Skote Slepoy Smith Smith? Smits So Soliman Soltani Sommer Spalart Spina Spiteri Squires Stalker Steinkamp Stolz Stoynov Stroud Sumitani Suzuki Svehla Sverak Swartz Sénéchal Tamano Tamura Taulbee Tavoularis Taylor Tchuenche Teitel Tenaud Terracol Thangada Thomas Tomkins Torres-Nieves Trefethen Trouvé Trucano Tsuji Tu Tyhtila Ulerich Unalmis Urbin Vahala Vallet Van Zee Vandromme Vasseur Vedula Vegesna Venugopal Visbal Viswanath Vreman Wagner Waleffe Wallace Wang Watmuff Waymire Webster Wendroff Wesseling Williams Winkel Winter Wolf Wray Wu Xenos Xinliang Xu Ya Yaglom Yakhot Yan Yanwen Yao Yeoh Yeung Zabaras Zampaglione Zandonade Zang Zeman Zeytounian Zhang Zhong Zhongmin Zhou Zimmermann Zolver Zumbrun de Freitas de Lange lio Vargas van Oudheusden van de Geijn van der Heul

Undergraduate research...

2009-12-22T23:50:00.002Z

...the gift to yourself that just keeps on giving.

Blue Diamond Mine

2009-12-17T02:36:00.001Z

Weird linkage stemming from the shuffle tonight. "Blue Diamond Mine" by the Texas Rubies sounds insanely like "Where Did You Sleep Last Night?" by Nirvana. Including all the lyrical structure. This isn't too crazy considering the former's genre and the known folk origin of the later, but it is odd when from nowhere a song you've never heard before strikes you as this uncannily familiar. I leave finding an MP3 of "Blue Diamond Mine" as an exercise for the consumer. It appears on Down to the Promised Land.

Getting 32-bit libstdc++.so.5 in {Karmic Koala,Lucid Lynx} on a 64-bit system

2009-11-01T05:14:00.006Z

The Intel compilers require a 32-bit version of libstdc++.so.5 to function. Getting a copy of that build of that library requires a couple of steps in the latest 64-bit Ubuntu {9.10,10.04}. These have been alluded to in other places. Just in case someone needs the details, here they are:

Download the i386 libstdc++ package for Jaunty: http://packages.ubuntu.com/jaunty/i386/libstdc++5/download
Unpack the .deb archive using ar vx ~/Desktop/libstdc++5_3.3.6-17ubuntu1_i386.deb
Unpack the resulting data.tar.gz using tar xzvf data.tar.gz
Confirm that you got the 32-bit version using file usr/lib/libstdc++.so.5.0.7
Install the library into /usr/lib32 using sudo install usr/lib/libstdc++.so.5.0.7 /usr/lib32
Change to /usr/lib32 using cd /usr/lib32
Create a symlink to finish the task via sudo ln -s libstdc++.so.5.0.7 libstdc++.so.5

Many thanks to everyone who provided the source material:

Updated 24 April 2010: Appears that this procedure also works on the Lucid Lynx release candidate.

Sugar Plum Plug

2009-10-28T19:19:00.000Z

Milady got the Nutcracker casting she's talked about wanting for several years. Go hit up the show on one of December 6, 12, 18, 19 @ 7:30pm, 21, or 23 to see her grinning from ear to ear.

BOOST_UBLAS_SHALLOW_ARRAY_ADAPTOR

2009-10-26T22:57:00.000Z

uBLAS has a shallow_array_adaptor that allows you to wrap an existing block of data with a uBLAS vector. This is useful if you want to wrap raw data provided by someone else. Exposing this functionality is questionable in terms of encapsulation and ownership, and so it isn't available by default. Consequently the documentation is pretty scarce. Here's an example of it's use in the context of a Boost.Test test case.

#define BOOST_TEST_MODULE $Id: test_complex.cc 1817 2009-03-18 19:01:25Z rhys $

// Before ublas #include, enable boost::numeric::ublas::shallow_array_adaptor<T>
#define BOOST_UBLAS_SHALLOW_ARRAY_ADAPTOR 1

#include <algorithm>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/test/included/unit_test.hpp>
#include <boost/foreach.hpp>
#include <complex>

namespace ublas = boost::numeric::ublas;

typedef std::complex<double>                            complex;
typedef ublas::shallow_array_adaptor<double>            shallow_adaptor_double;
typedef ublas::vector<double, shallow_adaptor_double>   shallow_vector_double;
typedef ublas::shallow_array_adaptor<complex>           shallow_adaptor_complex;
typedef ublas::vector<complex, shallow_adaptor_complex> shallow_vector_complex;

// Ensure we can use std::complex<double> as two consecutive doubles
BOOST_AUTO_TEST_CASE( shared_c_array )
{
  // Assumption must hold true for any of this scheme to work
  BOOST_REQUIRE_EQUAL( sizeof(complex), 2*sizeof(double) );

  const std::size_t N = 6;
  double carray_double[N] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };
  complex *carray_complex = reinterpret_cast<complex *>(carray_double);

  for (std::size_t i = 0; i < N/2; ++i)
    {
      BOOST_CHECK_EQUAL(carray_double[2*i],   carray_complex[i].real());
      BOOST_CHECK_EQUAL(carray_double[2*i+1], carray_complex[i].imag());
    }

  for (std::size_t i = 0; i < N/2; ++i)
    {
      carray_complex[i] += carray_complex[0];
    }

  for (std::size_t i = 0; i < N/2; ++i)
    {
      BOOST_CHECK_EQUAL(carray_double[2*i],   carray_complex[i].real());
      BOOST_CHECK_EQUAL(carray_double[2*i+1], carray_complex[i].imag());
    }
}

// Ensure we can use a C-array of doubles as a ublas::vector.
BOOST_AUTO_TEST_CASE( shallow_array_adaptor_double )
{
  const std::size_t N = 3;
  double raw[N];

  shallow_adaptor_double adaptor(N, raw);
  shallow_vector_double  vec(N, adaptor);

  BOOST_CHECK_EQUAL( vec.size(), N );

  std::fill(&raw[0], &raw[N], 1.0);
  BOOST_CHECK_EQUAL_COLLECTIONS(&raw[0], &raw[N], vec.begin(), vec.end());

  std::fill(vec.begin(), vec.end(), 2.0);
  BOOST_CHECK_EQUAL_COLLECTIONS(&raw[0], &raw[N], vec.begin(), vec.end());
}

// Ensure we can use a C-array of std::complex<double> as a ublas::vector.
BOOST_AUTO_TEST_CASE( shallow_array_adaptor_complex )
{
  const std::size_t N = 3;
  complex raw[N];

  shallow_adaptor_complex adaptor(N, raw);
  shallow_vector_complex  vec(N, adaptor);

  BOOST_CHECK_EQUAL( vec.size(), N );

  std::fill(&raw[0], &raw[N], complex(1.0, -1.0));
  BOOST_CHECK_EQUAL_COLLECTIONS(&raw[0], &raw[N], vec.begin(), vec.end());

  std::fill(vec.begin(), vec.end(), complex(2.0, -2.0));
  BOOST_CHECK_EQUAL_COLLECTIONS(&raw[0], &raw[N], vec.begin(), vec.end());
}

// Ensure we can use the same data as either double or complex ublas::vector
BOOST_AUTO_TEST_CASE( shallow_array_adaptor_shared )
{
  const std::size_t N = 6;
  double carray_double[N] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };
  complex *carray_complex = reinterpret_cast<complex *>(carray_double);

  shallow_adaptor_double adaptor_double(N, carray_double);
  shallow_vector_double  vec_double(N, adaptor_double);

  shallow_adaptor_complex adaptor_complex(N, carray_complex);
  shallow_vector_complex  vec_complex(N, adaptor_complex);

  for (std::size_t i = 0; i < N/2; ++i)
    {
      BOOST_CHECK_EQUAL(vec_double[2*i],   vec_complex[i].real());
      BOOST_CHECK_EQUAL(vec_double[2*i+1], vec_complex[i].imag());
    }

  for (std::size_t i = 0; i < N/2; ++i)
    {
      vec_complex[i] += vec_complex[0];
    }

  for (std::size_t i = 0; i < N/2; ++i)
    {
      BOOST_CHECK_EQUAL(vec_double[2*i],   vec_complex[i].real());
      BOOST_CHECK_EQUAL(vec_double[2*i+1], vec_complex[i].imag());
    }

  for (std::size_t i = 0; i < N; ++i)
    {
      vec_double[i] += vec_double[N-1];
    }

  for (std::size_t i = 0; i < N/2; ++i)
    {
      BOOST_CHECK_EQUAL(vec_double[2*i],   vec_complex[i].real());
      BOOST_CHECK_EQUAL(vec_double[2*i+1], vec_complex[i].imag());
    }
}

Suzerain poster from 2009 PECOS center review

2009-09-30T19:00:00.009Z

I spend some time assembling a state-of-my-research poster for my research center's annual review:

DFW IJ/And But So celebration this Saturday in Austin

2009-09-23T04:33:00.002Z

Partly because I stumbled across Infinite Summer shortly after David Foster Wallace died, and partly because I recommended the book to my older sister who bought it and dropped it and passed it on to my father who dropped it but still had the previous copy from several years back when I recommended it and he bought it, dropped it, and now gave his now double-dropped extra to me when I was visiting over the summer, I re-read Infinite Jest this July, August, and September. I finished this morning. Like last time, I'm worried that most anything I pick up will be again unsatisfying to read for a long while.

In commemoration of today's end of Infinite Summer, celebrations abound. In Austin, there's And But So this Saturday the 26th. Attendees bring their favorite Foster Wallace passages to read aloud. I will sadly miss it. If I could attend And But So, I would read this passage from "Lyndon" in David Foster Wallace's short story collection Girl With Curious Hair:

"The truth made the truth's usual quick circuit around the offices, the Building, the Hill. I was a homosexual. I had been a homosexual at Yale. In my last year before matriculating to the Business College, I met and became intimate with a Yale undergraduate, Jeffrey, a wealthy boy from Houston, Texas, who was beautiful, often considerate, wistful, but passionate, possessive, and a sufferer from periodic bouts of clinical depression so severe he had to be medicated. It was the medication, I discovered, that made him wistful.

"My lover Jeffrey ran with a group of synthetic but pleasant Texas socialites, one of whom was Margaret Childs, a tall, squarely built girl who eventually claimed, from unknown motive, to be in love with me. Margaret pursued me. I declined her in every sensitive way I knew. But Jeffrey grew inflamed. He revealed that his friends did not and must not know he was a homosexual. He pushed me to avoid Margaret altogether, which was hard: Margaret, gritty, bright enough to be chronically bored, had become puzzled, suspicious, of Jeffrey's (quite unsubtle) attempts to shield me from her. She smelled potential drama, and kept up the pursuit. Jeffrey become jealous as only the manic can. In my first year in Business, while I was shopping for my father's annual Christmas golf balls, Jeffrey and Margaret had it out, publicly, dramatically, in a Beat New Haven coffeehouse. Jeffrey put his foot through a doughnut counter. Certain information became public. Bits of this public information got back to my parents, who were close to the parents of two of my housemates. My parents came to me, personally, at Yale, on campus. It was snowing. At dinner with my parents and housemates, at Morty's, Jeffrey became so upset that he had to be taken to the men's room and calmed. My father swabbed Jeffrey's forehead with moist paper towels in a cold stall. Jeffrey kept telling my father what a kind man he was.

"Before my parents left—their hands literally on the handles of the station wagon's doors—my father, in the snow, asked me whether my sexual preferences were outside my own control. He asked me whether, were I to meet the right woman, I might be capable of heterosexual love, of marriage and a family and a pillar-type position in the community of my choice. These, my father explained, were his and my mother's great and only wishes for me, their one child, whom they loved without judgment. My mother did not speak. I remember a distanced interest in the steam of my own breath as I explained why I thought I could not and so would not do as my father wished, invoking Fifties' wisdom about deviancy, invoking a sort of god of glands as a shaman might blame vegetable spirits for a lost harvest. My father nodded continually throughout this whole very serious and civil conversation while my mother checked maps in the glove compartment. When I failed to present for next week's holiday, my father send me a card, my mother a check and leftovers in foil.

"I saw them only once more before my father dropped dead of something unexpected. I had left Jeffrey's company, and had been befriended in my upset by a still grimly determined Margaret Childs. Jeffrey unfortunately saw, in all this, cause to take his own life, which he did in an especially nasty way; and he left, on the table beneath the heating pipes from which he was found suspended, a note—a document—neatly typed, full enough of absolute truth concatenated with utter fiction that I was asked by the administration of the Business College to leave Yale University. Weeks after my father's wake I married Margaret Childs, under a mesquite tree, the blue stares of my mother and a Houston sky, and a system of vows, promises of strength, denial, trial, and compassion far beyond the Childs' Baptist minister's ritual prescriptions.

"The truth, to which there was really no more than that, and which made its way through the Senator's staff, the Dirksen, and Owen Buildings, and the Little Congress of the Hill's three-piece-suited infantry remarkably focused and unexaggerated, concluded with the fact that Margaret's father, Mr. Childs, less wealthy than outright powerful by the standards of 1958's Texas, had lines of political influence that projected all the way into the U.S. Senate, and that he, Mr. Childs, in a gesture that was both carrot and stick, slung his son-in-law on one line of that influence and had me hand-over-hand it into the offices of a risen and rising, uncouth and ingenious senior Senator, a possible Democratic candidate in the next Presidential election. Lyndon."

The skew-adjoint convective derivative operator

2009-09-14T22:59:00.001Z

Many CFD papers talk about using the skew-symmetric form of the convector operator in the Navier-Stokes equations. The math behind this terminology has confused me for some time (e.g. the operator's nonlinear in the Navier-Stokes equations), so finally I dug through the algebra to show why and when the linear operator is skew-Hermitian and/or skew-adjoint.

The skew-adjoint (skew-symmetric) form of the convective derivative

Printing arrays and gsl_matrix data in 2D form within gdb

2009-07-16T23:32:00.000Z

Because I'm lousy with indexing, I often have to watch loops unfurl live in gdb. Here are several quick user-defined gdb commands that make inspecting 2-dimensional data stored in a 1-dimensional array much easier.

Place the following commands in your ~/.gdbinit:

# Modified from http://sourceware.org/ml/gdb/2007-07/msg00140.html
define p_array2d_generic
    set $array=$arg0
    set $ibound=$arg1
    set $jbound=$arg2
    set $ilda=$arg3
    set $jlda=$arg4

    set $i = 0
    set $j = 0
    while $i<$ibound
        while $j<$jbound
            printf " "
            printf "%9.4g", $array[$i*$ilda+$j*$jlda]
            set $j = $j + 1
        end
        printf "\n"
        set $j = 0
        set $i = $i + 1
    end
end

document p_array2d_generic
Print an array in 2D form: p_array2d array ibound jbound ilda jlda
end

define p_array2d_cm
    p_array2d_generic $arg0 $arg1 $arg2 1 $arg1
end

document p_array2d_cm
Print an array in column-major 2D form: p_array2d array ibound jbound
end

define p_array2d_rm
    p_array2d_generic $arg0 $arg1 $arg2 $arg2 1
end

document p_array2d_rm
Print an array in row-major 2D form: p_array2d array ibound jbound
end

Now you can use p_array2d_cm and p_array2d_rm to print out column- and row-major double arrays, respectively.

Using the above and some information about the memory layout of gsl_matrix instances, it's now easy to write one line function to dump a gsl_matrix:

define p_gsl_matrix
    p_array2d_generic $arg0->data $arg0->size1 $arg0->size2 $arg0->tda 1
end

document p_gsl_matrix
Print a gsl_matrix in 2D form: p_gsl_matrix (gsl_matrix *)
end

Question for you: I would like to have the p_array2d_generic take an additional argument which is a format specifier. A check using $argc would allow a reasonable, double-based default but would allow the user to specify anything he or she desires. However, I cannot get code like

set $format="%3.2g"
printf $format, 123

to work the way I think it should. I always run into a gdb error like <<Bad format string, missing '"'.>>. Anyone have any ideas how to provide a printf template parameter from a set variable in gdb?

My favorite South Austin photo

2009-07-07T03:38:00.001Z

Quick! Find the following items in the South Austin picture to the right:

tiki torch
pina colada
high volume bangs
proto-dreads
man purse
brick in the middle of the grass
patriotic broach
aviator sunglasses
cheap Mexican beer
a chihuahua wearing a lei
a massive inflatable pool complete with pump

I need to cross the river (lake?) more often.

Skateboarding (Injury)

2009-07-01T04:45:00.003Z

Eight days ago I purchased a skateboard at the age of 27. Three days ago I fell at Mabel Davis park at the age of 27. I messed up an elbow badly enough that I've been in a sling for several days now. I see an orthopedic surgeon tomorrow morning to figure out why I cannot straighten my right arm completely. No harm though, I can still type and use a mouse.

I'm primarily annoyed that I can't roll around again until this injury is healed. I'm awful at the sport, but damn it is fun stuff. Behold shameless YouTubage of the guy who's name is on the bottom of my deck:

Frank Gabron

2009-06-26T16:00:00.009Z

The Penn State Mechanical and Nuclear Engineering department's last alumni newsletter had a small feature on Frank Gabron (article follows). I received one of his scholarships while I was an undergraduate, and it (along with my parents being very cool) was what allowed me to spend long enough on campus to pick up a double major with a double minor. Every year I wrote Gabron a thank you note, and every year he wrote back.

One fall his response said he was coming into State College and would like to have dinner. We ate at the Nittany Lion Inn and made small talk. He told some fantastic stories. After dinner we walked outside and shared that awkward I'm-going-to-smoke-dont-judge-me-oh-no-way-you-smoke-too moment. As we stood and talked he asked if I'd like to come to the Penn State game the following day.

I went back to the Inn for a pre-game pep rally with Graham Spanier firing up everyone. We rode over to the stadium and watched the game from the nicest seats I will ever have in Beaver Stadium—the president's box. At halftime we snuck out for a smoke at a place inside the stadium where he knew he could usually get away with one. He was telling me about his Helix Technology Corporation when a campus cop wandered into the outside area. Without missing a beat, the man said "Oh shit, here comes the gendarme!" and cupped his cigarette in his hand behind his back like a kid in a high school parking lot when the coach drives past.

After the game we rode back to the Inn. I fell asleep on the bus on the way back, woke up hurriedly, did a lousy job of saying thank you, and then ran off to flip burgers at Baby's for a closing shift. I don't remember who Penn State played that day, or even if we won. But I do and always will remember that as one of the two or three best weekends I spent on campus.

Update 3 June 2011: It seems Mr. Gabron passed away on the 17th of May at the age of 81. I still and always will smile every time I think of that football game and his generosity in helping me obtain my undergraduate degree.

PSU MNE Spring 2009 Newsletter: Frank Gabron Article

Trefethen's Gauss Quadrature vs Clenshaw-Curtis Paper

2009-06-18T14:43:00.006Z

From a quick look into B-spline Galerkin methods, I need to integrate B-spline basis functions against each other. Though I could extract piecewise polynomial representations and integrate them symbolically, it was more fun to see what fixed order Gauss quadrature FOSS existed. I stumbled across Pavel Holoborodko's excellent routines. I'm in the process of preparing a patch for the GNU Scientific Library that ports Pavel's code to GSL conventions (update: patch submitted). Hopefully they'll accept it (update: patch accepted).

While digging around, I ran across a beautiful paper by Trefethen which I just finished on the morning bus ride into campus. Absolutely fantastic.

Is Gauss Quadrature Better than Clenshaw–Curtis?
SIAM Rev. 50, 67 (2008)
http://link.aip.org/link/?SIREAD/50/67/1

The paper is so good that it makes me also want to try implementing Clenshaw-Curtis rules for GSL. Before coding, I need to read this follow up work which it seems may have some numerical details buried within it.

New Quadrature Formulas from Conformal Maps
SIAM J. Numer. Anal. 46, 930 (2008)
http://link.aip.org/link/?SJNAAM/46/930/1

Stack Overflow Dev Day in Austin

2009-06-13T21:50:00.002Z

Joel Spolsky and Jeff Atwood will be speaking at Stack Overflow Dev Day in Austin on October 14th. Student tickets to the event are a mere $10 and they provide lunch. Total no brainer.

Grrr... (Richardson Extrapolation)

2009-06-09T20:27:00.003Z

Many times, when banging my ahead against on what has to be a really simple problem, I waste a lot of time over-engineering my debugging tools. Behold: richardson.m

(But then, years later, some of this crap comes up in conversation and the tool is already polished enough to hand it to someone else. At least this is what I tell myself...)

le semaine prochaine avec mes parents

2009-05-24T04:39:00.000Z

There...

View Larger Map

...and back again

View Larger Map

Master's degree and a third patent grant

2009-05-18T00:55:00.001Z

Nice week all in all. I finished an M.S. in Computational and Applied Mathematics which has been four semesters in the making.

Longer than that in process, though markedly less work, was getting my third patent grant US 7,532,617:

A method, apparatus, and computer instructions are provided for a session initiation protocol application design, development, execution, and integration. A framework with a framework controller, an application data object, transition handlers, precondition handlers and post condition handlers is provided. When a message is received in a SIP application, the framework controller determines a state of application and parses the message to update application data by using corresponding handlers. The handlers for each state are specified in an input markup language model defined by a user externally via a Web service accessible component.

Relatively recently...

2009-04-26T02:37:00.002Z

...Michelle had some ankle surgery. She's doing well.
...I saw rms speak on U.T.'s campus about copyright law, where, among other things, he made the case that IP is an overgeneralization. Still digesting that.
...Michelle and I watched Cashback (excellent, too much skin though) and Special (excellent).
...I failed miserably at getting a simple, 1-D version of Durbin's V2F model working on a turbulent channel flow.
...despite it always being an obvious first choice, I keep forgetting to try to use geometric Brownian motion in my probability theory homework.
...I quit letting scope creep into the reading list for my turbulence term paper.
...I realized just how hard it is to achieve meaningful functional and system verification testing work, even when you as a tester understand the application and the underlying technology fairly well.

Turbulence/Martingale Theory Confusion

2009-03-27T02:42:00.000Z

Despite what I may think every Tuesday and Thursday when I walk from my probability theory class into my turbulence course, ν_T is an eddy viscosity and not the stopped viscosity process.

Nearly sufficently advanced technology

2009-03-24T04:16:00.003Z

Stephen Pope's Turbulent Flows is one of the go-to books for graduate courses in turbulence. It flows well, re-reading sections rewards the reader with new material, and the book attempts to develop results instead of just arriving with them from on high. All goodness. Still, every so often, Pope casually drops a bomb to check if you're paying attention. Most of the time I'm not, but tonight I saw one:

Oh, I thought, he's just using a Taylor expansion. No big deal. I kept on reading.

Pope then quietly announces that it was a stochastic Taylor expansion. At least that's my guess at what that object would be called.

I've never seen a “stochastic Taylor expansion” developed, nor do I know the technical constraints that say whether an expansion does or does not exist for random fields like u, v, and w. Any mathematical concept lacking a Wikipedia entry (as of today, anyhow) probably requires a bit more fanfare. A footnote? A reference? A statement declaring überformality? Not Pope's style it seems...

TAMMS 2009 Speakers

2009-03-19T17:07:00.001Z

I just accepted the final student speaker for the 2009 Texas Applied Mathematics Meeting for Students. I'm very happy with the breadth of schools and topics we'll have at the conference. Here are the talk titles firmly on the schedule:

Name	University	Talk Title
Yulia Hristova	Texas A&M University	Time reversal in thermoacoustic tomography - an error estimate
Anthony R. Kellems	Rice University	Dimension Reduction Techniques that Capture Nonlinear Behavior of Morphologically Accurate Neuronal Models
Linh Nguyen	Texas A&M University	On singularities and instability of reconstruction in thermoacoustic tomography
Yan LI	Texas A&M University	Local-Global Upscaling of Flow and Transport in heterogeneous porous media
Ryan Nong	Rice University	Numerical Solutions of Matrix Equations Arising in Dimension Reduction for Linear-Time-Invariant Systems in the Large-Scale Setting
Jay Raol	Rice University	Using Fast Activating Voltage Sensitive Calcium Channels as Voltage Sensors
Moritz Allmaras	Texas A&M University	Ultrasound Modulated Optical Tomography: Reconstructions for a Differential Model
Sean Hardesty	Rice University	Optimization of Shell Structure Acoustics
Dimitar Trenev	Texas A&M University	Approximating solutions of infinite domain Laplace and Helmholtz problems
Robert Rosenbaum	University of Houston	Correlation Propagation in Networks of Integrate-and-Fire Neurons

Using ACX_BLAS before AX_PATH_GSL with gfortran 4.3.2

2009-02-13T04:40:00.001Z

While autotooling a project that uses GSL, I wanted GSL to first attempt to detect and use a system BLAS. If none can be found, then I wanted GSL to fall back to its own gslcblas implementation. The rational being that a savvy user could eke out some performance, but that for everyone else the code just builds and runs without being too fussy.

First I tried this in configure.ac:

ACX_BLAS([
export GSL_CBLAS_LIB="${BLAS_LIBS}"
AC_SUBST(GSL_CBLAS_LIB)
],AC_MSG_WARN([Will use gslcblas if it is present. Try --with-blas=<lib>.]))
AX_PATH_GSL(1.12,,AC_MSG_ ERROR([Could not find required GSL version.]))

This works beautifully with Intel's compilers, but it dies for gcc/gfortran 4.3.2 during configure:

checking for GSL - version >= 1.12... no
*** Could not run GSL test program, checking why...
*** The test program failed to compile or link.

with the linking problems in config.log looking like

configure:7272: gcc -o conftest -I/org/centers/pecos/LIBRARIES/GSL/gsl-1.12-gcc-4.3.1-ubuntu-amd64/include -I/h2/rhys/include -L/h2/rhys/lib conftest.c -L/org/centers/pecos/LIBRARIES/GSL/gsl-1.12-gcc-4.3.1-ubuntu-amd64/lib -lgsl -lcblas -lf77blas -latlas -lm >&5
/usr/lib/../lib64/libf77blas.so: undefined reference to `_gfortran_st_write_done'

which a quick Google tells me is because -lgfortran isn't present in the linking line. I tried to obtain the required gfortran libraries in ${FLIBS} by using AC_FC_LIBRARY_LDFLAGS and then use them via

export GSL_CBLAS_LIB="${BLAS_LIBS} ${FLIBS}"

but later in the build I started hitting up against duplicate main issues because -lgfortranbegin winds up in ${FLIBS}.

More random Googling found the answer buried in lapack++'s configure script. The total solution looks like:

AC_PROG_SED
ACX_BLAS([
dnl Workaround for bogus FLIBS
FLIBS=`echo ${FLIBS} | ${SED} 's/-lgfortranbegin//'`
export GSL_CBLAS_LIB="${BLAS_LIBS} ${FLIBS}"
AC_SUBST(GSL_CBLAS_LIB)
],AC_MSG_WARN([Will use gslcblas if it is present. Try --with-blas=<lib>.]))
AX_PATH_GSL(1.12,,AC_MSG_ERROR([Could not find required GSL version.]))

where there's an explicit hack to remove -lgfortranbegin from FLIBS. Not pretty, but it seems to work.

Texas Applied Mathematics Meeting for Students: March 27-28th

2009-02-09T13:12:00.025Z

I've (finally) gotten to the point where I can announce:

The Austin Chapter of SIAM is proud to host the 3rd annual Texas Applied Mathematics Meeting for Students (TAMMS) on March 27–28th. Attendees will have an opportunity both to present their own research and to meet fellow students from other Texas institutions. Though targeted at SIAM members, this meeting is open to anyone with a mathematical leaning. Please see http://www.ices.utexas.edu/siam for more details.

Compared to the last time I wrote about the site, I've gotten the back button and search engine problems fixed. Or, I've got the back button working fine on Firefox. I'm unsure if dojo.back can be made to work on IE. Though not perfect, the site is usable sans JavaScript. I spent some time looking at a reCaptcha/PHP-based contact form but opted to forgo the contact form in the end.

Compressible Navier-Stokes formulation for a perfect gas

2009-02-09T03:48:00.016Z

It took me longer than I care to admit to track down and work through all the details/assumptions underneath the derivation of the conservative, compressible Navier-Stokes model for a perfect gas. This is on page two of every thesis I read, but is apparently such common knowledge that no one bothers to stick it in an appendix or give a good reference.

Compressible Navier-Stokes formulation for a perfect gas

Publish at Scribd or explore others: Research compressible NavierStokes

Particularly interesting was tracking down data to confirm that air viscosity goes like temperature to a goofy power like 0.666, at least over a reasonable temperature range (e.g. 200K - 5000K).

Avoiding tedious numerics using Boost Accumulators

2009-02-04T05:57:00.042Z

For a turbulence warm up assignment, we needed to compute some correlation coefficients for trajectory behavior on the Lorenz attractor. Specifically, compute the mean, variance, and cross-correlations of the three solution components. I'd skimmed through the Boost Accumulator Framework a couple of months back and thought it looked interesting:

#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/covariance.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/variance.hpp>
#include <boost/accumulators/statistics/variates/covariate.hpp>
#include <boost/format.hpp>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>

// Read an input streams and compute the relevant statistics
template<typename charT, typename traits>
int process(std::basic_ostream<charT, traits>& out,
            std::basic_istream<charT, traits>& in) {

    using namespace boost::accumulators;
    accumulator_set<double, stats<tag::mean,
                                  tag::variance,
                                  tag::covariance<double, tag::covariate1> > > x_acc;
    accumulator_set<double, stats<tag::mean,
                                  tag::variance,
                                  tag::covariance<double, tag::covariate1> > > y_acc;
    accumulator_set<double, stats<tag::mean,
                                  tag::variance,
                                  tag::covariance<double, tag::covariate1> > > z_acc;

    double t, x, y, z;
    std::basic_string<charT, traits> line;
    boost::format fmt(" %016e   %016e   %016e   %016e   %016e   %016e   %016e   %016e   %016e   %016e");
    while (in) {
        if (!std::getline(in, line)) break; // Read columns per line to
        std::istringstream in(line);       //    be a little defensive
        in >> t >> x >> y >> z;

        // Compute running statistics
        x_acc(x, covariate1 = y);
        y_acc(y, covariate1 = z);
        z_acc(z, covariate1 = x);

        // Output running statistics
        out << fmt % t
                   % mean(x_acc)
                   % mean(y_acc)
                   % mean(z_acc)
                   % variance(x_acc)
                   % variance(y_acc)
                   % variance(z_acc)
                   % covariance(x_acc)
                   % covariance(z_acc)
                   % covariance(y_acc)
            << std::endl;
    }

    return 0;
}

Compared with the alternative loops/bookkeeping/etc., this is a pretty slick alternative. The only downside was that the accumulator declarations took me awhile to decipher from the documentation.

As you probably noticed, I am reading/writing plain text data files and computing the statistics at every time step. This is overkill and slow. I'd probably have fixed that by now if it weren't for the fact that I ran across Pipe Viewer in the last two days. Quite a neat shell tool that.

Jet Plane Off To Jamaica

2009-01-26T00:09:00.000Z

Two weeks back Flip and I converged on Trapani's apartment in Jamaica, Queens. The three of us were together in Simmons Hall at Penn State, near Plateia Varnava in Athens, near Beaver and Atherton back at Penn State, and finally in Austin, Texas. Eight years.

Over the past year and a half both Flip and T left town for new opportunities. It's been hard to not see them nearly every Wednesday night like I had for the three years we were all in Austin. This trip we sat, talked, laughed, and did nothing of merit for three days. It was wonderful.

Austin Chapter of SIAM website redesign

2009-01-20T23:55:00.001Z

I finally got around to redesigning U.T.'s SIAM chapter website. I've been meaning to do it for months. The suitable effective kick in the pants was needing a site to house information for our upcoming student conference.

All in all I'm happy with it. I used

some simple Dojo loaded via Google's AJAX hosting,
a little Blueprint,
lifting/revising content from both Pearl Flath's older chapter site and the official SIAM pages, and finally
my feeble Inkscape logo design.

My three major complaints with my own work are:

The back button is broken.: I spent two hours futzing with dojo.back to get the back button to work. After I had it working, hitting back and then forward caused an infinite loop on page refresh. I gave up for now.
None of the content can be indexed by a search engine (I think).: Since all of the major content is loaded via dijit.layout.ContentPane.setHref fired with onClick events, I doubt a search engine can effectively index much of it. Something to check when I roll Google Analytics into the next update.
I have raw email addresses in mailto links.: I'm not sure how much server-side processing is permitted on this particular host. I need to move contact functionality into forms guarded by CAPTCHAs. Hopefully GMail's filter will save me in the meantime.

None of these are terrible to fix. They just take time. And classes began today...

DOE CSGF Research Statement

2009-01-14T17:01:00.000Z

I just finished applying for the DOE Computational Science Graduate Fellowship. In March they tell me if I got it. One of the short responses the CSGF application required was the first time I'd sat down and really thought through my research goals. With the addition of some background material links, here's the prompt and my response...

Essay prompt 1: Field of Interest

"Computational science" involves the innovative and essential use of high-performance computation, and/or the development of high performance computational technologies, to advance knowledge or capabilities in a scientific or engineering discipline. Please describe (in no more than 300 words) your specific research interest paying particular attention to how computational science will spur advances.

Essay response 1

My research uses direct numerical simulation (DNS) of the Navier-Stokes equations to investigate the compressible, turbulent boundary layers present in atmospheric reentry vehicle applications. My work is motivated by both its long-term value and an immediate need at the University of Texas' Center for Predictive Engineering and Computational Sciences (PECOS). PECOS is developing new uncertainty quantification techniques through the study of reentry vehicle thermal protection systems (TPS). Due to its order one impact on vehicle heating rate uncertainty, improving hypersonic turbulent boundary layer science is crucial for PECOS' success.

I have three interwoven research thrusts. First, by combining scalable spectral numerical techniques with a multiscale formulation, I will advance understanding of how TPS ablation products interact with turbulent boundary layers. Second, through a simulation-based characterization of the ablative conditions that may sustain turbulence, I will improve understanding of turbulence transition phenomena in transpired-wall compressible flows. Lastly, I will produce high quality data fields for use in compressible turbulence model development and calibration. This novel data will also fuel PECOS' Bayesian-inspired uncertainty quantification activities.

Improvements in computational science will aid all three of these endeavors. The scientific need to attain higher, more realistic Reynolds numbers and to obtain finer-grained, larger datasets drives DNS towards ever higher performance computing platforms. Increases in compute node availability, compute node memory, and interconnect speed will all boost the quality of the science I produce. Spectral turbulence simulations walk the line between being CPU and communications bound. Accordingly, my research stands to benefit from either improvements to low-level building block performance (e.g. BLAS) or algorithmic advances in high-level scalability (e.g. communication topology optimization). Finally, in turn, my application's stringent requirements will help me drive new computational developments for others to use.

UAW/Republican auto bailout insanity

2008-12-12T21:12:00.004Z

Consider the following snippet (emphasis mine):
U.A.W. in Center of Dispute Over Bailout Failure - NYTimes.com

Late Thursday, the Senate did not take up an assistance measure passed by the House, after hours of negotiations between Senate Republicans with the auto companies and the U.A.W. The sticking point apparently was the union’s refusal to agree to lower wage and benefit rates as soon as next year.

Representatives for the union, which had already accepted a series of cuts in its current contract, sought instead to push any more concessions back to 2011, when the U.A.W.’s contract with Detroit auto companies expires.

In a statement Thursday night, the union said it was “prepared to agree that any restructuring plan should ensure that the wages and benefits of workers at the domestic automakers should be competitive with those paid by the foreign transplants. But we also recognized that this would take time to work out and implement” using programs like buyouts and early retirement offers to bring in new workers at lower rates.

It can be (albeit simplistically) be translated into:

"Yeah yeah yeah. We know our non-market salaries are high enough that it's driving the company down. Yeah, the union will accept pay cuts. But we want bought out/paid off before any of us actually have to work for market wages. Oh, and by the way, millions of other people's jobs be damned."

Interesting that the Republicans killed it, and now that it's Bush who may come swooping in from left field. Especially funny/sad is that the executive action to save things will probably give the U.A.W. what they want. That's likely to overwhelm the domestic automakers yet again even if a current bailout can save them.

Nice Nutcracker article on Michelle in the Austin Decider

2008-12-09T12:45:00.000Z

The Austin Decider ran a fun piece on the role of snow within Ballet Austin's Nutcracker. Michelle was the chosen snowflake:
House Lights: Backstage at The Nutcracker

Fall 2008 :: Probability Theory I Concepts

2008-12-07T23:55:00.000Z

I'm a week away from finishing the first half of a year-long probability theory sequence at U.T. I (finally) managed to keep a concepts cheat sheet up to date for a course:

Probability Theory I: CAM 384K Concepts

Publish at Scribd or explore others: Mathematics Science probability statisti

Part of my motivation has been that while Durrett's text is useful, it is very difficult to find specific definitions and theorems within it while working problems.

Wasted Days

2008-11-25T04:53:00.000Z

Today I received a questionnaire asking for total number of contributions in the following categories over the past academic year:

Refereed Papers
Non-refereed Papers
Books
Book Chapters
Presentation at learned societies/symposia
Invention disclosures
Copyrighted Software
Patents Filed
Other Scholarly Works

to each of which I responded with a big, fat zero. It is early for resolutions, but I aim for '09 to be full of accomplishments that can be measured instead of my usual results (e.g. "Oh! Duh! That's what that equation means. I should'a known that already...").

A B-spline verification test case

2008-11-19T22:29:00.001Z

Update (16 Dec): The B-spline derivative code made it into GSL's 1.12 release

I'm trying to get B-spline derivative capabilities added to gsl. Though I could find lots of implementations of B-spline derivatives for fixed polynomial order, finding derivatives for arbitrary order B-splines turns out to be a gap in many libraries (with PPPACK being the notable exception).

It also turns out that finding a good, concrete verification test case is a pain too. Hopefully gsl will pick up my patch and no one will have to implement these again. But, if you should find yourself in that situation, here's some Mathematica that provides both symbolic and numeric results for quadratic B-splines and their derivatives on a particular set of knots:

A verification case for B-Spline implementations, including derivatives

Get your own at Scribd or explore others: Science Mathematics bspline pppack gsl m

Handcrafted memory layout: a case study in magic number mayhem

2008-11-12T18:18:00.001Z

Observe, in the code snippet below, that no fewer than 27 subtly different aliases have been created for the same chunk of memory...

Semantics/Syntax: The @ sign instructs the compiler to have the variables named varsc and vars share storage much like C's union construct. The symbolic keyword acts like a macro definition.

    dynamic complex(64) varsc[(nx+2)/2,ny,37] 
          @ real(64)    vars[     nx+2,ny,37]   
    
    symbolic mic[i,j,l]        = varsc[i,j,l]
           , rhoc[i,j]         = varsc[i,j,4]
           , engyc[i,j]        = varsc[i,j,5]
           , drho_dxic[i,j,l]  = varsc[i,j,l+5]
           , dengy_dxic[i,j,l] = varsc[i,j,l+8]
           , dmi_dxic[i,j,l,m] = varsc[i,j,11+(l-1)*3+m]
           , rhsc[i,j,l,m]     = varsc[i,j,25+(l-1)*3+m]
           , rhs_rhoc[i,j]     = varsc[i,j,21]
      
    symbolic mi[i,j,l]         = vars[i,j,l]
           , rho[i,j]          = vars[i,j,4]
           , engy[i,j]         = vars[i,j,5]
           , drho_dxi[i,j,l]   = vars[i,j,l+5]
           , dengy_dxi[i,j,l]  = vars[i,j,l+8]
           , dmi_dxi[i,j,l,m]  = vars[i,j,11+(l-1)*3+m]
           , rhs_rho[i,j]      = vars[i,j,21]
           , dui_dxi[i,j,l,m]  = vars[i,j,11+(l-1)*3+m]
           , pres[i,j]         = vars[i,j,22]
           , tmpt[i,j]         = vars[i,j,23]
           , mu[i,j]           = vars[i,j,24]
           , tau[i,j,l,m]      = vars[i,j,25+(l-1)*3+m]
           , q1[i,j]           = vars[i,j,9]
           , q3[i,j]           = vars[i,j,11]
           , rhs[i,j,l,m]      = vars[i,j,25+(l-1)*3+m]
           , dmu_dx2[i,j]      = vars[i,j,6]
           , dtmpt_dx2[i,j]    = vars[i,j,7]
           , d2tmpt_dx22[i,j]  = vars[i,j,8]
           , d2ui_dx22[i,j,l]  = vars[i,j,8+l]

Overnight trip to Mason, TX

2008-11-10T00:24:00.001Z

I just got back from an overnight motorcycle trip with my friend Paul.

We'd wanted to travel further west than Mason, but a late start and a large helping of pie in downtown Mason caused us to rethink our plan. Hanging out at The Old Peanut Mill Steakhouse and Fort McKavett Social Club filled most of Saturday evening. Nice place.

The roads that crisscross the Guadalupe river are beautiful this time of year:

The VFR turned over 50K miles on the trip. Paul, as always, did a great job with the route and navigating. Both days were a mix of unused high speed sweeping turns along with some intermittent choppy surface connections across ridges and into each successive valley.

Good? Bad? Over.

2008-11-05T04:38:00.000Z

Done deal:

I've got to admit I'm pleased with the outcome. Hopefully state and local races won't shock the system so badly that they inhibit smooth day-to-day operation.

Two more of my IBM patent applications published

2008-11-02T23:14:00.000Z

I just noticed that these two patent applications just appeared publicly...

(20080198748) System and Method for Burst Traffic Smoothing for SIP Processing Elements: A system and method for burst traffic smoothing for Session Initiation Protocol (SIP) processing elements are provided. A dispatch queue management engine determines whether a received packet is a TCP or UDP packet. If the packet is a TCP packet, the packet is automatically added to the dispatch queue. If the packet is a UDP packet, a value for a drop function f is generated and a random or pseudo-random number r is generated. If r has a predetermined relationship to f, then the UDP packet is added to the dispatch queue, otherwise the UDP packet is discarded. The value for f is based on the current dispatch queue load, the network quality, the retransmission rate, and the allowable drop rate. Thus, the determination as to whether to drop UDP packets or not is configurable by an administrator and also adaptable to the current network and dispatch queue conditions.
(20080229245) Multiple Sorting of Columns in a Displayed Table in a User Interactive Computer Display Interface Through Sequential Radial Menus: Multi-sorting of displayed columns representative of a tabular display is carried out by displaying a table of a plurality of columns, selectively displaying a first radial menu having a plurality of sectors, each sector representative of one of the plurality of columns, enabling a user to select one of the sectors, and responsive to a user selection of a sector for displaying a second radial menu of the plurality of sectors wherein the selected one sector is disabled, e.g. eliminated. This is continued through a sequence of subsequent radial menus until the user has selected the intended set of sequential columns for the multiple sorting.

The second one, because I worked with a friend whose last name starts with a W, is probably the first (and only) application where I'll ever see 'Ulerich et al.'

Compressible DNS within PECOS

2008-11-02T22:55:00.002Z

On Monday a bunch of DOE people are in town for the 1st annual review of the Predictive Engineering and Computational Sciences (PECOS) center at UT Austin. Because I'm partially funded through the center, another student and I put together a poster for this Monday afternoon:

Just in case embedded Scribd hiccups, here's a link to the poster.

U.T. needs to ditch driving on campus

2008-09-28T14:24:00.002Z

This past week U.T.'s Parking & Transportation Services proclaimed

...Beginning Wednesday, September 24th, there will be a Travis County Sheriff
Officer located at the intersection of Speedway and 24th Street. This
officer is being placed at this intersection to provide a secure crossing
environment for pedestrians and vehicular traffic....

As safety is our primary concern, all pedestrians in this area should heed
the warning signs and utilize sidewalks for travel, including the east side
of Speedway between Dean Keeton and 24th Street and the south side of 24th
Street between Speedway and Inner Campus Drive. Please do not walk in the
street or along the closed west sidewalk area of Speedway as construction
traffic in this area is very heavy and pedestrians create potentially
hazardous situations....

Having just spent two days at ACL with tens of thousands of pedestrians, I'm pretty sure the people on foot aren't the most preventable root cause of hazardous situations. If some smelly hungover dude accidentally staggers sideways and body checks you while you cross 24th street, you'll still make it to class. Promise.

The problem is that U.T. stupidly encourages driving across 24th street. Take a look. Those cars? Those aren't construction. The traffic consists of people parking and physical plant vehicles. The intersection becomes hazardous when impatient drivers nudge through a sea of thousands of students on foot mixed with bicycles. Guys with hard hats driving triple diesels aren't the problem vehicles. It's the regular cagers expecting to treat a college campus like any other street.

Students should, of course, be courteous to drivers. But cross campus traffic should be forced to travel north two blocks, across the four lane 26th street, and then two blocks south to San Jacinto. U.T. is not an urban campus, and it should capitalize on that fact by creating pedestrian-only areas where possible. 24th and Speedway would be a great place to start.