bugGNU Octave - Bugs: bug #52706, Rewrite xpow.cc using templates


bug #52706: Rewrite xpow.cc using templates

Submitted by:  Marco Caliari <caliari>
Submitted on:  Wed 20 Dec 2017 11:53:09 AM UTC  
Category:  Octave Function Severity:  2 - Minor
Priority:  3 - Low Item Group:  Feature Request
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  Open
Release:  dev Operating System:  GNU/Linux

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( Jump to the original submission)

Fri 22 Dec 2017 11:09:10 PM UTC, comment #8: 

I assume you mean things like this that appear in Array.h?

    octave_idx_type i;
    for (i = 0; i < len - 3; i += 4)
        octave_quit ();

        p[i] = fcn (m[i]);
        p[i+1] = fcn (m[i+1]);
        p[i+2] = fcn (m[i+2]);
        p[i+3] = fcn (m[i+3]);

    octave_quit ();

    for (; i < len; i++)
      p[i] = fcn (m[i]);

Instead of duplicating code like this many times, maybe we could do this loop unrolling job with a template as well?

Also, I would not worry too much about calling octave_quit every time through a loop if the loop does any substantial amount of work since it is a function that should be easily inlined to be a check on a global variable.  Nothing more should happen unless an signal has been caught.  Maybe that affects performance enough to notice, but if it does, then I suspect we should be unrolling the loop more than 4 times.

We could also probably gain something by properly using function parameters in templates instead of function pointers.  This may not apply to xpow.cc, but I think it does in other places like bsxfun.cc.

John W. Eaton <jwe>
Project Administrator
Fri 22 Dec 2017 10:50:18 PM UTC, comment #7: 

xpow.cc is a mess.  It seems to have been written before widespread use of templates and there is a lot of code duplication.

I applied your patch here (http://hg.savannah.gnu.org/hgweb/octave/rev/4827cbef0949).  In addition, I wrote a templated version of xisint and used Octave coding syntax for indexing.

The whole file should really be re-written with templates, and to make use of xelem () for faster access to array elements.  In addition, using octave_quit() in the innermost loop is overkill.  In liboctave there are examples where the inner loop runs 4 times before calling octave_quit().  That could be done here as well.

Rik <rik5>
Project Administrator
Thu 21 Dec 2017 07:12:00 AM UTC, comment #6: 

I attached a diff.

(file #42693)

Marco Caliari <caliari>
Project Member
Wed 20 Dec 2017 07:24:12 PM UTC, comment #5: 

@Marco: the static_cast is correct.  It just doesn't need to happen twice, nor be part of the conditional.


  if (static_cast<int> (b) == b)
      int btmp = static_cast<int> (b);


  if (xisint (b))
      int btmp = static_cast<int> (b);

Rik <rik5>
Project Administrator
Wed 20 Dec 2017 07:13:45 PM UTC, comment #4: 

@Rik, @jwe: but how do I cast to integer?

int btmp = static_cast<int> (b); => ???

Marco Caliari <caliari>
Project Member
Wed 20 Dec 2017 06:14:07 PM UTC, comment #3: 

@Marco: jwe supplied the answer about xisint: cruft happens.

Rik <rik5>
Project Administrator
Wed 20 Dec 2017 06:09:46 PM UTC, comment #2: 

Rik: I would bet that the static_cast to int is just an artifact from ancient times.

John W. Eaton <jwe>
Project Administrator
Wed 20 Dec 2017 05:45:14 PM UTC, comment #1: 

Interesting.  I suppose there is no objection for the time being.  Eventually Octave might implement in-place operators for Octave values in which case

result = result * atmp;
result *= atmp;  // with efficiency gain

This change might not be possible with the Matlab-equivalent expression

result = atmp * result;

Still, I'm not sure in-place operators will ever happen so we shouldn't delay things based on an uncertain future.  Can you leave a comment in the code about why the ordering is important for Matlab compatibility?

I assume this only applies to matrices and not to scalars or the elem_xpow functions?

Also, while you're in xpow.cc, can you figure out why this code exists

  if (static_cast<int> (b) == b)
      int btmp = static_cast<int> (b);

It seems like it would be simpler to use xisint (b) which is also used throughout the code and is clearer about what is going on

  if (a < 0.0 && ! xisint (b))

Rik <rik5>
Project Administrator
Wed 20 Dec 2017 11:53:09 AM UTC, original submission:  

It is easy to check that, e.g., A^7 is computed in octave as

A2 = A * A;
A4 = A2 * A2;
(A * A2) * A4

while in matlab as

A2 = A * A;
A4 = A2 * A2;
A4 * (A2 * A)

In exact arithmetic it is obviously the same. The modification would be

result = result * atmp; -> result = atmp * result;

in xpow.cc. Any objections?

Marco Caliari <caliari>
Project Member


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

Attach Files:

Attached Files
file #42693:  xpow.diff added by caliari (5KiB - 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 jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-12-22 rik5 Priority5 - Normal => 3 - Low
        Item GroupMatlab Compatibility => Feature Request
        StatusPatch Submitted => Confirmed
        Summarympower compatibility behavior => Rewrite xpow.cc using templates
    2017-12-21 caliari Attached File- => Added xpow.diff, #42693
        StatusIn Progress => Patch Submitted
    2017-12-20 rik5 StatusNone => In Progress

    Back to the top

    Powered by Savane 3.5