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

Sandbo

Mon, 06/04/2018 - 10:02 am

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?

June 4, 2018 at 10:21 am - Permalink

Thanks Igor for the prompt reply,

I tried the following code, where the old way was commented out.

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:

•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

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

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

Have you tried the Multithread keyword as in

`Multithread momentWave=(ScWave^n)*(SWave^m)`

?June 4, 2018 at 11:50 pm - Permalink

(1)

(2)

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

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

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

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