Re: Understanding how octave works...

[David Bateman]

Consider the code in src/xdiv.cc

Matrix
xdiv (const Matrix& a, const Matrix& b)
{
  if (! mx_div_conform (a, b))
    return Matrix ();

  Matrix atmp = a.transpose ();
  Matrix btmp = b.transpose ();

  int info;
  if (btmp.rows () == btmp.columns ())
    {
      double rcond = 0.0;

      Matrix result
        = btmp.solve (atmp, info, rcond, solve_singularity_warning);

      if (result_ok (info))
        return Matrix (result.transpose ());
    }

  int rank;
  Matrix result = btmp.lssolve (atmp, info, rank);

  return result.transpose ();
}


If the matrix is square we call the "solve" routine based on an LU
decomposition of the matrix (dgetrf and dgetrs). If LAPACK flags the
result as being doubtful with the info flag then we fall back on
"lssolve" that uses an SVD (dgelss).

Maybe you should do the same...

D.



According to Joerg Frochte <address@hidden> (on 01/26/04):
> Hello,
> 
> I often use octave to control results of a program that uses lapack.
> 
> I have got to solve the following :
> 
> A=[
>  0.00e+00  1.00e-01  0.00e+00  5.00e-03  0.00e+00  0.00e+00  0.00e+00
>  0.00e+00  1.67e-04  ;
>  0.00e+00 -1.00e-01  0.00e+00  5.00e-03 -0.00e+00  0.00e+00 -0.00e+00
>  0.00e+00 -1.67e-04  ;
> -7.09e-02  4.48e-02  2.51e-03  1.01e-03 -3.18e-03 -5.94e-05  1.13e-04
> -7.13e-05  1.50e-05  ;
> -7.09e-02 -5.52e-02  2.51e-03  1.52e-03  3.91e-03 -5.94e-05 -1.39e-04
> -1.08e-04 -2.80e-05  ;
> -3.55e-02  7.24e-02  6.29e-04  2.62e-03 -2.57e-03 -7.43e-06  4.55e-05
> -9.30e-05  6.33e-05  ;
>  0.00e+00  5.00e-02  0.00e+00  1.25e-03  0.00e+00  0.00e+00  0.00e+00
>  0.00e+00  2.08e-05  ;
> -3.55e-02 -7.76e-02  6.29e-04  3.01e-03  2.75e-03 -7.43e-06 -4.88e-05
> -1.07e-04 -7.78e-05  ;
>  0.00e+00 -5.00e-02  0.00e+00  1.25e-03 -0.00e+00  0.00e+00 -0.00e+00
>  0.00e+00 -2.08e-05  ;
> -3.55e-02 -2.76e-02  6.29e-04  3.80e-04  9.78e-04 -7.43e-06 -1.73e-05
> -1.35e-05 -3.50e-06  ;
> -3.55e-02  2.24e-02  6.29e-04  2.51e-04 -7.95e-04 -7.43e-06  1.41e-05
> -8.91e-06  1.88e-06  ;
> -7.09e-02 -5.16e-03  2.51e-03  1.33e-05  3.66e-04 -5.94e-05 -1.30e-05
> -9.44e-07 -2.29e-08  ;
> ]
> 
> b=[  0.00e+00 0.00e+00 7.29e-02 8.51e-02 3.46e-02 0.00e+00 4.39e-02 0.00e+00
> 4.10e-02 3.80e-02 7.92e-02 ]
> 
> 
> octave:16> A\b'
> ans =
> 
>    -1.1090e+00
>     3.7483e-04
>    -1.6496e-02
>    -4.5344e-04
>     1.7691e+00
>     1.0526e-03
>     8.8226e-01
>     2.5739e+00
>    -2.6411e-01
> 
> The result in octave is quite OK, but the result my own program computes
> is very unpleasing.
> 
> ans = [ -7.12e+11-3.70e-05-6.03e+13 1.82e-03 1.74e+00-1.70e+15 1.48e-01
> 2.79e+00 6.10e-02 ]
> 
> This has something to do with the fact that A is very close to be singular.
> 
> octave:21> cond(A)
> ans =  2.9251e+18
> 
> Nevertheless, the result of octave is more pleasing.
> How does octave deal with such a situation?
> I have the source-code but I am unable to find a point to start my analysis,
> because I do not know how octave is designed.
> 
> Could you tell me how octave works in this situation or give me a file and
> line where to start in the octave code?
> 
> Thanks very much,
> 
> Joerg Frochte
> 
> 
> 
> 
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
> 
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------

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

The information contained in this communication has been classified as: 

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



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------