Introduction
This document describes a two-level C interface to LAPACK, consisting of a high-level interface and a middle-level interface. The high-level interface handles all workspace memory allocation internally, while the middle-level interface requires the user to provide workspace arrays as in the original FORTRAN interface. Both interfaces provide support for both column-major and row-major matrices. The prototypes for both interfaces, associated macros and type definitions are contained in the header file lapacke.h.
Naming Schemes
The naming scheme for the high-level interface is to take the FORTRAN LAPACK routine name, make it lower case, and add the prefix LAPACKE_. For example, the LAPACK subroutine DGETRF becomes LAPACKE_dgetrf.
The naming scheme for the middle-level interface is to take the FORTRAN LAPACK routine name, make it lower case, then add the prefix LAPACKE_ and the suffix _work. For example, the LAPACK subroutine DGETRF becomes LAPACKE_dgetrf_work.
Complex Types
Complex data types are defined by the macros lapack_complex_float and lapack_complex_double, which represent single precision and double precision complex data types respectively. It is assumed throughout that the real and imaginary components are stored contiguously in memory, with the real component first. The lapack_complex_float and lapack_complex_double macros can be either C99 _Complex types, a C struct defined type, C++ STL complex types, or a custom complex type. See lapacke.h for more details.
Array Arguments
Arrays are passed as pointers, not as a pointer to pointers. All the LAPACKE routines that take one or more 2D arrays as a pointer receive a single extra parameter of type int. This argument must be equal to either LAPACK_ROW_MAJOR or LAPACK_COL_MAJOR which are defined in lapacke.h, specifying whether the arrays are stored in row-major or column-major order. If a routine has multiple array inputs, they must all use the same ordering.
Note that using row-major ordering may require more memory and time than column-major ordering, because the routine must transpose the row-major order to the column-major order required by the underlying LAPACK routine.
Each 2D array argument in a FORTRAN LAPACK routine has an additional argument that specifies its leading dimension. For row-major 2D arrays, elements within a row are assumed to be contiguous and elements from one row to the next are assumed to be a leading dimension apart. For column-major 2D arrays, elements within a column are assumed to be contiguous and elements from one column to the next are assumed to be a leading dimension apart.
Aliasing of Arguments
Unless specified otherwise, only input arguments (that is, scalars passed by values and arrays specified with the const qualifier) may be legally aliased on a call to C interface to LAPACK.
INFO Parameters
The LAPACKE interface functions set their lapack_int return value to the value of the INFO parameter, which contains information such as error and exit conditions. This differs from LAPACK routines, which return this information as a FORTRAN integer parameter.
In LAPACKE, INFO is used exactly as it is in LAPACK. If INFO returns the row or column number of a matrix using 1-based indexing in FORTRAN, the value is not adjusted for zero-based indexing.
NaN Checking
The high-level interface includes an optional, on by default, NaN check on all matrix inputs before calling any LAPACK routine. This option affects all routines. If the inputs contain any NaNs, the input parameter corresponding matrix will be flagged with an INFO parameter error. For example, if the fifth parameter is found to contain a NaN, the function will return with the value -5.
The NaN check, as well as other parameters, can be disabled by defining LAPACK_DISABLE_NAN_CHECK macro in lapacke.h. The middle-level interface does not contain the NaN check.
Integers
Variables with the FORTRAN type integer are converted to lapack_int in LAPACKE. This conforms with modifiable integer type size, especially given ILP64 programming model: re-defining lapack_int as long int (8 bytes) will be enough to support this model, as lapack_int is defined as int (4 bytes) by default, supporting LP64 programming model.
Logicals
FORTRAN logicals are converted to lapack_logical, which is defined as lapack_int.
Memory Management
All memory management is handled by the functions LAPACKE_malloc and LAPACKE_free. This allows users to easily use their own memory manager instead of the default by modifying their definitions in lapacke.h.
This interface should be thread-safe to the extent that these memory management routines and the underlying LAPACK routines are thread-safe.
New Error Codes
Since the high level interface does not use work arrays, error notification is needed in the event of a user running out of memory.
If a work array cannot be allocated, LAPACK_WORK_MEMORY_ERROR is returned by the function; if there was insufficient memory to complete a transposition, LAPACK_TRANSPOSE_MEMORY_ERROR is returned.
Function List
This section list the currently available LAPACK subroutines that are available in the LAPACKE C interface. The LAPACK base names are given below; the corresponding LAPACKE function name is LAPACKE_xbase or LAPACKE_xbase_work where x is the type: s or d for single or double precision real, c or z for single or double precision complex, with base representing the base name. Function prototypes are given in the file lapacke.h. See the LAPACK documentation for detailed information about the routines and their parameters.
Real Functions
The following LAPACK subroutine base names are supported for single precision (s) and double precision (d+), in both the high-level and middle-level interfaces:
bdsdc bdsqr disna gbbrd gbcon gbequ gbequb gbrfs gbrfsx gbsv gbsvx gbsvxx gbtrf gbtrs gebak gebal gebrd gecon geequ geequb gees geesx geev geevx gehrd gejsv gelqf gels gelsd gelss gelsy geqlf geqp3 geqpf geqrf geqrfp gerfs gerfsx gerqf gesdd gesv gesvd gesvj gesvx gesvxx getrf getri getrs ggbak ggbal gges ggesx ggev ggevx ggglm gghrd gglse ggqrf ggrqf ggsvd ggsvp gtcon gtrfs gtsv gtsvx gttrf gttrs hgeqz hsein hseqr opgtr opmtr orgbr orghr orglq orgql orgqr orgrq orgtr ormbr ormhr ormlq ormql ormqr ormrq ormrz ormtr pbcon pbequ pbrfs pbstf pbsv pbsvx pbtrf pbtrs pftrf pftri pftrs pocon poequ poequb porfs porfsx posv posvx posvxx potrf potri potrs ppcon ppequ pprfs ppsv ppsvx pptrf pptri pptrs pstrf ptcon pteqr ptrfs ptsv ptsvx pttrf pttrs sbev sbevd sbevx sbgst sbgv sbgvd sbgvx sbtrd sfrk spcon spev spevd spevx spgst spgv spgvd spgvx sprfs spsv spsvx sptrd sptrf sptri sptrs stebz stedc stegr stein stemr steqr sterf stev stevd stevr stevx sycon syequb syev syevd syevr syevx sygst sygv sygvd sygvx syrfs syrfsx sysv sysvx sysvxx sytrd sytrf sytri sytrs tbcon tbrfs tbtrs tfsm tftri tfttp tfttr tgevc tgexc tgsen tgsja tgsna tgsyl tpcon tprfs tptri tptrs tpttf tpttr trcon trevc trexc trrfs trsen trsna trsyl trtri trtrs trttf trttp tzrzf
Complex Functions
The following LAPACK subroutine base names are supported for complex single precision (c) and complex double precision (z), in both the high-level and middle-level interfaces:
bdsqr gbbrd gbcon gbequ gbequb gbrfs gbrfsx gbsv gbsvx gbsvxx gbtrf gbtrs gebak gebal gebrd gecon geequ geequb gees geesx geev geevx gehrd gelqf gels gelsd gelss gelsy geqlf geqp3 geqpf geqrf geqrfp gerfs gerfsx gerqf gesdd gesv gesvd gesvx gesvxx getrf getri getrs ggbak ggbal gges ggesx ggev ggevx ggglm gghrd gglse ggqrf ggrqf ggsvd ggsvp gtcon gtrfs gtsv gtsvx gttrf gttrs hbev hbevd hbevx hbgst hbgv hbgvd hbgvx hbtrd hecon heequb heev heevd heevr heevx hegst hegv hegvd hegvx herfs herfsx hesv hesvx hesvxx hetrd hetrf hetri hetrs hfrk hgeqz hpcon hpev hpevd hpevx hpgst hpgv hpgvd hpgvx hprfs hpsv hpsvx hptrd hptrf hptri hptrs hsein hseqr pbcon pbequ pbrfs pbstf pbsv pbsvx pbtrf pbtrs pftrf pftri pftrs pocon poequ poequb porfs porfsx posv posvx posvxx potrf potri potrs ppcon ppequ pprfs ppsv ppsvx pptrf pptri pptrs pstrf ptcon pteqr ptrfs ptsv ptsvx pttrf pttrs spcon sprfs spsv spsvx sptrf sptri sptrs stedc stegr stein stemr steqr sycon syequb syrfs syrfsx sysv sysvx sysvxx sytrf sytri sytrs tbcon tbrfs tbtrs tfsm tftri tfttp tfttr tgevc tgexc tgsen tgsja tgsna tgsyl tpcon tprfs tptri tptrs tpttf tpttr trcon trevc trexc trrfs trsen trsna trsyl trtri trtrs trttf trttp tzrzf ungbr unghr unglq ungql ungqr ungrq ungtr unmbr unmhr unmlq unmql unmqr unmrq unmrz unmtr upgtr upmtr
Mixed Precision Functions
The following LAPACK subroutine base names are supported only for double precision (d) and complex double precision (z):
sgesv sposv
Examples
This section contains examples of calling LAPACKE functions from a C program.
Calling DGEQRF
Suppose you wish to call the function DGEQRF, which computes the QR factorization of a double precision real rectangular matrix in LAPACK, and you wish to have the LAPACKE interface handle the necessary work space memory allocation for you.
The base name for this function is geqrf, which is included in the list of real functions above. The LAPACKE function name is then constructed by prepending LAPACK_ followed by d to the base name: LAPACKE_dgeqrf.
We will assume that our matrix is stored in column-major order in the m-by-n array a, which has a leading dimension of lda. The variable declarations should be as follows:
lapack_int m, n, lda, info; double *a, *tau;
The LAPACKE function call is then:
info = LAPACKE_dgeqrf( LAPACK_COL_MAJOR, m, n, a, lda, tau );
Calling CUNGQR
Suppose you wish to call the function CUNGQR, which generates Q from the results of a QR factorization of a single precision complex rectangular matrix, and you wish to provide the required workspace.
The base name for this function is ungqr, which is included in the list of complex functions above. The LAPACKE function name is then constructed by prepending LAPACK_ followed by c to the base name, with the suffix _work to indicate that the user will supply the work space: LAPACKE_cungqr_work.
We will assume again that our matrix is stored in column-major order in the m-by-n array a, which has a leading dimension of lda. From the LAPACK documentation, the work space array work must have a length of at least n; the length of work is given in lwork. The variable declarations should be as follows:
lapack_int m, n, k, lda, lwork, info; lapack_complex_float *a, *tau, *work;
The LAPACKE function call is then:
info = LAPACKE_cungqr_work( LAPACK_COL_MAJOR, m, n, k, a, lda, tau, work, lwork );
Calling DGELS
In this example, we wish solve the least squares problem min_x || B - Ax || for two right-hand sides using the LAPACK routine DGELS. For input we will use the 5-by-3 matrix
( 1 1 1 )
( 2 3 4 )
A = ( 3 5 2 )
( 4 2 5 )
( 5 4 3 )
and the 5-by-2 matrix
( -10 -3 )
( 12 14 )
B = ( 14 12 )
( 16 16 )
( 18 16 )
We will first store the input matrix as a static C two-dimensional array, which is stored in row-major order, and let LAPACKE handle the work space array allocation. The LAPACK base name for this function is gels, and we will use double precision (d), so the LAPACKE function name is LAPACKE_dgels.
thus lda=3 and ldb=2. The output for each right hand side is stored in b as consecutive vectors of length 3. The correct answer for this problem is the 3-by-2 matrix
( 2 1 ) ( 1 1 ) ( 1 2 )
A complete C program for this example is given below. Note that when the arrays are passed to the LAPACK routine, they must be dereferenced, since LAPACK is expecting arrays of type double *, not double **.
/* Calling DGELS using row-major order */ #include <stdio.h> #include <lapacke.h> int main (int argc, const char * argv[]) { double a[5][3] = {1,1,1,2,3,4,3,5,2,4,2,5,5,4,3}; double b[5][2] = {-10,-3,12,14,14,12,16,16,18,16}; lapack_int info,m,n,lda,ldb,nrhs; int i,j; m = 5; n = 3; nrhs = 2; lda = 3; ldb = 2; info = LAPACKE_dgels(LAPACK_ROW_MAJOR,'N',m,n,nrhs,*a,lda,*b,ldb); for(i=0;i<n;i++) { for(j=0;j<nrhs;j++) { printf("%lf ",b[i][j]); } printf("\n"); } return(info); }
Alternatively, we can use column-major ordering for the matrices in this example, as shown below. Here, the matrices are stored as static one-dimensional C arrays. These arrays have a leading dimension that is equal to the number of rows.
/* Calling DGELS using column-major order */ #include <stdio.h> #include <lapacke.h> int main (int argc, const char * argv[]) { double a[5*3] = {1,2,3,4,5,1,3,5,2,4,1,4,2,5,3}; double b[5*2] = {-10,12,14,16,18,-3,14,12,16,16}; lapack_int info,m,n,lda,ldb,nrhs; int i,j; m = 5; n = 3; nrhs = 2; lda = 5; ldb = 5; info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m,n,nrhs,a,lda,b,ldb); for(i=0;i<n;i++) { for(j=0;j<nrhs;j++) { printf("%lf ",b[i+ldb*j]); } printf("\n"); } return(info); }
Calling CGEQRF and the CBLAS
In this example, we will call the LAPACK routine CGEQRF to compute the QR factorization of a matrix. We then call CUNGQR to construct the Q matrix and then use the CBLAS routine CGEMM to compute QH Q - I to check that Q is Hermitian. The error ||QH Q - I|| is printed at the end of the program.
In the first version, given below, we let LAPACKE handle the memory allocation for the workspace internally. The second version, shown below, uses the workspace query facility for both CGEQRF and CUNGQR to obtain the optimal size for the parameter lwork, which we use to allocate our own workspace in the array work.
/* Calling CGEQRF and CUNGQR to compute Q without workspace querying */ #include <stdio.h> #include <stdlib.h> #include <lapacke.h> #include <cblas.h> int main (int argc, const char * argv[]) { lapack_complex_float *a,*tau,*r,one,zero; lapack_int info,m,n,lda; int i,j; float err=0.0; m = 10; n = 5; lda = m; one = lapack_make_complex_float(1.0,0.0); zero= lapack_make_complex_float(0.0,0.0); a = calloc(m*n,sizeof(lapack_complex_float)); r = calloc(n*n,sizeof(lapack_complex_float)); tau = calloc(m,sizeof(lapack_complex_float)); for(j=0;j<n;j++) for(i=0;i<m;i++) a[i+j*m] = lapack_make_complex_float(i+1,j+1); info = LAPACKE_cgeqrf(LAPACK_COL_MAJOR,m,n,a,lda,tau); info = LAPACKE_cungqr(LAPACK_COL_MAJOR,m,n,n,a,lda,tau); for(j=0;j<n;j++) for(i=0;i<n;i++) r[i+j*n]=(i==j)?-one:zero; cblas_cgemm(CblasColMajor,CblasConjTrans,CblasNoTrans, n,n,m,&one,a,lda,a,lda,&one,r,n ); for(i=0;i<n;i++) for(j=0;j<n;j++) err=MAX(err,cabs(r[i+j*n])); printf("error=%e\n",err); free(tau); free(r); free(a); return(info); }
/* Calling CGEQRF and CUNGQR to compute Q with workspace querying */ #include <stdio.h> #include <stdlib.h> #include <lapacke.h> #include <cblas.h> int main (int argc, const char * argv[]) { lapack_complex_float *a,*tau,*r,*work,one,zero,query; lapack_int info,m,n,lda,lwork; int i,j; float err; m = 10; n = 5; lda = m; one = lapack_make_complex_float(1.0,0.0); zero= lapack_make_complex_float(0.0,0.0); a = calloc(m*n,sizeof(lapack_complex_float)); r = calloc(n*n,sizeof(lapack_complex_float)); tau = calloc(m,sizeof(lapack_complex_float)); for(j=0;j<n;j++) for(i=0;i<m;i++) a[i+j*m] = lapack_make_complex_float(i+1,j+1); info = LAPACKE_cgeqrf_work(LAPACK_COL_MAJOR,m,n,a,lda,tau,&query,-1); lwork = (lapack_int)query; info = LAPACKE_cungqr_work(LAPACK_COL_MAJOR,m,n,n,a,lda,tau,&query,-1); lwork = MAX(lwork,(lapack_int)query); work = calloc(lwork,sizeof(lapack_complex_float)); info = LAPACKE_cgeqrf_work(LAPACK_COL_MAJOR,m,n,a,lda,tau,work,lwork); info = LAPACKE_cungqr_work(LAPACK_COL_MAJOR,m,n,n,a,lda,tau,work,lwork); for(j=0;j<n;j++) for(i=0;i<n;i++) r[i+j*n]=(i==j)?-one:zero; cblas_cgemm(CblasColMajor,CblasConjTrans,CblasNoTrans, n,n,m,&one,a,lda,a,lda,&one,r,n); err=0.0; for(i=0;i<n;i++) for(j=0;j<n;j++) err=MAX(err,cabs(r[i+j*n])); printf("error=%e\n",err); free(work); free(tau); free(r); free(a); return(info); }