unravel calculation

[dkokron]

I'm hoping some blitz++ expert can help me unravel the following calculation.

int ub1=46+1, ub2=18+1, ub3=5+1, ub4=6+1;
Array<double, 4> rjac(ub1,ub2,ub3,ub4);
Array<double, 2> fjac(ub2,ub4);
Array<double, 3> pv(ub1,ub2,ub3);
rjac += fjac(i2,i4) * pv(i1,i2,i3);

My attempt to convert this into a simple loop nest gives the wrong results. Can anyone tell me what I'm doing wrong?

const double * restrict pvp = pv.data();
const double * restrict fjacp = fjac.data();
double * restrict rjacp = rjac.data();
for (int o=0; o<oub; ++o) { // 0:46
for (int n=0; n<nub; ++n) { // 0:18
for (int m=0; m<mub; ++m) { // 0:5
for (int l=0; l<lub; ++l) { // 0:6
rjacp[o*(lubmubnub) + n*(lubmub) + mlub + l] += fjacp[nlub + l] * pvp[o(nubmub) + nmub + m];
}
}
}
}

[slayoo]

Please provide an example that compiles. In the code you provided, i1, i2, i3, i4 are undefined.
You can use backticks (https://docs.github.com/en/free-pro-team@latest/github/writing-on-github/creating-and-highlighting-code-blocks) to paste C++ code conveniently, so Github will preserve indentation and do syntax highlighting. E.g.,.

#include <blitz/array.h>

int main()
{
    using namespace blitz;
    ...
}

[dkokron]

Any hints as to what I might be doing wrong?

[dkokron]

ForBlitzDev.tar.gz
This version reads inputs from the real application. The #if switch controls whether the code goes through my loop nest or the blitz code.

The trouble is that both paths give the same wrong answer (rjac.out holds the correct result from the real application).

cmp rjac.out rjac.out2
rjac.out rjac.out2 differ: byte 345, line 1

cmp rjac.out rjac_dkokron.out
rjac.out rjac_dkokron.out differ: byte 345, line 1

Many thanks for taking a look at this.
Dan

[slayoo]

Dan, IIUC, both the blitz-based and hand-coded-loops versions give the same result, but it is not what you expect? right?

[dkokron]

Correct. I was hoping there was some obvious error in the reproducer that
you could spot. If not, then I'll need to dig deeper into the real
application to figure out what is going on.

[dkokron]

forgot to include the blitz patch in the tarball
blitz++.patch.gz

[slayoo]

Note that what is in this patch has already been corrected in the repo long time ago (2016):
c07fec9
dbb5fef

There have been a few releases since 0.9, with latest being 1.0.2 from 2019: https://github.com/blitzpp/blitz/releases

[dkokron]

The real application gives the correct results using blitz-0.9 plus the patch. If I replace the blitz computation with my loop-nest in the real application, it gives the wrong results. The same wrong results as both blitz and the loop-nest give from the reproducer. That's why I'm confused.

[dkokron]

Correct. I was hoping there was some obvious error in the reproducer that you could spot. If not, then I'll need to dig deeper into the real application to figure out whats going on.

…

On Fri, Dec 11, 2020, 6:20 PM Sylwester Arabas ***@***.***> wrote: Dan, IIUC, both the blitz-based and hand-coded-loops versions give the same result, but it is not what you expect? right? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#167 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACODV2CRASCLFEYXIBWGPG3SUKZMDANCNFSM4UVF27JQ> .

[slayoo]

Perhaps there is some Array.reindex() or .reindexSelf() call in between in the "real application"?
See https://github.com/blitzpp/blitz/blob/master/testsuite/reindex.cpp