Fast way to raise a complex array to power n (element wise)

Share icon

Hi,

I am wondering if there is a way to speed up the operation:

SWaveN = Swave^n


I tried to look into FastOp and matrixOp, but it seems they do not have the ^ operation available (sorry if I missed).
Would there be another way?
Igor wrote:
Check out MatrixOP under powR() or PowC().

Thanks Igor for the prompt reply,

I tried the following code, where the old way was commented out.
    Variable t1=startmstimer
    for (n=0;n<size;n=n+1)
        for (m=0;m<size;m=m+1)
            //momentWave=(ScWave^n)*(SWave^m)
           
            MatrixOp/o ScWaveN = powC(ScWave,n)
            MatrixOp/o SWaveM = powC(SWave,m)
            MatrixOp/o momentWave = ScWaveN*SWaveM
           
            MatrixOp/o theMean=mean(momentWave)
            momentMatOneTrig[n][m]=theMean[0]
            killwaves/z theMean
            endfor
        endfor
    print stopMsTimer(t1)


They turned out slower, plus the CPU utilization stays around 15% on my i7-4790:
//Without using MatrixOp, "momentWave=(ScWave^n)*(SWave^m)"
•fullMomentRecovery(1, "IQDataDG1_2058_0004", 11, 2.7637e6,1e6, 4.2044e9, -13.9,skip_u=1,startWave_u=4)
  5.79424e+07

//Using MatrixOp as shown above
•fullMomentRecovery(1, "IQDataDG1_2058_0004", 11, 2.7637e6,1e6, 4.2044e9, -13.9,skip_u=1,startWave_u=4)
  1.34118e+08

Did I apply it wrongly, or is there something I am missing?

June 4, 2018 at 11:06 am - Permalink

First, to be clear, the two MatrixOP functions work on an element by element basis. This is not a matrix to a power but the equivalent to the wave expression that you provided.

Not all MatrixOP functions have automatic multithreading so you need to be tricky in using it. How many points are in your wave?

June 4, 2018 at 12:49 pm - Permalink

Igor wrote:
First, to be clear, the two MatrixOP functions work on an element by element basis. This is not a matrix to a power but the equivalent to the wave expression that you provided.

Not all MatrixOP functions have automatic multithreading so you need to be tricky in using it. How many points are in your wave?

Yes, I understand and I intended to do element by element multiplication of two arrays. Thanks for the note and I just updated the title.

The waves are 5e6 points each.
Also, I am on Igor 7 Version: 7.0.9.1 (Build 31885) (Nightly build downloaded earlier today).

June 4, 2018 at 01:33 pm - Permalink

Sandbo wrote:
The waves are 5e6 points each.


Have you tried the Multithread keyword as in Multithread momentWave=(ScWave^n)*(SWave^m)?

June 4, 2018 at 11:50 pm - Permalink

At large number of n and m, equation(2) will be faster than equation(1).

(1)
momentWave = (ScWave^n)*(SWave^m)

(2)
SWave = r2polar(SWave)
sRadii = real(SWave)
sRadii = ln(sRadii)
sTheta = imag(SWave)

ScWave = r2polar(ScWave)
scRadii = real(ScWave)
scRadii = ln(scRadii)
scTheta = imag(ScWave)

//for(n=0; n<size; n=n+1)
scRadiiN = (n) * scRadii
scRadiiN = exp(scRadiiN)
scThetaN = (n) * scTheta

//for (m=0; m<size; m=m+1)
sRadiiM = (m) * sRadii
sRadiiM = exp(sRadiiM)
sThetaM = (m) * sTheta
momentWaveR = scRadiiN * sRadiiM
momentWaveT = scThetaN + sThetaM
momentWave = cmplx(momentWaveR, momentWaveT)
momentWave = p2rect(momentWave)
//endfor

//endfor

SWave = p2rect(SWave)
ScWave = p2rect(ScWave)

June 5, 2018 at 05:41 am - Permalink

Sandbo wrote:
Igor wrote:
Check out MatrixOP under powR() or PowC().

Thanks Igor for the prompt reply,

I tried the following code, where the old way was commented out.
    Variable t1=startmstimer
    for (n=0;n<size;n=n+1)
        for (m=0;m<size;m=m+1)
            //momentWave=(ScWave^n)*(SWave^m)
           
            MatrixOp/o ScWaveN = powC(ScWave,n)
            MatrixOp/o SWaveM = powC(SWave,m)
            MatrixOp/o momentWave = ScWaveN*SWaveM
           
            MatrixOp/o theMean=mean(momentWave)
            momentMatOneTrig[n][m]=theMean[0]
            killwaves/z theMean
            endfor
        endfor
    print stopMsTimer(t1)


They turned out slower, plus the CPU utilization stays around 15% on my i7-4790:
//Without using MatrixOp, "momentWave=(ScWave^n)*(SWave^m)"
•fullMomentRecovery(1, "IQDataDG1_2058_0004", 11, 2.7637e6,1e6, 4.2044e9, -13.9,skip_u=1,startWave_u=4)
  5.79424e+07

//Using MatrixOp as shown above
•fullMomentRecovery(1, "IQDataDG1_2058_0004", 11, 2.7637e6,1e6, 4.2044e9, -13.9,skip_u=1,startWave_u=4)
  1.34118e+08

Did I apply it wrongly, or is there something I am missing?


It looks to me like you only have two source waves ScWave and SWave that you are raising to the power of 1, 2, 3, 4, 5 etc... If that is true couldn't you use the previous result for ScWave^n to calculate the next value ScWave^(n+1) simply as the old value times ScWave? Wouldn't that be faster, especially for large values of n?

It also looks like you are recalculating ScWave^n for each value of SWave^m. I would suggest to first calculate all values of ScWave^n then all values of SWave^m and then do the multiplication of the two.

June 5, 2018 at 05:35 am - Permalink

olelytken wrote:
Sandbo wrote:
Igor wrote:
Check out MatrixOP under powR() or PowC().

Thanks Igor for the prompt reply,

I tried the following code, where the old way was commented out.
    Variable t1=startmstimer
    for (n=0;n<size;n=n+1)
        for (m=0;m<size;m=m+1)
            //momentWave=(ScWave^n)*(SWave^m)
           
            MatrixOp/o ScWaveN = powC(ScWave,n)
            MatrixOp/o SWaveM = powC(SWave,m)
            MatrixOp/o momentWave = ScWaveN*SWaveM
           
            MatrixOp/o theMean=mean(momentWave)
            momentMatOneTrig[n][m]=theMean[0]
            killwaves/z theMean
            endfor
        endfor
    print stopMsTimer(t1)


They turned out slower, plus the CPU utilization stays around 15% on my i7-4790:
//Without using MatrixOp, "momentWave=(ScWave^n)*(SWave^m)"
•fullMomentRecovery(1, "IQDataDG1_2058_0004", 11, 2.7637e6,1e6, 4.2044e9, -13.9,skip_u=1,startWave_u=4)
  5.79424e+07

//Using MatrixOp as shown above
•fullMomentRecovery(1, "IQDataDG1_2058_0004", 11, 2.7637e6,1e6, 4.2044e9, -13.9,skip_u=1,startWave_u=4)
  1.34118e+08

Did I apply it wrongly, or is there something I am missing?


It looks to me like you only have two source waves ScWave and SWave that you are raising to the power of 1, 2, 3, 4, 5 etc... If that is true couldn't you use the previous result for ScWave^n to calculate the next value ScWave^(n+1) simply as the old value times ScWave? Wouldn't that be faster, especially for large values of n?

It also looks like you are recalculating ScWave^n for each value of SWave^m. I would suggest to first calculate all values of ScWave^n then all values of SWave^m and then do the multiplication of the two.


That's a good catch, I totally missed that and spent a lot of extra computation time. Thanks for pointing it out.
I should probably create a higher dimension wave and compute the ScWave^n and SWave^m recursively first, then call them in the ScWaveN*SWaveM loop.

June 5, 2018 at 08:29 am - Permalink

thomas_braun wrote:
Sandbo wrote:
The waves are 5e6 points each.


Have you tried the Multithread keyword as in Multithread momentWave=(ScWave^n)*(SWave^m)?

Thanks for the info, I didn't know I should do this. Now the CPU utilization goes up to 50% most of the time and they run much faster.

June 5, 2018 at 08:33 am - Permalink