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


bug #52706: Rewrite xpow.cc using templates

Submitter:  Marco Caliari <caliari>
Submitted:  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
* Mandatory Fields

Add a New Comment (Rich Markup)

Comment Type & Canned Response:


( 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)

    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 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.9