From help-octave-request at bevo dot che dot wisc dot edu Mon Jan 26 08:33:14 2004 Subject: Re: Understanding how octave works... From: David Bateman To: Joerg Frochte Cc: help-octave at bevo dot che dot wisc dot edu Date: Mon, 26 Jan 2004 15:27:18 +0100 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 (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 David dot Bateman at motorola dot com 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 -------------------------------------------------------------