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.

- Find the LAPACK routine you want.

The University of Colorado Computer Science department provides a search tool for finding the routine you want. Click on "driver routines". Your browser needs to be able to handle Java applets for this tool. See unlockingJava.html if your browser won't run this tool because of security concerns.Netlib.org, which hosts LAPACK provides another handy search tool . To use it, in the left panel, click on "LAPACK->Modules->LAPACK". You will get a graphic. For a plain, garden-variety matrix, click on the "General Matrices" box. Then follow the other boxes until you get a list of routines. You'll need to scan down until you find the one appropriate for your problem.

The Lighthouse project is developing a search tool for a wide variety of scientific software tools. You may access it as a "guest", but they like you to register.

Alternatively, you may use the LAPACK manual on the netlib site. Find the table of contents and go to the "Index of Driver and Computational Routines". Note that to get the names of the equivalent double precision routines, replace the first S with a D. To get complex single precision, use C, and complex double precision, use Z.

- Understand what is inside the black box

Be sure you understand the numerical methods used in the package routine, and know their limitations.

- Find out what arguments are needed for the routine you have
chosen. Use the Netlib documentation listed above.
There you will see the Fortran specification. You will
need to construct a prototype declaration for C++ based on what you
see and insert the declaration at the top of your code. Instructions
are given in the Fortran binding
handout. An example is given below.

- Compiling and linking

If you are compiling without also linking, you do as usual:g++ -O -c myprog.cc

To link your code with the LAPACK library, do this:g++ myprog.o -o myprog -L/usr/local/lib64 -llapack -lblas -lgfortran -lm

*Be sure to put all of the g++ options on one line*: Of course compilation and linking can be done in a single step withg++ -O myprog.cc -o myprog -L/usr/local/lib64 -llapack -lblas -lgfortran -lm

again,*be sure to put all of the g++ options on one line*.

- Be skeptical

It is always a good idea to check your results against problems with known solutions before trusting that you have understood how to make the black box work.

Suppose you wanted to solve a general system of linear equations

man dgesvgives 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

Armed with this information we write our code as follows:

#include <iostream>Note that we are not required to use the same names as given in the prototype.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; }

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

AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=jThese instructions are appropriate for Fortran arrays and must be translated for use in C++ and C. The rules are simple:

- Subtract 1 from each subscript
- Take the transpose.

ab[j-1][ka+i-j] = a[j-1][i-1] for max(1,j-ka)<=i<=jWe 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

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

ab[j][ka+i-j] = a[i][j] for max(0,j-ka)<=i<=j.The array

float ab[MAX][ldab]where

Physics Department Home Page

*This page is maintained by:
*