LA Pack: Linear Equation package

What is LA Pack

LA Pack is a collection of Fortran 77 functions that can help you solve Linear Algebra related problems
The library has been converted to Fortran 90

From the LA Pack Internet Interface and Search engine website:

LAPACK

systems of simultaneous linear equations

least-squares solutions of linear systems of equations

eigenvalue problems

singular value problems

The associated matrix factorizations (LU, Cholesky, QR, SVD, Schur, generalized Schur) are also provided, as are related computations such as reordering of the Schur factorizations and estimating condition numbers.

Dense and banded matrices are handled, but not general sparse matrices. In all areas, similar functionality is provided for real and complex matrices, in both single and double precision.

LAPack = Linear Algebra PACKage

LA Pack Resources on the Internet:

LAPack homepage : http://www.netlib.org/lapack/ ( click here )
LA Pack online user guide: click here
Man page s for all LAPack routines (gzip tar archive): Click here
After download, use gunzip manpages.tgz and tar xvf manpages.tar to extract the files.
The man pages are written in "ROFF". To format the man pages, use "psroff -man xxxx" or "nroff -man xxxx"
LAPACK Internet Interface and Search Engine (VERY USEFUL): Click here
A quick overview (where I got the demo programs): click here (link is gone now :-()

Types of subroutines in LAPack

There are three different types of LA PACK routines.

Simple Driver Routines: click here
- Simple driver routines solve a complete problem
- Examples:

Expert Driver Routines: click here

Expert driver routines do the same things as simple driver routines, but will provide more options or information to the user

Example:

SGESV is used to solve linear systems
SGESVX not only solves the linear system, but will also provides an estimate of the condition number of the input matrix and provides error bounds on the solutionx.

Computational Routines
- These routines are mainly for internal use by LA Pack itself.
  These routines are called by the driver routines to perform a specific computational task
- Examples:
- These routines are only necessary when you cannot find any driver routine in LA Pack to solve your problem.
  You must implement your own driver program (thus extending LA Pack to solve your problem)
  
  The computational routines are organized into 3 categories :
  - systems of linear equations : click here
  - eigenvalue problems: click here
  - orthogonal factorization and singular value decomposition: click here

Naming schemes for LAPack routines

LA Pack routines still uses the F77 naming constrain
F77 restricts names to 6 characters
The name of each LAPACK routine is a coded specification of what kind of function the routine will perform.

An LAPack function has the name of the form:

XYYZZZ

The encoding scheme used is:

X = type of problem that the routine solves S Single precision real D Double precision real C Single complex Z Double precision complex YY = Matrix types BD bidiagonal DI diagonal GB general band GE general (i.e., unsymmetric, in some cases rectangular) GG general matrices, generalized problem (i.e., a pair of general matrices) GT general tridiagonal HB (complex) Hermitian band HE (complex) Hermitian HG upper Hessenberg matrix, generalized problem (i.e a Hessenberg and a triangular matrix) HP (complex) Hermitian, packed storage HS upper Hessenberg OP (real) orthogonal, packed storage OR (real) orthogonal PB symmetric or Hermitian positive definite band PO symmetric or Hermitian positive definite PP symmetric or Hermitian positive definite, packed storage PT symmetric or Hermitian positive definite tridiagonal SB (real) symmetric band SP symmetric, packed storage ST (real) symmetric tridiagonal SY symmetric TB triangular band TG triangular matrices, generalized problem (i.e., a pair of triangular matrices) TP triangular, packed storage TR triangular (or in some cases quasi-triangular) TZ trapezoidal UN (complex) unitary UP (complex) unitary, packed storage ZZZ = indicate the computation performed e.g., SV = Solve, SVX = Solve Expert

Example: list of LAPack Routines for solving linear equations

(See: click here)

Type of matrix and storage scheme	Operation	Single precision		Double precision
		real	complex	real	complex
general	simple driver	SGESV	CGESV	DGESV	ZGESV
	expert driver	SGESVX	CGESVX	DGESVX	ZGESVX
general band	simple driver	SGBSV	CGBSV	DGBSV	ZGBSV
	expert driver	SGBSVX	CGBSVX	DGBSVX	ZGBSVX
general tridiagonal	simple driver	SGTSV	CGTSV	DGTSV	ZGTSV
	expert driver	SGTSVX	CGTSVX	DGTSVX	ZGTSVX
symmetric/Hermitian	simple driver	SPOSV	CPOSV	DPOSV	ZPOSV
positive definite	expert driver	SPOSVX	CPOSVX	DPOSVX	ZPOSVX
symmetric/Hermitian	simple driver	SPPSV	CPPSV	DPPSV	ZPPSV
positive definite (packed storage)	expert driver	SPPSVX	CPPSVX	DPPSVX	ZPPSVX
symmetric/Hermitian	simple driver	SPBSV	CPBSV	DPBSV	ZPBSV
positive definite band	expert driver	SPBSVX	CPBSVX	DPBSVX	ZPBSVX
symmetric/Hermitian	simple driver	SPTSV	CPTSV	DPTSV	ZPTSV
positive definite tridiagonal	expert driver	SPTSVX	CPTSVX	DPTSVX	ZPTSVX
symmetric/Hermitian	simple driver	SSYSV	CHESV	DSYSV	ZHESV
indefinite	expert driver	SSYSVX	CHESVX	DSYSVX	ZHESVX
complex symmetric	simple driver		CSYSV		ZSYSV
	expert driver		CSYSVX		ZSYSVX
symmetric/Hermitian	simple driver	SSPSV	CHPSV	DSPSV	ZHPSV
indefinite (packed storage)	expert driver	SSPSVX	CHPSVX	DSPSVX	ZHPSVX
complex symmetric	simple driver		CSPSV		ZSPSV
(packed storage)	expert driver		CSPSVX		ZSPSVX

Obtaining details on LAPack routines
- Method 1:
  - Use this internet serach engine for LAPack: click here
  - Click on the name of the routine to see its man page entry
- Method 2:
  - Use the LAPack man pages from: click here
  - Uncompress and extract the archive using:
  - Following instructions in README file within the archive
- NOTE: while CS561 is in session, the LAPack man pages are available in:
  Try:

How to use LAPack

Preliminary note:

Before

what the LAPack functions do (i.e., solve the right problem)
the meaning (how to use) of each parameter (i.e., use the right function in the correct manner)

"Learning" about LAPack is not the object of this course.

I leave that to your exploratory endeavor.

But I will show you how to use LAPack so you can use the libray when you find the need for it.

I will illustrate how to use LAPack by using a program that calls SGESV
(I am relying on my limited recollection of numerical analysis that I did as a graduate student in Holland: click here)

Obtaining information on the SGESV function:

nroff -man /home/cs561000/man/manl/sgesv.l | more

Summary:

SUBROUTINE SGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) INTEGER N INTEGER NRHS REAL A(:,:) INTEGER LDA INTEGER IPIV(:) REAL B(:,:) INTEGER LDB INTEGER INFO

N = (input) INTEGER. The number of linear equations, i.e., the order of the matrix A. N >= 0.
NRHS = (input) INTEGER. The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.
A = (input/output) REAL array, dimension (LDA,N).
On entry, it is the N-by-N coefficient matrix A.
On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.
LDA = (input) INTEGER. The leading dimension of the array A. LDA >= solution matrix X. (PS: please don't ask me what "leading dimension" mean :-))
IPIV = (output) INTEGER array, dimension (N). The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i) --- i.e., IPIV gives the pivot rows used to perform the LU decomposition.
B = (input/output) REAL array, dimension (LDB,NRHS)
On entry, it is the N-by-NRHS matrix of right hand side matrix B.
On exit, if INFO = 0, the N-by-NRHS solution matrix X.
LDB = (input) INTEGER The leading dimension of the array B. LDB >= max(1,N).
INFO = (output) INTEGER
INFO = 0: successful exit.
INFO < 0: if INFO = -i, the i-th argument had an illegal value.
INFO > 0: if INFO = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.

Example usage:

PROGRAM Main IMPLICIT NONE REAL :: A(3,3), b(3) INTEGER :: i, j, pivot(3), ok ** Initialize A(:,:) ** Initialize b(:) CALL SGESV(3, 1, A, 3, pivot, b, 3, ok) PRINT *, b !! Print solution END PROGRAM

Example Program: (Demo above code)

Prog file: click here

Compile using the following commands:

f90 -c Solve01.f90
f90 -xlang=f77 Solve01.o -L/home/cs561000/lib -lLAPack -lBLAS
- The flag "-xlang=f77" includes the appropriate runtime libraries and insure the proper runtime environment for legacy Fortran 77....
- I have stored the LAPack libraries in /home/cs561000/lib

Using LA Pack in C-program

Recall that:
See: click here

Therefore, we must transpose the matrix before passing the matrix to a subroutine of LA Pack

Example:

#define SIZE 3 float A[SIZE][SIZE]; // A matrix float AT[SIZE*SIZE]; // A transposed matrix for (i=0; i < SIZE; i++) { for(j=0; j < SIZE; j++) AT[j+SIZE*i]=A[j][i]; }

There is one important caveat that you must be aware of
Let me illustrate the caveat with the subroutine SGESV

Recall the Fortran example :

PROGRAM Main IMPLICIT NONE REAL :: A(3,3), b(3) INTEGER :: i, j, pivot(3), ok ** Initialize A(:,:) ** Initialize b(:) CALL SGESV(3, 1, A, 3, pivot, b, 3, ok) PRINT *, b !! Print solution END PROGRAM

A naive translation into C++ would result in:

int main(int argc, char **argv) { int i, j , pivot[SIZE], ok; float A[SIZE][SIZE], b[SIZE]; float AT[SIZE*SIZE]; // Initialize A(:,:) // Initialize b(:) for (i=0; i < SIZE; i++) // Convert to Fortran storage format { for(j=0; j < SIZE; j++) AT[j+SIZE*i]=A[j][i]; } sgesv_(3, 1, AT, 3, pivot, b, 3, &ok) }

Example Program: (Demo above code)

Prog file: click here

Compile using the following commands:

cc -c Solve02.c

You will get compilation error
What is wrong ???

Recall also that parameters are passed by reference to non-recursive routines in Fortran (all routines in LAPack are non-recursive)
Caveat

Fortran passes constant values by by-reference

Example:

call SGESV(3, 1, A, 3, pivot, b, 3, ok)

passes the value 3 as follows:

integer dummy1 = 3 call SGESV(dummy1, 1, A, 3, pivot, b, 3, ok)

This makes sense if you think throught it a little bit:

The function cannot know if a variable or a constant has been passed as a parameter.

Variables are pass by-reference

The only way to accommodate a constant is to pass it "by-reference" also...

Example using LAPack from C-program:

int main(int argc, char **argv) { int i, j , c1, c2, pivot[size], ok; /* Used to store constants */ float A[size][size], b[size]; /* Input problem */ float AT[size*size]; /* Help variable - stores column major */ // Initialize A() // Initialize b()... /* ********************************* */ /* Convert row-major to column major */ /* ********************************* */ for (i=0; i < size; i++) { for (j=0; j < size; j++) AT[j+size*i]=A[j][i]; } c1=size; /* Put all CONSTANTS we want to pass */ c2=1; /* to the routine in VARIABLES */ /* Find solution using LAPACK routine SGESV. */ /* 1. all the arguments have to be passed by reference */ /* 2. you have to add an underscore to the routine name */ sgesv_ (&c1, &c2, AT, &c1, pivot, b, &c1, &ok); for (j=0; j < size; j++) printf("%f\n", b[j]); /* print solution vector */ }

Notice that you must store a constant in a variable in order to pass the constant value by-reference

Example Program: (Demo above code)

Prog file: click here

Compile using the following commands:

cc -c Solve03.c
f90 -xlang=f77 Solve03.o -L/home/cs561000/lib -lLAPack -lBLAS