Re: command \

[Christophe Prud'homme]

Hi Dave 

thank you so much for the detailed info !!

basically we have the same functionalities except that :

o octave (2.1.72, 2.9.4) does not support symbolic calculation
o octave-2.1.72 is not interfaced with umfpack


[ Wednesday 8 March 2006 15:17 ]
| 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

does that apply to some extent also for octave 2.9.4 or only to CVS ?






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

-- 
Christophe Prud'homme
EPFL SB IACS CMCS
MA B2 534 (Bâtiment MA)
Station 8
CH-1015 Lausanne
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
-------------------------------------------------------------