Boost logo

Ublas :

From: Karl Meerbergen (Karl.Meerbergen_at_[hidden])
Date: 2006-04-14 02:39:36


Hi,

Usually, factorization methods do not crash, but deliver a complete
factorization. The singularity is usually pivoted to the bottom of the
matrix. Some information is passed on to the user. See the manual for
LAPACK that I added below.

Karl

      *SYNOPSIS*

     SUBROUTINE DGETRF(M, N, A, LDA, IPIVOT, INFO)

     INTEGER M, N, LDA, INFO
     INTEGER IPIVOT(*)
     DOUBLE PRECISION A(LDA,*)

     SUBROUTINE DGETRF_64(M, N, A, LDA, IPIVOT, INFO)

     INTEGER*8 M, N, LDA, INFO
     INTEGER*8 IPIVOT(*)
     DOUBLE PRECISION A(LDA,*)

      *PURPOSE*

     dgetrf computes an LU factorization of a general M-by-N
     matrix A using partial pivoting with row interchanges.

     The factorization has the form
        A = P * L * U
     where P is a permutation matrix, L is lower triangular with
     unit diagonal elements (lower trapezoidal if m > n), and U
     is upper triangular (upper trapezoidal if m < n).

     This is the right-looking Level 3 BLAS version of the algo-
     rithm.

      *ARGUMENTS*

     M (input) The number of rows of the matrix A. M >= 0.

     N (input) The number of columns of the matrix A. N >= 0.

     A (input/output)
               On entry, the M-by-N matrix to be factored. On
               exit, the factors L and U from the factorization A
               = P*L*U; the unit diagonal elements of L are not
               stored.

     LDA (input)
               The leading dimension of the array A. LDA >=
               max(1,M).

     IPIVOT (output)
               The pivot indices; for 1 <= i <= min(M,N), row i
               of the matrix was interchanged with row IPIVOT(i).

     INFO (output)
               = 0: successful exit
               < 0: if INFO = -i, the i-th argument had an ille-
               gal value
> 0: if INFO = i, U(i,i) is exactly zero. The
               factorization has been completed, but the factor U
               is exactly singular, and division by zero will
               occur if it is used to solve a system of equa-
               tions.

Sohail Somani wrote:

>Im sure it can, but perhaps the author assumed you knew whether or not
>the matrix was singular (I tend not to assume this)
>
>
>
>>-----Original Message-----
>>From: ublas-bounces_at_[hidden]
>>[mailto:ublas-bounces_at_[hidden]] On Behalf Of Manoj Rajagopalan
>>Sent: Thursday, April 13, 2006 12:53 PM
>>To: ublas mailing list
>>Subject: Re: [ublas] Solving linear equations of the form Ax = b ?
>>
>>I can confirm that this code will crash on singular matrices
>>:-) Can't
>>ublas be modified to throw an exception or return an error value like
>>INFO in the original FORTRAN version?
>>
>>
>>Sohail Somani wrote:
>>
>>
>>>Warning: I think the wiki says that this code will crash on
>>>
>>>
>>a singular
>>
>>
>>>matrix.
>>>
>>>Just fyi.
>>>
>>>
>>>
>>>
>>>>-----Original Message-----
>>>>From: ublas-bounces_at_[hidden]
>>>>[mailto:ublas-bounces_at_[hidden]] On Behalf Of Manoj
>>>>
>>>>
>>Rajagopalan
>>
>>
>>>>Sent: Thursday, April 13, 2006 7:06 AM
>>>>To: ublas mailing list
>>>>Subject: Re: [ublas] Solving linear equations of the form Ax = b ?
>>>>
>>>>you can solve Ax=b using three lines of ublas code:
>>>>
>>>>permutation_matrix<> piv;
>>>>lu_factorize(A, piv);
>>>>lu_substitute(A, piv, x);
>>>>
>>>>Then x will contain the solution and A is overwritten with its LU
>>>>decomposition. i.e., the original values of A and x are
>>>>destroyed so the
>>>>above lines are actually the inplace_solve(A, x) for a dense,
>>>>generic A.
>>>>
>>>>BTW, I've come across this in many places - what does
>>>>"boilerplate code"
>>>>mean and why is it called so?
>>>>
>>>>cheers!
>>>>Manoj
>>>>
>>>>
>>>>
>>>>John Maddock wrote:
>>>>
>>>>
>>>>
>>>>>I gather from a quick web search that ublas does have LUP
>>>>>
>>>>>
>>>>decomposition
>>>>
>>>>
>>>>
>>>>>routines, however a grep of the docs failed to find any
>>>>>
>>>>>
>>>>mention of the
>>>>
>>>>
>>>>
>>>>>routine in question, and how they should be used.
>>>>>
>>>>>So.... does anyone have the boiler plate code needed to
>>>>>
>>>>>
>>>>solve for x given
>>>>
>>>>
>>>>
>>>>>matrix A and vector b ?
>>>>>
>>>>>Many thanks,
>>>>>
>>>>>John Maddock.
>>>>>
>>>>>_______________________________________________
>>>>>ublas mailing list
>>>>>ublas_at_[hidden]
>>>>>http://lists.boost.org/mailman/listinfo.cgi/ublas
>>>>>
>>>>>
>>>>_______________________________________________
>>>>ublas mailing list
>>>>ublas_at_[hidden]
>>>>http://lists.boost.org/mailman/listinfo.cgi/ublas
>>>>
>>>>
>>>>
>>>_______________________________________________
>>>ublas mailing list
>>>ublas_at_[hidden]
>>>http://lists.boost.org/mailman/listinfo.cgi/ublas
>>>
>>>
>>_______________________________________________
>>ublas mailing list
>>ublas_at_[hidden]
>>http://lists.boost.org/mailman/listinfo.cgi/ublas
>>
>>
>>
>_______________________________________________
>ublas mailing list
>ublas_at_[hidden]
>http://lists.boost.org/mailman/listinfo.cgi/ublas
>
>
>

-- 
==============================================
Look at our unique training program and
Register on-line at http://www.fft.be/?id=35
----------------------------------------------
Karl Meerbergen
Free Field Technologies
NEW ADDRESS FROM FEBRUARY 1st ONWARD:
Axis Park Louvain-la-Neuve
rue Emile Francqui, 1
B-1435 Mont-Saint Guibert - BELGIUM
Company Phone:  +32 10 45 12 26
Company Fax:    +32 10 45 46 26
Mobile Phone:   +32 474 26 66 59
Home Phone:     +32 2 306 38 10
mailto:Karl.Meerbergen_at_[hidden]
http://www.fft.be
============================================