bugGNU Octave - Bugs: bug #67741, Add "funm()" function to...

 
 

bug #67741: Add "funm()" function to Octave

Submitter:  Dmitri A. Sergatskov <dasergatskov>
Submitted:  Thu 27 Nov 2025 05:40:25 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  In Progress Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  11.1.0 (current default) Planned Release:  11.1.0 (current default)
* Mandatory Fields

Post a Comment

Add a New Comment Rich Markup
   

Discussion

Jump to the original submission

Fri 05 Dec 2025 04:31:28 PM UTC, comment #22: 

@Dmitri, comment #18:
OK so I'll do the "more cosmetic things" first and once that is up the "innards" can be treated, OK?

As to comment #19, at the time it occurred to me that maybe we should have first added the funm.m from Mudit Sharma's repo as-is immediately followed by the fixes needed to get it working more reliably.  That is also what I usually did/do with contributions for the packages I maintain.
But in hindsight the Dec, 1th deadline + all the AI bot-induced repo problems interfered quite a bit with calm thinking.

If so required the vector optimization changes of the innermost loop can be put up here in a comment.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 05 Dec 2025 02:31:12 PM UTC, comment #21: 

207: if (isequal (A, triu (A)))

I think the code here is meant "if (isdiag(A))"
I looked at it shortly; it looks to me that the check is overly strict and should be allowed to have some noise on the off-diagonal elements...

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 05 Dec 2025 01:55:26 PM UTC, comment #20: 

I added functions trexc and funm to the NEWS file and to Octave's manual, as well as updating the seealso links for schur and ordschur to have an entry for trexc.  See https://hg.savann ... /rev/6b17f22b9526.


Rik <rik5>
Group administrator
Fri 05 Dec 2025 12:55:29 PM UTC, comment #19: 

@Rik,

The code definitely could be improved, but I wanted just to do minimal things to get things working so we have a base version
with as much "original" code as possible, so the following changes can be tracked in HG etc... At this moment I even kind-of regret that I did this vector optimization of the innermost loop since this change is not tracked anywhere.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 05 Dec 2025 12:48:13 PM UTC, comment #18: 

Philip,

I could help, but at this moment you should lead the effort.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 05 Dec 2025 11:57:33 AM UTC, comment #17: 

Okay, I'll leave it to you and Dmitri.  Thanks.

Rik <rik5>
Group administrator
Fri 05 Dec 2025 11:18:36 AM UTC, comment #16: 

First thing to be done is push the cset for NEWS entries (see comment #9), AFAICS that still has to happen.

As to "where to continue": most of all it is better to keep it all in one location, both on Discourse and Savannah isn't very practical.
Personally I'd prefer following up here, but lets hear what Dmitri wants (he who did most of the work).
If we keep going here, should Status be back to "In progress"?

Thanks for the list of actions.  Unless you can't wait I'd say let Dmitri and me pick up on further improvements.

As I see it, funm.m as it stood was a mere collection of separate functions written by Higham and its team.  Mudit Sharma had the first throw at it to glue it all into a working funm.m, but we're unsure if he got it working completely (probably not).  Dmitri finally tied up the loose ends.
So little wonder that .g., input validation is left in subfuncs.

@Dmitri - if you agree I'll try to do as much as I can w.r.t. style and input validation + input validation BISTs. I'd rather let LAPACK stuff for you, that's not quite my proficiency.
Most of all, we should avoid to do double work independently of each other.

I'm busy between the scenes with adding newer Matlab-compatible input, i.e., an options struct rather than position dependent params etc.  I wrote "adding" because older scripts will continue to work then and it is easy to implement the newer input options. Most of that is done, I got distracted by private things, but this weekend I hope to finish it.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 05 Dec 2025 10:48:59 AM UTC, comment #15: 

Do we want to continue further development here?  Or on Discourse thread (https://octave.di ... n-of-funm/7040/48)?

Now that it is in Mercurial we have a common point to work off of.  I just took a brief glance and there are things to improve.  On Octave coding conventions:
1) Standard copyright block at top of file
2) parentheses around argument to switch statement
3) No use of "otherwise" if unused
4) Prefer rows(), columns(), or numel() to length () for clarity
5) Demo %! blocks immediately follow code (before BISTS)
6) Input validation BISTs at end
7) Use '##' for full-line comments
8) Generally, all function code should end with ';'.  I noticed lines with disp() that did not.  Also, *return statement should be followed by semicolon ';'

Actual code improvements
1) Use istriu/istril
Here is one example, but there are others
```
207: if (isequal (A, triu (A)))
->
207: if (istriu (A))
```

2) Is it better to use LAPACK routine for solving triangular Sylvester equation?
The m-file defines function sylv_tri, but Octave already has a general purpose function sylvester which solves Sylvester equation.

3) Eliminate internal input validation on subfunctions.  Input validation should be done for the top-most function, and then it should arrange never to call subfunctions in an incorrect manner.  As one example, consider funm_atom which has this code

function [F, n_terms] = funm_atom (T, fun, tol, prnt)

  if (nargin < 3 || isempty (tol))
    tol = eps;
  endif
  if (nargin < 4)
    prnt = 0;
  endif

Just arrange to always call funm_atom with 4 arguments.

I can make some of the changes, but I don't know the algorithm well enough to have an opinion on, for example, replacing sylv_tri.


Rik <rik5>
Group administrator
Mon 01 Dec 2025 03:43:48 PM UTC, comment #14: 

I think we can continue to improve on the funm code further.
It definitely dows not look "finished", there are part of code in BLOCKING that does not appear to be used at all.

But it all looks like a significant refactoring / re-writing.
Not a quick fix.

But at the moment we have a working code that can be used as a "base" to which one compare new and improved code.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 01 Dec 2025 03:33:19 PM UTC, comment #13: 

Thanks for the explanation about the LAPACK stuff.  So trexc is much more useful (for those who know) than it looks at first sight :-)

Just beating on the dead horse:
AFAICS ordschur was called with a 2-element vector rather than a 3-element vector that it expected / needed; IMO that was a fault of funm as it stood. I found the part of the code responsible for creating that input vector (funm's internal 'swapping' function) to be a snakehole, and maybe it even got wrong info in the first place from another previously called internal function 'blocking' that perhaps didn't do its job properly. I'm not that bad in my linalg but I completely lost track of what was happening.
Maybe the funm.m code we got from Mudit Sharma's github repo was simply in a transitional if not unfinished state and you connected the missing bits.

All in all, good that you circumvented and fixed the lot.

Philip Nienhuis <philipnienhuis>
Group Member
Mon 01 Dec 2025 02:40:01 PM UTC, comment #12: 

No, it is not that simple. Lapack has 3 different Schur ordering function that do different things. `trsen` (that is what `ordschur` is using), `trexc`, and `laexc`.

`trsen` is the highest level that uses `trexc` which in turn  uses `laexc` (the lowest level).

For the current code in funm, the `trexc` is the exactly what it supposed to be used for. But, technically speaking `trsen` can do its own blocking/swapping, so yes, `funm` can be-re-written to use `trsen` (and hence, `ordeschur` once it has a complete interface), but it is a significant code re-write (essentially replacing BLOCKING and SWAPPING and may be something else) and perhaps that was the initial idea. I could not figure it all out in the time allotted.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 01 Dec 2025 02:27:49 PM UTC, comment #11: 

Only (tiny) gripe I might have is that AFAICS trexc essentially does what ordschur was supposed to do in the first place, or do I get that wrong? (I still believe that not ordschur was the culprit but rather that funm was calling it erroneously.)

(I know very well, the code freeze + continuous AI bot attacks left all of us with little time for contemplation, so you acted flexibly and did what could be done ~ "if it can't be done the way it should, it should be done the way it can".)

[OT]
BTW there's a patch for ordschur's missing input argument options, see bug #67740, I just built Octave with it and will test & review a bit the next days. If you would have a chance to do that as well, please do; you know the background of those functions better than me. I'll put you in the cc there, hope you don't mind.

Philip Nienhuis <philipnienhuis>
Group Member
Mon 01 Dec 2025 12:44:09 PM UTC, comment #10: 

It looks fine.

Thanks!

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 01 Dec 2025 12:35:39 PM UTC, comment #9: 

Good idea.
New cset attached.

@Dmitri, please check your name etc (+ apologies for misspelling it in previous csets). I copied it from Savannah users info

(file #57884)

Philip Nienhuis <philipnienhuis>
Group Member
Mon 01 Dec 2025 12:17:34 PM UTC, comment #8: 

Probably you should add `trexc` as well.
Though it is written to be used primarily in `funm` it is a user-callable function that can be used in other (most likely Schur-related) code.

Dmitri.
--


Dmitri A. Sergatskov <dasergatskov>
Mon 01 Dec 2025 12:13:08 PM UTC, comment #7: 

NEWS entry cset (overlooked in previous patch)

(file #57883)

Philip Nienhuis <philipnienhuis>
Group Member
Sun 30 Nov 2025 07:37:31 PM UTC, comment #6: 
Philip Nienhuis <philipnienhuis>
Group Member
Sun 30 Nov 2025 03:58:39 PM UTC, comment #5: 

Uploaded a cset with proper bug number ...

(file #57880)

Philip Nienhuis <philipnienhuis>
Group Member
Sun 30 Nov 2025 03:55:24 PM UTC, comment #4: 

Added cset based on funm_add_02.diff

(file #57879)

Philip Nienhuis <philipnienhuis>
Group Member
Sat 29 Nov 2025 05:09:27 PM UTC, comment #3: 

Cleaned up files a little more.
trexc.cc: Fix !%tests: Seeded randn(), adjusted tolerance

funm.m: Added extra input validation, convert input matrix to double explicitly (so that all internal calculations are done with double precision), but return single for single input, return double for int input. Added more tests.
Removed stray TABs.

file (funm_add_02.diff)

Dmitri.
--


(file #57877)

Dmitri A. Sergatskov <dasergatskov>
Thu 27 Nov 2025 08:14:37 PM UTC, comment #2: 

Where do we follow up, here or in patch #10067?

Provisionally here then:

I just built Octave 11 on Fedora 42 with the funm patch and got 1 FAIL for trexc:

>>>>> processing /home/philip/devel/octdev/oct1100+/libinterp/corefcn/trexc.cc-tst
***** test
 A = randn (4) + 1i * randn (4);
 [U, T] = schur (A, "complex");
 M = [1, 3; 2, 4];
 [U2, T2] = trexc (U, T, M);
 assert (U2 * T2 * U2', A, 20*eps);
!!!!! test failed
ASSERT errors for:  assert (U2 * T2 * U2',A,20 * eps)

  Location  |  Observed  |  Expected  |  Reason
   (3,3)     2.7737+1.5639i 2.7737+1.5639i   Abs err 4.6204e-15 exceeds tol 4.4409e-15 by 2e-16

so, just a tolerance thing.

Later tonight I'll make a crossbuild for Windows

Philip Nienhuis <philipnienhuis>
Group Member
Thu 27 Nov 2025 07:28:03 PM UTC, comment #1: 

Actually a duplicate of patch #10067.

Philip Nienhuis <philipnienhuis>
Group Member
Thu 27 Nov 2025 05:40:25 PM UTC, original submission:  

Currently octave is missing "funm()" function.
The attached patch implements it.

See
https://octave.di ... n-of-funm/7040/43
for details.

The attached patch is based on an old script written by Mudit Sharma and Marco Caliari, which in-turn based on work by Awad H. Al-Mohy and Nicholas J. Higham.

See:

Davies, P. I. and N. J. Higham, “A Schur-Parlett algorithm for computing matrix functions,” SIAM J. Matrix Anal. Appl., Vol. 25, Number 2, pp. 464-485, 2003.

for theoretical background.

The original script by Mudit Sharma has been re-styled to fit Octave's coding standards. There have been a number of bug fixes and improvements made. In particular a completely wrong code path that was using "ordschur" is replaced by a code calling Lapack's "ZTREXC" function for Schur re-ordering. The wrapper for this function, that provides a function "trexc" is included in the patch.

(file funm_add.diff)

Code compiles and "make check" passes on my computer (MacOS/M4/Homebrew).

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>

 

Attached Files

Attached Files
file #57884:  news_funm_trexc.cset added by philipnienhuis (607B - application/octet-stream)
file #57883:  news_funm.cset added by philipnienhuis (579B - application/octet-stream)
file #57880:  bug67741_patch10067_funm.cset added by philipnienhuis (35KiB - application/octet-stream)
file #57879:  bug64770_patch10067_funm.cset added by philipnienhuis (35KiB - application/octet-stream)
file #57877:  funm_add_02.diff added by dasergatskov (35KiB - text/x-patch)
file #57873:  funm_add.diff added by dasergatskov (31KiB - application/octet-stream)

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

Attach Files:
   
   
Comment:
   

 

Dependencies

This item does not depend on any other items.

No items depend on this one.

 

Mail Notification Carbon-Copy List

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

    Votes

    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.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

    History

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2025-12-05 rik5 StatusFixed In Progress
        Fixed ReleaseNone 11.1.0 (current default)
        Planned ReleaseNone 11.1.0 (current default)
    2025-12-01 philipnienhuis Attached File- Added news_funm_trexc.cset, #57884
    2025-12-01 philipnienhuis Attached File- Added news_funm.cset, #57883
    2025-11-30 philipnienhuis StatusNone Fixed
    2025-11-30 philipnienhuis Attached File- Added bug67741_patch10067_funm.cset, #57880
    2025-11-30 philipnienhuis Attached File- Added bug64770_patch10067_funm.cset, #57879
    2025-11-29 dasergatskov Attached File- Added funm_add_02.diff, #57877
    2025-11-27 dasergatskov Attached File- Added funm_add.diff, #57873

    Back to the top

    Powered by Savane 3.16-a7ba.
    Corresponding source code