Physics 6720: Using LAPACK with C++ or C with CentOS Linux

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. See instructions above for finding the one you want.

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 dgesv
gives the Fortran definition and explains the meaning of the variables. We don't reproduce the explanation here. 
       SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )

           INTEGER       INFO, LDA, LDB, N, NRHS

           INTEGER       IPIV( * )

           DOUBLE        PRECISION A( LDA, * ), B( LDB, * )
Following instructions in the Fortran binding handout, we write the C++ prototpye declaration as follows: 
extern "C" {
     void dgesv_(int *n, int *nrhs,  double *a,  int  *lda,  
           int *ipivot, double *b, int *ldb, int *info) ;
}
The man page 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), a is a double precision array of dimension (lda,n), and lda is the "leading" dimension of the array a. 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 b 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>

extern "C" {
     void dgesv_(int *n, int *nrhs,  double *a,  int  *lda,  
           int *ipivot, double *b, int *ldb, int *info) ;
}



#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


The documentation in our man 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. Notice that we didn't subtract 1 from the inequality expression. Subtracting 1 applies only to what is inside the subscript brackets []. 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
(Notice that we make this replacement everywhere, including the inequality, because we are redefining the dummy variables). 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 29 March 2017