A sparse matrix calculator

For the past few months I've been working quite hard on both the ctraj and "libpetey" software.  I hope to have new releases of both within the next month.

One of the most significant additions is a sparse matrix calculator.  This tool allows you to perform operations on sparse matrices much as one might do to simple floating point numbers on a pocket calculator or a calculator app.  This addition is so significant, in fact, that I'm consider creating a whole new project to contain it.

Two of the "peteysoft" projects already output sparse matrices.  The Gaspard-Rice simulator, gr_scat, can output a sparse matrix representing a discretization of the Laplace or Helmholtz equation which models the wave scattering behaviour of the system.  Also, the ctraj tracer simulation outputs a series of sparse matrices that represent the tracer dynamics.  To model a tracer, the initial tracer configuration is multiplied through with these matrices.  Storing them as an array is both more space and time efficient than multiplying them through.

I have already used the matrices output from ctraj in a quite sophisticated analysis: by correlating the singular vectors with a series of sparse measurements (e.g. from a satellite), they can be used in a type of dynamical tracer reconstruction I call "PC proxy."  I submitted a paper about PC proxy analysis to Atmospheric Environment--it can be found here.  The genesis of the sparse calculator project comes about because the reviewers were interested in learning in more detail how the method works, what makes it better than existing methods and how better to parametrize it.  Since I don't really understand the method all that well myself, I decided to write this software tool to conduct a series of numerical experiments.

Sparse matrics can be used to model a large class of partial differential equations.  In the past, the process of matrix multiplication would be coded directly in Fortran.  A similar method is used to compute the tangent model for generating singular vectors from numerical weather models. Each time a matrix multiplication is desired, to be fed back to the eigenvalue subroutines, the model is re-run from scratch. Porting them to a data structure instead is not only faster, it is far more versatile.  Since the matrix is typically quite sparse, the extra storage is not particularly burdensome.

The basic idea behind the tool is that files generated from other programs are imported and used as variables.  Because these files can become quite large, the matrices are only loaded into RAM when they are used in a calculation.  Except for scalars, variable names always correspond to a file name.

Using the ctraj package, we can generate a series of sparse matrices which represent the tracer dynamics for a given range of dates, say May 1, 2000 to May 5, 2000, and a given vertical level, say 500 K.  First, you will need to download the NetCDF files containing NCEP temperature and wind data: air.2000.nc, uwnd.2000.nc, vwnd.2000.nc. Next, extract and convert the desired velocity data using the nc2ds command:

  > nc2ds -i 2000/5/1 -f 2000/5/6 -T 500 vfield_n.ds
  > nc2ds -- -i 2000/5/1 -f 2000/5/6 -T 500 vfield_s.ds

Next, we run the tracer model on this data:

  > ctraj_tracer vfield_s.ds vfield_n.ds tracer.dat

Now we can find the top 10 singular vectors of the tracer dynamics using sparse_calc and put them in a file called sv.vec:

  > echo "sv.vec=svd(sparse_array(tracer.dat)*transpose(tracer.dat), 10)"\
      | sparse_calc -A 50

Finally, plot and display the first singular vector:

  > plot_frame.sh 0 sv.vec sv0.ps

  > gv sv0.ps

Now lets say you want to do the same thing, but only for the first day's
worth of data:

  > echo \
"sv_1day.vec=svd(sparse_array(tracer.dat(0:4))*transpose(tracer.dat(0:4)), 10)"\
| sparse_calc -A 50

The sparse matrix calculator will be available in the next release of libpetey , currently a part of the msci project.  The tool can handle generalized sparse matrices, sparse matrix arrays, full matrices and simple floating point calculations.  There are routines for matrix solution, inversion and eigenvalue calculation.  In addition, there is a generous online help function as well as a man page to quickly get you up to speed using it.