Lapack
From charlesreid1
Contents
Calling Lapack functions
Calling from C
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
extern void sgetri_ ( int* n, const void* A, int* lda,
int* ipiv, const void* work, int* lwork, int* info );
extern void sgetrf_ ( int* m, int* n, const void* a, int* lda,
int* ipiv, int* info );
int main (void) {
// A matrix (the one being inverted)
float A[] = {1,3,2,4}; // Matrix [1 2
// 3 4]
// (matrices have to be provided to Lapack as 1-D arrays)
int i,j,n;
int lda, arows, acols;
int ipiv[4];
int info;
// Set the dimension of A
n=2;
arows=n;
acols=n;
// Set the leading dimension of A
lda=2;
// Set the dimension of the "workspace" array WORK
// see http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#workspace
int lwork=100;
float work[lwork];
printf ("Matrix A before any calls:\n");
for (i=0;i<n;i++) {
for (j=0;j<n;j++)
printf ("A[%d][%d]:%f\t", i,j, A[n*j+i]); // this looks backwards, but it's b
printf ("\n");
}
printf ("\n");
// Allocate memory for arrays:
// Pivot indices array
memset (ipiv, 0, sizeof(int)*4);
// Work (workspace) array
memset (work, 0, sizeof(float)*100);
// -------------------------------------------
// Calculate the LU decomposition of A (and store it in A)
// see http://www.netlib.org/lapack/single/sgetrf.f
sgetrf_ ( &arows, &acols, A, &lda, ipiv, &info );
printf ("Matrix A after sgetrf\n");
for (i=0;i<2;i++) {
for (j=0;j<2;j++) {
printf("A[%d][%d]:%f\t", i,j, A[n*j+i]); // this looks backwards, but it's be
}
printf("\n");
}
if( info != 0 ) {
printf("sgetrf_ returned an error! \n");
}
printf ("\n");
// -------------------------------------------
// Calculate the inverse of A using the LU decomposition of A
// http://www.netlib.org/lapack/single/sgetri.f
sgetri_ (&n, A, &lda, ipiv, work, &lwork, &info);
printf ("Matrix A after sgetri\n");
for (i=0;i<2;i++) {
for (j=0;j<2;j++) {
printf("A[%d][%d]:%f\t", i,j, A[n*j+i]); // this looks backwards, but it's be
}
printf("\n");
}
if( info != 0 ) {
printf("sgetri_ returned an error! \n");
}
printf ("\n");
return(0);
}
Calling from C++
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
extern "C" void sgetri_ ( int* n, const void* A, int* lda,
int* ipiv, const void* work, int* lwork, int* info );
extern "C" void sgetrf_ ( int* m, int* n, const void* a, int* lda,
int* ipiv, int* info );
int main (void) {
// A matrix (the one being inverted)
float A[] = {1,3,2,4}; // Matrix [1 2
// 3 4]
// (matrices have to be provided to Lapack as 1-D arrays)
int i,j,n;
int lda, arows, acols;
int ipiv[4];
int info;
// Set the dimension of A
n=2;
arows=n;
acols=n;
// Set the leading dimension of A
lda=2;
// Set the dimension of the "workspace" array WORK
// see http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#workspace
int lwork=100;
float work[lwork];
printf ("Matrix A before any calls:\n");
for (i=0;i<n;i++) {
for (j=0;j<n;j++)
printf ("A[%d][%d]:%f\t", i,j, A[n*j+i]); // this looks backwards, but it's b
printf ("\n");
}
printf ("\n");
// Allocate memory for arrays:
// Pivot indices array
memset (ipiv, 0, sizeof(int)*4);
// Work (workspace) array
memset (work, 0, sizeof(float)*100);
// -------------------------------------------
// Calculate the LU decomposition of A (and store it in A)
// see http://www.netlib.org/lapack/single/sgetrf.f
sgetrf_ ( &arows, &acols, A, &lda, ipiv, &info );
printf ("Matrix A after sgetrf\n");
for (i=0;i<2;i++) {
for (j=0;j<2;j++) {
printf("A[%d][%d]:%f\t", i,j, A[n*j+i]); // this looks backwards, but it's be
}
printf("\n");
}
if( info != 0 ) {
printf("sgetrf_ returned an error! \n");
}
printf ("\n");
// -------------------------------------------
// Calculate the inverse of A using the LU decomposition of A
// http://www.netlib.org/lapack/single/sgetri.f
sgetri_ (&n, A, &lda, ipiv, work, &lwork, &info);
printf ("Matrix A after sgetri\n");
for (i=0;i<2;i++) {
for (j=0;j<2;j++) {
printf("A[%d][%d]:%f\t", i,j, A[n*j+i]); // this looks backwards, but it's be
}
printf("\n");
}
if( info != 0 ) {
printf("sgetri_ returned an error! \n");
}
printf ("\n");
return(0);
};
Compiling with Lapack
Compiling C/C++
C and C++ code using Lapack can be compiled with any standard C/C++ compiler; in this example, gcc and g++.
Linux
On Linux, C code using Lapack can be compiled by linking to Lapack libraries. In most cases, the Lapack libraries are located in /usr/lib:
$ ls /usr/lib | grep lapack liblapack-3.a liblapack-3.so liblapack.a liblapack.so
To link to the Lapack libraries, the -l and -L flags must be used:
$ gcc lapack.c -llapack -L/usr/lib -o bin.x
For more details on libraries, see C++, or the Summer Workshop on Scientific Computing (Number 2, Software) at the Presentations page.
Likewise, C++ code can be compiled using:
$ g++ lapack.cpp -llapack -L/usr/lib -o bin.x
Mac
Using Lapack on a Mac is more straightforward, since Macs come pre-built with an ATLAS version of Lapack, Blas, etc. These are packaged in a vector library framework. The compiler can be pointed to these using:
$ gcc lapack.c -framework vecLib -o bin.x
Similarly, C++ code can be compiled with the same -framework compiler flag:
$ g++ lapack.cpp -framework vecLib -o bin.x
Useful Lapack functions
http://www.netlib.org/lapack/double/
Matrix inversion
Eigenvalues
External References
- Using Lapack from C: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=171
- Using Lapack from C++: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=2058
 |
Scientific Computing Topics in scientific computing.
Numerical Software: Lapack · Sundials · Matlab · Octave · FFTW Petsc · Example Petsc Makefile · Trilinos · Hypre · Ginac · Gnuplot
Python: Numpy · Scipy · Pandas · Matplotlib · Python Sundials · Py4Sci Scikit-learn: Sklearn · Skimage
|
|
Linear Algebra Topics in linear algebra.
Matlab · Octave · Sundials · Trilinos
|