LAPACK

LAPACK

General

According to my numerical analysis book, LAPACK has supplanted LINPACK and EISPACK for linear algebra routines, as well as using BLAS as much as possible.

Naming, etc.

Major 'whole-problem' solutions in LAPACK are called 'drivers', while specific linear algebra algorithms are called 'algorithms'. The naming scheme for routines is usually a 6-letter name of the form:

XYYZZZ

where

X - the data type ('S'ingle, 'D'ouble, 'C'omplex or 'Z' - complex*16) YY - the matrix type, ie 'SY'mmetric ZZZ - the routine, for example EVD for divide-and-conquer eigenvalue/vector decomposition

CLAPACK

I couldn't get code calling the LAPACK library routines directly to work, so I downloaded and installed CLAPACK, which is just the fortran code run through f2c, a fortran-to-c compiler. Read the documentation, you'll need to make the BLAS libraries too (this takes a while). A test is on my machine in code/Projects/Utah/test/clapack/dsyev_example.cpp, it is a modified version of this example.

dsyev_example.cpp

#include <cstdlib>
#include <stdio.h>
#include <cmath>

extern "C" {

#include "blaswrap.h"
#include "f2c.h"
#include "clapack.h"

}

#define NDIM 3 

int main (){
  
  int row, col, info2;
  char JOBZ, UPLO;
  long N, LDA, LWORK, INFO;
  double *A;
  double *W;
  double *WORK;

  // symmetric 3x3 matrix
  // will hold eigenvectors upon completion
  N = NDIM;
  A = (double*) malloc(N*N*sizeof(double));
  // ouput array for eigenvalues
  W = (double*) malloc(N*sizeof(double));
  LWORK = (64+2)*N; // size of WORK array
  WORK = (double*) malloc(LWORK*sizeof(double));
  
  // matrix must be in column-major order

  // row 0
  A[0] = 1.0;
  A[3] = -1.0;
  A[6] = 2.0;

  // row 1
  A[1] = 0.0;
  A[4] = 3.0;
  A[7] = -2.0;

  // row 2
  A[2] = 0.0;
  A[5] = 0.0;
  A[8] = -1.0;

  for (row=0;row<N;row++){
    printf("[");
    for (col=0;col<N;col++) {
      printf(" %f ",A[N*col+row]);
    }
    printf("]\n");
  }
  
  JOBZ = 'V'; // compute eigenvectors as well as values
  UPLO = 'U'; // data stored in upper part of matrix
  // N=3, 3x3 matrix
  LDA = N; // Leading dimension of A = 3
  // W is resulting eigenvalues in ascending order
  // WORK is the workspace array
  // LWORK is the size of the WORK array
  // INFO is where error codes are returned
  
  dsyev_(&JOBZ, &UPLO, &N, A, &LDA, W, WORK, &LWORK, &INFO);
  
  printf("info %d \n",INFO);

  // eigenvectors are the columns of A

  printf("--- A ---\n");
  for (row=0;row<N;row++){
    printf("[");
    for (col=0;col<N;col++) {
      printf(" %f ",A[N*col+row]);
    }
    printf("]\n");
  }
  
  printf("--- W ---\n");
  for (row=0;row<N;row++)
    printf("   %f \n",W[row]);
}

And the Makefile:

#       
# Makefile for building a program that calls a routine from
# the CLAPACK Library
#
# Dr. Andrew J. Pounds
# Departments of Chemistry and Computer Science
# Mercer University
#
# November 2004
#

CC = g++
ROOTPATH = /wasabi-local/jsam/src/CLAPACK
INCDIRS = -I$(ROOTPATH)/SRC -I$(ROOTPATH)
F2CDIR  = $(ROOTPATH)/F2CLIBS
LDLIBS  = $(ROOTPATH)/lapack_LINUX.a \
          $(ROOTPATH)/blas_LINUX.a \
          $(F2CDIR)/libF77.a $(F2CDIR)/libI77.a -lm 

dgesv_example: dgesv_example.o
   $(CC) dgesv_example.o  $(LDLIBS)  -o dgesv_example 

dgesv_example.o: dgesv_example.c
   $(CC) dgesv_example.c  $(INCDIRS) -c 

dsyev_example: dsyev_example.o
   $(CC) dsyev_example.o  $(LDLIBS)  -o dsyev_example 

dsyev_example.o: dsyev_example.cpp
   $(CC) dsyev_example.cpp  $(INCDIRS) -c

Other Libraries

I also tried the LAPACK++ and TNT libraries (LAPACK++ says it was supplanted by TNT, but a note in the later versions of LAPACK++ says TNT never took off and development continued on LAPACK++), but LAPACK++ didn't seem to have the routines I needed, and TNT wasn't immediately working for me, though I didn't try too hard. Also, I checked and ATLAS didn't support the library routines I needed from LAPACK.

LAPACK

General

Naming, etc.

CLAPACK

Other Libraries

Links