Vocabulary/LAPACK
Pre-configuring
The math/lapack2 addon uses a Fortran 77 reference implementation of the BLAS (Basic Linear Algebra Subprograms). It was written to de ne the basic operations and do not employ the standard tricks for optimizing Fortran code, so it will not perform as well as a specifically tuned implementation. Therefore its use is discouraged (see LAWNS 41, 81).
However, it's possible to force addon to use another BLAS implementation. All you need is to set a global noun pointing to your favorite library before loading addon:
liblapack_jlapack2_=: 'c:/j9.6/addons/math/lapack2/lib/libopenblas.dll' NB. use OpenBLAS; replace string by 'libopenblas.so.0' for OS Linux
If not has been set, then it will be initialized by default as:
liblapack_jlapack2_=: 'c:/j9.6/addons/math/lapack2/lib/liblapack3.dll' NB. use reference LAPACK; replace string by 'liblapack.so.3' for OS Linux
Loading
First, go to Package Manager and make sure you have the latest version of the math/lapack2 addon.
To use it:
load 'math/lapack2'
Identifying a LAPACK version
ver_jlapack2_ ''
OpenBLAS 0.3.26 DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=64
NB. to find out the version of built-in LAPACK library
ilaver_jlapack2_ (,0);(,0);(,0)
+-+-+--+-+
|0|3|12|0|
+-+-+--+-+
So, LAPACK addon uses OpenBLAS library of version 0.3.26 with built-in LAPACK library of version 3.12.0.
Finding a LAPACK function
The LAPACK documentation is here. Browse through it to find the function you need, and go to the manual page to find the meaning of the arguments.
Examples of use
QR decomposition
You have a matrix a and you want to find the QR decomposition.
]a=: 3 3 $12. _51 4 6 167 _68 _4 24 _41 NB. Floating-point matrix
12 _51 4
6 167 _68
_4 24 _41
Decide which LAPACK function to use
Go through the LAPACK library and decide which function you want to use. This may take some time, as there are often several different routines depending on what you know about your inputs (are they symmetric, upper-triangular, etc?) and what form of output you want. Consult whatever references you need.
For QR factorization we find GEQRF, listed in some references as ?GEQRF. This needs one more letter to complete the name of the function.
The first letter of a LAPACK function indicates the argument precision, the rest indicate the operation. J floating-point values are double-precision floating-point, so they begin with the letter D.
The routine we will use is DGEQRF in LAPACK. The j verb is named dgeqrf.
Understand the argument layout
Starting from the page at http://www.netlib.org/lapack/explore-html/modules.html, we expand Modules>LAPACK>General Matrices>Computational Routines>double>dgeqrf to get to the page on [dgeqrf]. Read and understand the parameters.
The parameters match those given in the J function header created by the math/lapack2 addon:
load 'math/lapack2'
dgeqrf_jlapack2_
'"c:/j9.6/addons/math/lapack2/lib/liblapack3.dll" dgeqrf_ + n &i &i *d &i *d *d &i *i'&cd
The part of this you need to understand is the argument descriptors n &i &i *d &i *d *d &i *i. (For the full description of the last line see here.) The first field describes the result; the value n indicates that the function does not return a result. The other 8 fields correspond to the definition on the dgeqrf page. i signifies an integer, d floating-point.
Note that every argument is passed by a pointer, indicated by either the * or & in each argument descriptor. This is the FORTRAN call-by-reference standard. What this means to you is that every argument whose descriptor contains either * or & must be an array, not an atom.
Also be aware that the LAPACK routines, being FORTRAN routines, expect matrices to be stored in column-major order rather than J's row-major order. Thus, you will need to transpose arrays going into and coming out of LAPACK routines unless you know they are symmetric.
Write the interface routine
We will build up the routine for QR factorization step by step in an explicit definition.
Reading the argument description, we see that the argument WORK is an array whose size must be determined by a preliminary call to dgeqrf. We'll start with that call:
require 'math/lapack2'
NB. QR factorization of y
NB. Result is Q;R such that Q is unitary, R is upper-triangular, and y -: Q +/ . * R
dgeqrf =: 3 : 0
assert. 2 = #@$ y NB. y must be a matrix
'r c' =. ,"0 $ y NB. get # rows and columns, as lists
NB. preliminary call to find best LWORK
dgeqrf_jlapack2_ c;r;(|:y);(1>.c);((r<.c)$0.);(1$0.);(,_1);,_1
)
dgeqrf a
+-+-+-+-----------+-+-----+--+--+-+
|0|3|3| 12 6 _4|3|0 0 0|96|_1|0|
| | | |_51 167 24| | | | | |
| | | | 4 _68 _41| | | | | |
+-+-+-+-----------+-+-----+--+--+-+
It worked. Good return (the 0 at the end), and the WORK argument has been filled in with the value to use for LWORK. We will extend our program to use that value of LWORK to do the work:
If you pass a name into a LAPACK routine whose descriptor contains *, J makes a copy of the name before calling LAPACK. Thus, the name itself is unchanged, even if the LAPACK documentation says that the field is modified. To see the result of the LAPACK function you must look in the boxed result of the call (ret here).
dgeqrf =: 3 : 0
assert. 2 = #@$ y NB. y must be a matrix
'r c' =. ,"0 $ y NB. get # rows and columns, as lists
NB. preliminary call to find best LWORK
ret =. dgeqrf_jlapack2_ c;r;(|:y);(1>.c);((r<.c)$0.);(1$0.);(,_1);,_1
lwork =. <. (6;0) {:: ret NB. extract best value for LWORK
dgeqrf_jlapack2_ c;r;(|:y);(1>.c);((r<.c)$0.);(lwork$0.);(,lwork);,_1
)
dgeqrf a
+-+-+-+----------------------+-+-----------------+----------------+--+-+
|0|3|3|_14 0.230769 _0.153846|3|1.85714 1.99385 0|3 _5.38462 0 ...|96|0|
| | | |_21 _175 0.0555556| | | | | |
| | | | 14 70 _35| | | | | |
+-+-+-+----------------------+-+-----------------+----------------+--+-+
We have the result of dgeqrf, but where are Q and R? After examination of the description of dgeqrf we realize that Q and R are packed into the result matrix, with the main diagonal and above being R and the part below the diagonal representing a sequence of Householder transformations that identify Q. To realize Q itself we need another search through our resources, and we finally land on the function dorgqr]. It takes the result of dgeqrf and expands the product to produce Q. This makes the whole program:
dgeqrf =: 3 : 0
assert. 2 = #@$ y NB. y must be a matrix
'r c' =. ,"0 $ y NB. get # rows and columns, as lists
NB. preliminary call to find best LWORK
ret =. dgeqrf_jlapack2_ c;r;(|:y);(1>.c);((r<.c)$0.);(1$0.);(,_1);,_1
lwork =. <. (6;0) {:: ret NB. extract recommended size of workarea
ret =. dgeqrf_jlapack2_ c;r;(|:y);(1>.c);((r<.c)$0.);(lwork$0.);(,lwork);,_1
R =. |: ({."0 1~ #\) 3 {:: ret NB. Remember transposed order! take up through the main diagonal of each row, then transpose
ret =. dorgqr_jlapack2_ (1 2 2 3 4 5 6 7{ret),<,_1 NB. Realize Q. Reuse workareas
(|: 4{::ret);R NB. Return Q and R, remembering to transpose the Q result
)
dgeqrf a
+------------------------------+------------+
|_0.857143 0.394286 0.331429|_14 _21 14|
|_0.428571 _0.902857 _0.0342857| 0 _175 70|
| 0.285714 _0.171429 0.942857| 0 0 _35|
+------------------------------+------------+
+/ . *&>/ dgeqrf a
12 _51 4
6 167 _68
_4 24 _41
Estimate condition number of a symmetric matrix
Searching resources we find the LAPACK call DSYCON which requires first calculating the 1-norm of the matrix and then factoring it using DSYTRF. We can calculate the 1-norm in J directly. The program is
require 'math/lapack2'
NB. Program to estimate the reciprocal condition number (in the 1-norm) of symmetric real matrix y.
NB. Result is reciprocal condition number
dsycon =: 3 : 0
assert. 2 = #@$ y
n =. ,#y NB. # rows, as vector
ret =. dsytrf_jlapack2_ (,'U');n;y;n;(n$00);(1$0.);(,_1);,_1 NB. no need to transpose symmetric matrix
lwork =. <. (6;0) {:: ret NB. extract size of workarea
ret =. dsytrf_jlapack2_ (,'U');n;y;n;(n$00);(lwork$0.);(,lwork);,_1 NB. Call again with workarea allocated
ret =. dsycon_jlapack2_ (1 2 3 4 5{ret),(, >./ +/"1 | y);(1$0.);((2*n)$0.);(n$01);,_1
7{::ret
)
dsycon 2 2 $ 2 0 0 1
0.5
Verbs defined in math/lapack2
The lapack2 addon consists mainly of a set of cover verbs for the LAPACK functions. These verbs call the LAPACK DLL. You need to understand how to pass parameters to a DLL. The cover verbs are already set up to give your argument the correct datatype.
All these verbs are defined in the jlapack2 locale.