File Exchange

image thumbnail

rook pivoting

version 1.0 (35.6 KB) by

LURP: Gaussian elimination with rook pivoting using efficient Fortran and MATLAB mex code



View License

The LURP package has efficient Fortran code to compute
LU factorizations of dense, unsymmetric matrices using
Gaussian elimination with rook pivoting. Also the package
includes a MATLAB mex interface so that the routines
can be called directly from MATLAB as well as MATLAB
code to install and test the package.
Gaussian elimination with rook pivoting produces an
LU factorization of a matrix A: PAQ = LU = LDV, where
P and Q are permutation matrices, L is unit lower
triangular or trapezoidal matrix, U is an upper triangular
or trapezoidal matrix, D is a diagonal matrix and V is an
unit upper triangular or trapezoidal matrix. The trapezoidal
matrices appear when A is not square. In practice rook
pivoting produces a rank revealing factorization which can
be used to construct fundamental subspaces of a matrix,
to solve systems of equations involving rank deficient
matrices, for basis repair in optimization algorithms and
for other uses. (See Nicholas Higham's text "Accuracy and
Stability of Numerical Algorithms", 2nd edition, pp.
159-160, 170, 188, 219-220 and its references for additional
The Basic Linear Algebra Subroutines (BLAS) are sets
of linear algebra primitives that are highly optimized
for specific computer architectures. In developing
efficient linear algebra software it is critical to
incorporate the BLAS and, in particular the level three
BLAS (BLAS-3), as much as possible. Code using BLAS-3
routines is more efficient than code using BLAS-1 and
BLAS-2 routines. The Fortran routines supplied with the
LURP package use BLAS-3 routines whenever possible.
(See Higham's text pp. 578-579 for more information
about the BLAS.)
Standard rook pivoting results in a factorization with
the largest magnitude elements in L and in V no larger
than one. The code also includes threshold rook
pivoting which produces factorizations with the largest
magnitude elements in L and in V no larger than a user
selected pivoting tolerance TOL.
Here is comparison of speed of the LURP routine with
alternatives for factoring a 4000 by 4000 random matrix
on 3 computers (A, B and C):
rook pivoting, [L,U]= lurp(A) ------- A: 8.8s, B: 3.3s, C: 3.7s
threshold rook, [L,U] = lurp(A,2) -- A: 7.8s, B: 2.9s, C: 3.1s
partial pivoting, [L,U] = lu(A) ------- A: 5.8s, B: 2.5s, C: 2.0s
QR factorization, [Q,R,P] = qr(A) -- A: 69s, B: 62s, C: 24s
singular value, [U,D,V] = svd(A) --- A: 163s, B: 126s, C: 38s
unblocked rook pivoting code ------ A: 109s, B: 79s, C: 28s
Computer A: 2 cores, Intel SU4100 processor, MATLAB 7.14
Computer B: 4 cores, Intel Xeon E5404 processor, MATLAB 7.14
Computer C: 8 cores, Intel Xeon E5620 processor, MATLAB 7.12
The results show that MATLAB's partial pivoting code is 30%
to 85% faster than rook pivoting for these matrices and
computers. Note that partial pivoting does not reliably
reveal rank. The results also show that rook pivoting
is faster, by more than a factor of 6, than other MATLAB
dense matrix routines that produce rank revealing
factorizations in practice. The unblocked algorithm
uses BLAS-1 and BLAS-2 but not BLAS-3 routines. It is
much slower than the BLAS-3 rook pivoting code.

To install the package:

   download and uncompress the zip file containing the
   start MATLAB and move from inside MATLAB to the folder
       containing the uncompressed files
   set the proper compiler using MATLAB's mex utility: type
       mex -setup
       and following the instructions
Lurp_install requires either gfortran for Linux computers
or Intel Fortran on Windows computers. See lurp_install
("help lurp_install") for more discussion of the

The files in this package include:
  lurp_install.m -- installs the LURP package
  lurp.m -- "help lurp" describes the use of LURP
  lurp.F -- the Mex gateway routine for LURP package
  lurpf.F -- a Mex "computational routine" called by lurp.F
  dgerpf.F -- Fortran code implementing rook pivoting
  dgerp2.F, dgerp3.F, dlaswq.F - called by dgerpf.F
  ilaenv.F -- Fortran version of LAPACK code that returns
     the block size. Included since the Fortran version
     returns a good block size whereas the version of
     ILAENV that comes with MATLAB's LAPACK library may not
     return a good block size.
  lurp_test.m, lurp_test_SJ.m, lurp_test_UF.m -- test routines
  readme.txt -- describes the installation and use of the
The LURP package was created by L. Foster, Department of
Mathematics, San Jose State University. 7-31-2012, copr.

Comments and Ratings (1)

If you do not have a Fortran compiler, compiled binaries -- lurp.mexw32, lurp.mexw64 and lurp.mexa64 -- needed to run the rook factorization routine, lurp, from MATLAB are at Leslie Foster

MATLAB Release
MATLAB 7.12 (R2011a)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video