Physics 6720: Using LAPACK with C++ or C with SUN systems

A beginning programmer can extend his/her computational power considerably by learning how to use scientific subroutine packages effectively. LAPACK is an example of such a public domain package. It contains mostly linear algebra routines, so is especially useful for solving eigenvalue problems, solving linear systems of equations by direct methods, and doing LU decompositions, singular value decompositions, etc. There are other packages, many of them proprietary and more versatile or especially tuned to a particular architecture, including IMSL, NAG, ESSL.

Learning to use the package requires a little effort, especially for C++ programmers, because the code was originally written in Fortran. It requires some care in interpreting and passing arguments. But it is worth the effort.

There are hundreds of routines to choose from in the LAPACK packages. Our computer laboratories have online documentation in the form of man pages for each of them. The documentation is intended mainly for Fortran users, but there is some help for C and C++ programmers.

We will use the SUN repackaging of LAPACK in our laboratories. Here are the steps with the current installation:

Example


Suppose you wanted to solve a general system of linear equations Ax = b, where A is an n x n square matrix and x and b are n-element column vectors. You decide to use the routine dgesv. The man page abstract obtained with 
     man -M /opt/SUNWspro/man 3p dgesv
shows the C(C++) declaration: 
     void dgesv(int n, int  nrhs,  double  *da,  int  lda,  int
                 *ipivot, double *db, int ldb, int *info) ;
The abstract explains that n is the number of linear equations, nrhs is the number of columns of the rhs matrix B (so you want nrhs = 1), da is a double precision array of dimension (lda,n), and lda is the "leading" dimension of the array da. In C(C++) the subscript convention is reversed, so lda will be the second dimension, as shown below. Continuing with the arguments, ipivot is an output array indicating the rows of the matrix that were interchanged upon pivoting. If this information is of no interest to you, the array must still be supplied, but the result can be ignored. Next db is a double precision array of dimension (ldb,nrhs) specifying the rhs upon entry and containing the solution upon return, ldb is its "leading" dimension, and info tells whether the solution was successful.

Armed with this information we write our code as follows: 

#include <iostream>
#include <sunperf.h>
#define MAX 10

int main(){
   // Values needed for dgesv
   int n;
   int nrhs = 1;
   double a[MAX][MAX];
   double b[1][MAX];
   int lda = MAX;
   int ldb = MAX;
   int ipiv[MAX];
   int info;
   // Other values
   int i,j;

   // Read the values of the matrix
   std::cout << "Enter n \n";
   std::cin >> n;
   std::cout << "On each line type a row of the matrix A followed by one element of b:\n";
   for(i = 0; i < n; i++){
     std::cout << "row " << i << " ";
     for(j = 0; j < n; j++)std::cin >> a[j][i];
     std::cin >> b[0][i];
   }

   // Solve the linear system
   dgesv(n, nrhs, &a[0][0], lda, ipiv, &b[0][0], ldb, &info);

   // Check for success
   if(info == 0)
   {
      // Write the answer
      std::cout << "The answer is\n";
      for(i = 0; i < n; i++)
        std::cout << "b[" << i << "]\t" << b[0][i] << "\n";
   }
   else
   {
      // Write an error message
      std::cerr << "dgesv returned error " << info << "\n";
   }
   return info;
}
Note that we are not required to use the same names as given in the prototype.

Fortran to C(C++) Array Conversion


While the man pages provided in the SUNWspro release of LAPACK give the sunperf.h prototype you need to determine how to call the subroutine from C++ or C, the documentation in these pages is based on Fortran conventions for arrays. The storage and subscripting conventions in Fortran are different from those in C++ and C. These conventions must be taken into account when loading data into the arrays. For example, to use the SSBGV routine for computing eigenvalues and eigenvectors of a generalized symmetric-definite banded matrix, we start with a symmetric matrix A(i,j) and must load a compact matrix AB with its diagonal and off diagonal bands. The packing instructions for AB say 
AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j
These instructions are appropriate for Fortran arrays and must be translated for use in C++ and C. The rules are simple:
  1. Subtract 1 from each subscript
  2. Take the transpose.
This transformation results in 
ab[j-1][ka+i-j] = a[j-1][i-1] for  max(1,j-ka)<=i<=j
We have changed to the C++ and C bracket subscript notation and use lower case names to emphasize that we are now speaking of C++ and C arrays. The subscripts i and j are dummy variables, so we are allowed to redefine them by adding 1 to each everywhere they appear. Making that replacement then gives 
ab[j][ka+i-j] = a[j][i] for  max(1,j-ka+1)<=i+1<=j+1
Of course, since a is symmetric, we could also write a[i][j]. With a little algebra the inequalities can also be rewritten more neatly, giving the final result: 
ab[j][ka+i-j] = a[i][j]  for  max(0,j-ka)<=i<=j.
The array ab must be declared with 
   float ab[MAX][ldab]
where MAX >= N and ldab >= ka. It is passed to the subprogram as the pointer &ab[0][0].
Physics Department Home Page

This page is maintained by:

Carleton DeTar Mail Form Last modified 1 November 2007