Re: [ublas] Solving linear equations of the form Ax = b ?

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:  [phone]
Company Fax:    [phone]
Mobile Phone:   [phone]
Home Phone:     [phone]
mailto:Karl.Meerbergen_at_[hidden]
http://www.fft.be
============================================