bugGNU Octave - Bugs: bug #53440, Wrong factorization with ilutp type

 
 

bug #53440: Wrong factorization with ilutp type

Submitter:  Marco Caliari <caliari>
Submitted:  Fri 23 Mar 2018 01:31:47 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Tue 27 Mar 2018 10:25:04 PM UTC, comment #4: 

@Marco: There is an easier way to compute the transposed, permuted identity matrix.  Instead of creating the speye matrix, permuting it via indexing, and then calculationg the transpose of the result, it is possible te just transpose the permuting index and then use that in a single indexing step.

Here is the code:


-            retval(0) = L + speye.index (perm, idx_vector::colon);
+            retval(0) = L + speye.index (idx_vector::colon, perm);


I benchmarked the original code construction and this method for speye(1e3) and it was only 18% faster, but that is still worth something.

I checked in the patch here (http://hg.savannah.gnu.org/hgweb/octave/rev/acb7094aebb3).

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Tue 27 Mar 2018 09:33:43 PM UTC, comment #3: 

As far I can see Marco is totally right and this patch should go in the 4.4.0 release.

The documentation error was accidentally introduced right at the beginning with https://hg.savannah.gnu.org/hgweb/octave/rev/76a6ba7d65d0 (don't want to blame anyone, I did not check carefully enough myself).

@Marco: At the moment I do not have the time to really look into this.  Your approach improves the situation until the function gets a makeover some day.

Kai Torben Ohlhus <siko1056>
Group Member
Tue 27 Mar 2018 05:23:14 PM UTC, comment #2: 

@Marco: Should this patch be applied for 4.4.0 release or is there a risk associated with checking it in now?

Rik <rik5>
Group administrator
Fri 23 Mar 2018 02:29:31 PM UTC, comment #1: 

Well, it was quite easy. Is there a better way to compute the transposed permuted identity?

(file #43649)

Marco Caliari <caliari>
Group Member
Fri 23 Mar 2018 01:31:47 PM UTC, original submission:  

According to the documentation (and matlab behavior) the following code


A=sprand(10,10,1/2);opts.type="ilutp";[L,U]=ilu(A,opts);norm(L*U-A,inf)


should give something like eps. In fact, droptol is zero, which means complete LU factorization. The explicit pivoted version


A=sprand(10,10,1/2);opts.type="ilutp";[L,U,P]=ilu(A,opts);norm(L*U-P*A,inf)


works fine.

Marco Caliari <caliari>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #43649:  ilu.diff added by caliari (3KiB - application/x-tex)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Submitted the item)
  • -email is unavailable- added by caliari
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-03-27 rik5 StatusPatch Reviewed Fixed
        Open/ClosedOpen Closed
    2018-03-27 siko1056 StatusPatch Submitted Patch Reviewed
    2018-03-23 caliari StatusNone Patch Submitted
    2018-03-23 caliari Attached File- Added ilu.diff, #43649
    2018-03-23 caliari Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code