Re: Using LAPACK in octave

[David Bateman]

David Bateman wrote:
> John W. Eaton wrote:
>   
>> On 11-Sep-2007, David Bateman wrote:
>>
>> | David Bateman wrote:
>> | > 
>> | > I've often felt that Octave should change to use the same solver as
>> | > Matlab, especially now as the sparse matrices in Octave use a QR solver
>> | > with pivoting and don't give exactly the same results as the dense 
>> solvers..
>> | 
>> | If you want another reason to change consider Figure 3.3 on the page
>> | 
>> | http://www.netlib.org/lapack/lug/node71.html
>> | 
>> | for N = 500 dgelss is about 35 times slower than dgelsy
>>
>> Switching is fine with me if the results are approximately as good.  I
>> chose dgelss because it was available and I knew about it.
>>     
>
> If we want to stay with an SVD solver we should at least change to
> dgelsd that is significantly faster than dgelss and gives essentially
> the same results.
>
> D.
>
>   
In fact xGELS, xGELSY, xGELSD and xGELSS all return the same minimum
norm solutions. There is therefore no reason not to change to xGELSS..
Consider the attached patch, before this addition I had the following

octave:1> A = rand(500,499); t = cputime(); b = A \ ones(500,1);
cputime() - t
ans =  3.1035

and after the patch I get

 A = rand(500,499); t = cputime(); b = A \ ones(500,1); cputime() - t
ans =  0.23296

with ATLAS 3.6.0 installed. As for the accuracy, both versions returned.

octave:2> norm(A*b)
ans =  22.361

One thing that might be done in addition is to remove the non longer
needed xGELSS code from LAPACK, however care needs to be taken to not
remove a dependency of another function used in Octave, so I decided not
to do that..

Regards
David

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

2007-09-25  David Bateman  <address@hidden>

        * lapack/dgelsy.f,  lapack/dlatrz.f,  lapack/zlarz.f,
        lapack/dgeqp3.f,  lapack/dtzrzf.f,  lapack/zlarzt.f,
        lapack/dlaic1.f,  lapack/zgelsy.f,  lapack/zlatrz.f,
        lapack/dlaqp2.f,  lapack/zgeqp3.f,  lapack/ztzrzf.f,
        lapack/dlaqps.f,  lapack/zlaic1.f,  lapack/zunmr3.f,
        lapack/dlarzb.f,  lapack/zlaqp2.f,  lapack/zunmrz.f,
        lapack/dlarz.f,   lapack/zlaqps.f, lapack/dlarzt.f,  
        lapack/zlarzb.f: New files
        * lapack/Makefile.in (FSRC): Add the new files.

2007-09-25  David Bateman  <address@hidden>

        * dMatrix.cc (lssolve): Replace the use of xGELSS with xGELSY with
        is much faster and no less accurate.
        * CMatrix.cc (lssolve): ditto.

patch8.gz
Description: application/gzip