nwh/lusol

Matlab interface to LUSOL

LUSOL

LUSOL maintains LU factors of a square or rectangular sparse matrix.

This repository provides LUSOL source code and a Matlab interface.

The code is distributed under the terms of the MIT License or the BSD License.

Contents

• gen/: code generation scripts and specification files
• matlab/: Matlab interface code
• src/: LUSOL Fortran code
• LICENSE: [Common Public License][CPL]
• makefile: GNU Make file to build interface

Download and install

Pre-built downloads are available on the Github release page.

Installation simply requires adding the matlab subdirectory to your Matlab path. This may be done with Matlab's addpath function.

If the interface has not been built, please follow the directions below.

Basic usage

In Matlab:

% get L and U factors
[L U P Q] = lusol(A);

See >>> help lusol.

Advanced usage

In Matlab:

% create lusol object
mylu = lusol_obj(A);

% solve with lusol object (ie x = A\b)
x = mylu.solveA(b);

% update factorization to replace a column
mylu.repcol(v,1);

% solve again with updated factorization
x1 = mylu.solveA(b);

See >>> help lusol_obj.

Build

Environment

The build environments as of 2016-01-26 are:

• Fedora 21 & Matlab 2013b
• Mac OS X 10.11 & Matlab 2015b

Building the LUSOL interface in other environments may require modification of makefile and matlab/lusol_build.m.

Requirements

Linux:

• make
• gcc
• gfortran
• Matlab

Mac:

• Xcode for clang and make
• gfortran (possibly via Homebrew)
• Matlab

Notes:

• The matlab binary must be on the system PATH.
• python3 is required to generate the interface code. However, the interface code is pre-generated and included in the repository.
• It is necessary to launch Xcode and accept the license agreement before building the interface.
• The smaller Xcode Command Line Tools package does not appear to work with Matlab 2015b. The full Xcode install is required.

Setup mex

Matlab's mex compiler driver must be configured to use the appropriate C compiler. This can be achieved by executing mex -setup from the Matlab prompt or operating system terminal. On Linux the selected compiler must be the correct version of gcc. On Mac OS X 10.9 the selected compiler must be clang. It is a good idea to match the compiler suggested on the Mathworks supported compilers page. See this note on Matlab compatibility with Xcode 7.

Install gfortran on Mac OS X

1. Install Homebrew
2. $ brew install gcc

Steps

From the base directory:

$ make
$ make matlab

Test:

$ make matlab_test

See <NOTES.md> for example build output.

Notes

The basic requirements to build LUSOL are GNU make, gfortran, a C compiler, and Matlab. The build works with gcc on Linux and clang on Mac OS X. It may be necessary to modify the compiler variables in the makefile (CC, F90C, and F77C) depending on the operating system and environment.

The matlab executable must be on the system path. On Mac OS X with Matlab R2015b this is achieved with:

$ export PATH=/Applications/MATLAB_R2015b.app/bin:$PATH

The makefile may have to be modified on Mac OS X depending on versions of Matlab and gfortran.

Authors

• LUSOL Fortran code: Michael Saunders
• Matlab interface: Nick Henderson

lusol's Issues

Question about Lusol_obj.repcol

I met a problem in testing lusol_obj.repcol().

I replace one column of matrix A with itself, then recalculate sum of elements in the matrix. Result1 and result2 are always same. But sometimes result2 becomes NaN. When scale of A is bigger as 1000, result will never be correct.

Here is my code in Matlab:

n = 20;

A = sprand(n,n,0.2)/10+speye(n);

lu_options = lusol_obj.luset('nzinit',200000,'pivot','TCP');

one_vector = ones(n,1);
mylu = lusol_obj(A,lu_options);
result1 = one_vector'*mylu.mulA(one_vector)

mylu.repcol(A(:,10),10);   %replace col 10 with himself
result2 = one_vector'*mylu.mulA(one_vector)

result:

test2

result1 =

54.5849

result2 =

54.5849

test2

result1 =

53.1647

result2 =

8.4224e+252

Version of MATLAB, and make matlab_test:

                       < M A T L A B (R) >
              Copyright 1984-2012 The MathWorks, Inc.
                 R2012b (8.0.0.783) 64-bit (maci64)
                          August 22, 2012

test_addcol_lhr02: passed
test_addrow_lhr02: passed
test_delcol_lhr02: passed
test_delrow_lhr02: passed
test_factorize_lhr02: passed
test_mulA_lhr02: passed
test_mulAt_lhr02: passed
test_r1mod_lhr02: passed
test_repcol_lhr02: passed
test_reprow_lhr02: passed
test_solveA_lhr02: passed
test_solveAt_lhr02: passed

and I also tested it in ubuntu 13.04, matlab R2012b. Same problem exist.

$ gfortran --version
GNU Fortran (GCC) 4.8.1
Copyright (C) 2013 Free Software Foundation, Inc.

Rank Updates Incorrect for Matrices that are Not Fullrank

When the matrix being factored was not full rank, it seemed the updates were not being fully applied.

It looks like the problem comes from the following code in the lu8mod function:

  120 if (nrank .lt. m) then
         nrank  = nrank + 1
         jelm   = 0
         call lu7elm( m, n, jelm, v,
     $                lena, luparm, parmlu,
     $                lenl, lenu, lrow, nrank,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
         if (inform .eq. 7) go to 970
         if (inform .eq. 0) nrank = nrank - 1
      end if

I don't think nrank should be incremented before the call to lu7elm. lu7elm already increments nrank. Incrementing it beforehand causes nrank to be incremented twice which causes lu7elm to miss an element when performing the elimination.

I believe the problem can be fixed by changing the above code to

  120 if (nrank .lt. m) then
         jelm   = 0
         call lu7elm( m, n, jelm, v,
     $                lena, luparm, parmlu,
     $                lenl, lenu, lrow, nrank,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
         if (inform .eq. 7) go to 970
         if (inform .eq. 1) nrank = nrank + 1
      end if 

So that nrank is updated after the call to lu7elm.

generation of non upper triangular U

I've come across a problem in LUSOL where the U matrix is not upper triangular.
This can happen when nzinit is too low, increasing it will generate correct results.
The inform flag is always 1 (successful factorization but rank deficient) and not 7 (LUSOL needs more storage).
Here is a script that should reproduce the problem:

m = 800;
n = 900;
A = sparse(m, n);
nfill = 10;
for i = 1:m
    for j = 1:nfill
        c = i + j - nfill / 2;
        if (c >= 1 && c < n) 
            A(i, c) = i;
        end
    end
end

options = lusol_obj.luset('pivot','TRP',...
                          'maxcol', 5,...
                          'keepLU', 1,...
                          'Ltol1', 1.5,...
                          'Ltol2',1.5,...
                          'nzinit', 10000);
mylu = lusol_obj(A, options);

fprintf('LUSOL inform: %i\n', mylu.inform());
if mylu.inform() == 0
    fprintf('factorization successful (full rank)\n');
elseif mylu.inform() == 1
    fprintf('factorization successful (rank deficient)\n');
else
    fprintf('factorization FAILED\n');
end

[U, P, Q] = mylu.U('matrix');

fprintf('min lena: %i\n', mylu.stats().minlen);

fprintf('is U upper triangular: ');
if istriu(U)
    fprintf('yes\n');
else
    fprintf('no\n');
end

The output is

LUSOL inform: 1
factorization successful (rank deficient)
min lena: 11390
is U upper triangular: no

Changing nzinit to 100000 will generate the output:

LUSOL inform: 1
factorization successful (rank deficient)
min lena: 11390
is U upper triangular: yes

How could I get updated L factors using lusol??

Hello, I'm trying to use lusol to solve my problem. In my problem, I need measure the number of non-zeros of L factor after some updates. lusol provides a function to get L0, which is the initial L factor, but I couldn't find a function what I want. How could I get updated L factors using lusol??

And, Is it possible to use only updating functions for L, U factors which are computed by MATLAB?

factorization fails for rook pivoting with zero columns

Factorization using rook pivoting provides a bad U matrix when there are columns with zeros.
Here is a simple example:

A = sparse([1,0,0,0,0;...
            0,0,0,1,0;...
            0,0,0,0,1]);

options = lusol_obj.luset('pivot', 'TRP');
mylu = lusol_obj(A, options);
[U, P, Q] = mylu.U('matrix');
[L0] = mylu.L0('matrix');

fprintf('U matrix: \n'); disp(full(U));
fprintf('L0 matrix: \n'); disp(full(L0));

In my system (Ubuntu 14.04, 64bit, gfortran 4.8.2) I get the following output:

U matrix: 
     1     0     0     0     0
     0     0     0     0     0
     0     0     1     0     0

L0 matrix: 
     1     0     0
     0     1     0
     0     0     1

the second row of U is not OK.

Changing pivoting to partial pivoting will produce the correct results.
Also all tests in lusol_test.m pass.

No changes were performed to the source code.
The only change was the path to matlab in the makefile.

Calling delrow when there are linear dependencies

When using addrow causes some rows to be linearly dependent rows, the information related to one row seems to be lost and, hence, delrow will not provide the expected results.
Here is a simple script to exemplify the issue:

A = sparse([1,2,3,4,0,0;...
            5,6,0,0,0,0;...
            0,0,0,0,7,8]);
[m, n] = size(A);

% factorize A
options = lusol_obj.luset('pivot','TRP',...
                          'maxcol', 5,...
                          'keepLU', 1,...
                          'Ltol1', 2.0,...
                          'Ltol2',2.0);

mylu = lusol_obj(A, options);
assert(mylu.inform() == 1); % 1: rank defecient

% add [5, 0, 0, 0, 0, 0] 
x = zeros(1,n);
x(1) = 5;
mylu.addrow(x);
[U, P, Q] = mylu.U('matrix'); fprintf('U = \n'); disp(full(U));

% add [0, 6, 0, 0, 0, 0]
x = zeros(1,n);
x(2) = 6;
mylu.addrow(x);
[U, P, Q] = mylu.U('matrix'); fprintf('U = \n'); disp(full(U));

% remove [5, 0, 0, 0, 0, 0] 
mylu.delrow(4);
[U, P, Q] = mylu.U('matrix'); fprintf('U = \n'); disp(full(U));

The result is:

U = 
     3     0     2     1     0     4
     0     7     0     0     8     0
     0     0     6     5     0     0
     0     0     0     5     0     0

U = 
     3     0     2     1     0     4
     0     7     0     0     8     0
     0     0     6     0     0     0
     0     0     0     5     0     0
     0     0     0     0     0     0

U = 
     3     0     2     1     0     4
     0     7     0     0     8     0
     0     0     6     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     0

How should one handle this situation (considering that one might not know which rows are linearly dependent)?

Build error

How to add new function into lusol

As the first step of creating a new function in lusol, I copied everything from "lu1fac" and rename it to "lu1cop" in lusol.f90. I claim the function in clusol.c like this:

void __lusol_MOD_lu1cop(
int64_t* m,
int64_t* n,
int64_t* nelem,
int64_t* lena,
int64_t* luparm,
double* parmlu,
double* a,
int64_t* indc,
int64_t* indr,
int64_t* p,
int64_t* q,
int64_t* lenc,
int64_t* lenr,
int64_t* locc,
int64_t* locr,
int64_t* iploc,
int64_t* iqloc,
int64_t* ipinv,
int64_t* iqinv,
double* w,
int64_t* inform);

void clu1cop(
int64_t* m,
int64_t* n,
int64_t* nelem,
int64_t* lena,
int64_t* luparm,
double* parmlu,
double* a,
int64_t* indc,
int64_t* indr,
int64_t* p,
int64_t* q,
int64_t* lenc,
int64_t* lenr,
int64_t* locc,
int64_t* locr,
int64_t* iploc,
int64_t* iqloc,
int64_t* ipinv,
int64_t* iqinv,
double* w,
int64_t* inform) {
__lusol_MOD_lu1cop(m,n,nelem,lena,luparm,parmlu,a,indc,indr,p,q,lenc,lenr,locc,locr,iploc,iqloc,ipinv,iqinv,w,inform);
}

And I added similar thing to clusol.h.

void clu1cop(
int64_t* m,
int64_t* n,
int64_t* nelem,
int64_t* lena,
int64_t* luparm,
double* parmlu,
double* a,
int64_t* indc,
int64_t* indr,
int64_t* p,
int64_t* q,
int64_t* lenc,
int64_t* lenr,
int64_t* locc,
int64_t* locr,
int64_t* iploc,
int64_t* iqloc,
int64_t* ipinv,
int64_t* iqinv,
double* w,
int64_t* inform);

In makefile:

INTERFACE_FILES :=
gen/interface.py
gen/interface_files.org
gen/lu1fac.org
gen/lu1cop.org \

Then, when it is compiled, system returns that

Undefined symbols for architecture x86_64:
"___lusol_MOD_lu1cop", referenced from:
_clu1cop in clusol.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [src/libclusol.dylib] Error 1

This error happens in Mac OS. And it makes me confused that it could pass in linux. Is it problem of compiler?

Use of uninitialized value

As noted in #12, lu8mod consults a variable named jrep on line 759 while deciding if a matrix is singular. There is no parameter or variable named jrep in lu8mod, so an uninitialized value is used. The upshot is that, if execution reaches line 759, the decision as to whether the matrix is singular or not is made arbitrarily, independent of whether the matrix really is singular or not. If building with -Wall, this warning is produced:

src/lusol8b.f:761:21: warning: 'jrep' may be used uninitialized [-Wmaybe-uninitialized]
  761 |             if (.not. singlr) then
      |                     ^
src/lusol8b.f:759:33: note: 'jrep' was declared here
  759 |             singlr = j1 .ne. jrep
      |                                 ^