Re: command \

[David Bateman]

Hi Christophe,

Your description below is pretty much the goal of the octave, though I
think Matlab also special cases at least three other matrices you don't
include, these being

1) If A is symbolic, matlab uses maple for symbolic factorization
2) If A is sparse and diagonal (or permuted diagonal), the result is
just a scaled version of the input
3) If A is dense Hessenberg matrix, reduce to upper triangular and solve.

I haven't gotten to implementing the dense upper/lower triangular and
cholesky solvers yet (this is on my todo list for octave 3.0), but once
that is done the only missing things from matlab relative to octave will
be the symbolic and Hessenberg factorization. That is at the moment
dense matrix are only treated by LU or QR factorization. The current
situation for sparse matrices in the octave CVS tree is the following

  1. If the matrix is diagonal, solve directly and goto 8

  2. If the matrix is a permuted diagonal, solve directly taking into
     account the permutations. Goto 8

  3. If the matrix is square, banded and if the band density is less
     than that given by `spparms ("bandden")' continue, else goto 4.

       a. If the matrix is tridiagonal and the right-hand side is not
          sparse continue, else goto 3b.

            1. If the matrix is hermitian, with a positive real
               diagonal, attempt       Cholesky factorization using
               LAPACK xPTSV.

            2. If the above failed or the matrix is not hermitian with
               a positive       real diagonal use Gaussian elimination
               with pivoting using       LAPACK xGTSV, and goto 8.

       b. If the matrix is hermitian with a positive real diagonal,
          attempt       Cholesky factorization using LAPACK xPBTRF.

       c. if the above failed or the matrix is not hermitian with a
          positive       real diagonal use Gaussian elimination with
          pivoting using       LAPACK xGBTRF, and goto 8.

  4. If the matrix is upper or lower triangular perform a sparse forward
     or backward substitution, and goto 8

  5. If the matrix is a upper triangular matrix with column permutations
     or lower triangular matrix with row permutations, perform a sparse
     forward or backward substitution, and goto 8

  6. If the matrix is square, hermitian with a real positive diagonal,
     attempt sparse Cholesky factorization using CHOLMOD.

  7. If the sparse Cholesky factorization failed or the matrix is not
     hermitian with a real positive diagonal, and the matrix is square,
     factorize using UMFPACK.

  8. If the matrix is not square, or any of the previous solvers flags
     a singular or near singular matrix, find a minimum norm solution
     using CXSPARSE.

Note that I could probably use a dmperm on the permuted triangular
matrix to allow both row and column permutations. I didn't have this
function when I wrote the code though. Also note that at step 8 above
I'm currently using CXSPARSE's QR solver directly, that uses householder
reflections, whereas Matlab uses given's rotations. I'm in the process
of modifying step 8 so that it performs a Dulmange-Mendelsohn (dmperm)
on the matrix and reduces it to 1 leading over-determined problem, N
well determined problems and one trailing under-determined problem, that
should result in significant acceleration and stability over the
existing sparse QR in octave.

Regards
David



Christophe Prud'homme wrote:

> Hello
>
> I am part of the team that writes "Scientific Computing with MATLAB"
> (*). We
> are now working on a new edition that will "work" with both matlab and
> octave
> with mentions of some differences or incompatibilities.
> At some point we describe the command \ fpor Matlab, does it follow
> the same
> logic in octave or is it different ? here is what we have for matlab
>
> -----------------------------------------------------------------------
> It is useful to know that the specific algorithm used by MATLAB when
> the \ command is invoked depends upon the structure of the matrix A.
> To determine the structure of A and select the appropriate algorithm,
> MATLAB follows this precedence:
>
> 1. if A is sparse and banded, then banded solvers are used (like the
> Thomas algorithm of Section 5.4). We say that a matrix A \220 Rm×n
> (or in Cm×n) has lower band p if aij = 0 when i > j + p and upper
> band q if aij = 0 when j > i + q. The maximum between p and q is
> called the bandwidth of the matrix.
>
> 2. If A is an upper or lower triangular matrix (or else a permutation of
> a triangular matrix), then the system is solved by a backsubstitution
> algorithm for upper triangular matrices, or by a forward substitution
> algorithm for lower triangular matrices. The check for triangularity
> is done for full matrices by testing for zero elements and for sparse
> matrices by accessing the sparse data structure.
>
> 3. If A is symmetric and has real positive diagonal elements (which does
> not imply that A is positive definite), then a Cholesky factorization
> is attempted (chol). If A is sparse, a preordering algorithm is applied
> first.
>
> 4. If none of previous criteria are fulfilled, then a general
> triangular fac-
> torization is computed by Gaussian elimination with partial pivoting
> (lu).
>
> 5. If A is sparse, then the UMFPACK library is used to compute the
> solution of the system.
>
> 6. If A is not square, proper methods based on the QR factorization
> for undetermined systems are used (for the overdetermined case, see
> Section 5.5).
>
> (*)http://www.amazon.com/gp/product/3540443630/sr=8-3/qid=1141824041/ref=pd_bbs_3/104-1361323-2155125?%5Fencoding=UTF8
>
> -- 
> Christophe Prud'homme
> Office MA B2 534
> CMCS-IACS EPFL
> CH-1015 Lausanne, Switzerland
> Tel: +41 (0)21 693 25 47
> Fax: +41 (0)21 693 43 03
>
>
>
> -------------------------------------------------------------
> 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 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



-------------------------------------------------------------
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
-------------------------------------------------------------