D. Adding to vnl_algo

The strategy adopted for converting and wrapping the fortran files is a little involved. Some routines are simple to do, others very tricky. The general procedure is as follows. These steps are elaborated upon in the example below.

D.1 Overview

1. Use GAMS to find the module name, in SLATEC if possible, although CMLIB and TOMS routines are also public domain and good.
2. Convert the fortran to C using f2c
3. Add the routine to the Imakefile in the netlib library.
4. Encapsulate the routine in a class in vnl, after determining a suitable interface.
5. Read the module documentation and determine the calling sequence.
6. In the calling method, create all necessary workspace arrays and temporary variables that the call requires, call the external routine, and convert the results into the classes that VXL expects.
7. After the call, interpret the error code, and handle accordingly.

D.2 Problems

There are a few potential sources of difficulty, mostly in item 5, but in general I find that gritting one's teeth and guessing is a surprisingly good strategy. The main points to remember are:

1. All scalar variables are passed by reference. This means that you need to store all constants in variables and pass their addresses or declare the routines as accepting references. I do the latter for input variables, and the former for outputs.
2. Fortran arrays start from 1 rather than 0. This is actually a non-problem, as f2c generates code which interfaces zero-based to one-based arrays using the Numerical Recipes trick of decrementing the pointer, but is mentioned here for the benefit of fortran programmers.
3. Fortran arrays are stored column-wise rather than row-wise. Class vnl_fortran_copy provides an easy and efficient way to transpose matrices before calling.

In addition to these fortran specifics, it is important to be aware of the sorts of design patterns seen in numerical code. Many routines are coded for maximum generality and efficiency, which can make reading the descriptions heavy going. Common conventions are:

• An array is passed with three dimensions: number of rows in the physical array, number of rows to use for computation and number of columns. This allows the routines to be used on any submatrix of a larger matrix.
• Output results overwrite the input matrix.
• Output results are stored in some compact form, which must be decoded before use. Note however that in many cases routines are supplied to perform further computations using the encoded representation for time and space efficiency. The new QR class will demonstrate how to take advantage of this.
• The results of pivoting are generally returned in vectors of integers (say ipvt), where ipvt[i] is the position to which row/column i has been moved. These permutations which must be applied to the results in order to complete the process.

D.3 Example conversion - QR decomposition

Given the need for an algorithm that is not yet included in the vnl package, say a routine to compute the QR decomposition, your first stop is the GAMS decision tree. Class "D" is Linear Algebra, and class "D5" is QR decomposition. The SLATEC implementation is called DQRDC (Double precision QR DeComposition). Download the source, or obtain it from a local SLATEC distribution. Convert it to a C source file and a prototype file using

   f2c -P dqrdc.f

and from the prototype file dqrdc.P we find that the function prototype is

   int dqrdc_(doublereal *x, integer *ldx, integer *n, integer *p,
              doublereal *qraux, integer *jpvt, doublereal *work,
              integer *job);

At this point, the header of the fortran file dqrdc.f is examined in order to determine the meaning of the parameters. Considering parameter X, we find

   X     DOUBLE PRECISION(LDX,P), where LDX .GE. N.
         X contains the matrix whose decomposition is to be computed.

This means that X is a LDX row by P column matrix, and that we require a decomposition of the first N rows. This is a common convention in fortran programs which allows computation on subblocks of matrices. In general, we will assume that we wish to work on the full matrix, and therefore that LDX = N. To create the required transformed copy of the matrix, use class vnl_fortran_copy:

   vnl_fortran_copy Xtranspose(X);

Now, the function may be called as

   int n = X.rows();
   int p = X.columns();
   vnl_vector<int> jpvt(p);
   jpvt.fill(0); // Mark all columns as pivotable
   vnl_vector<double> work(p);
   int do_pivoting = 1;
   vnl_vector<double> qraux(p);
   dqrdc_(Xtranspose, &n, &n, &p,
          qraux.data_block(), jpvt.data_block(), w.data_block(),
          &do_pivoting);

This document was generated on May, 1 2013 using texi2html 1.76.