| Matrix operations with IDL: Avoiding for loops [message #92487] |
Tue, 29 December 2015 16:57  |
vince33600
Messages: 3
Registered: December 2015
|
Junior Member |
|
|
Dear all,
I was trying to improve the performance of some pieces of code that are taking forever to run.
Basically, I'm trying to multiply a set of n matrix (3x3) by a set of n vectors (3x1) without using any for loops. The results of these operations should give me a set of n vectors (3x1).
Let's take a simplified example where n=2. Therefore, I have 2 matrixes (let's call them a and b) that needs to be multiplied to 2 vector (let's call them u and v).
I figured out that the operation could be done by reshaping (using rebin and reform for instance) the matrixes into a bigger array (let's call it M) whose diagonal elements are the a and b matrixes, so that:
M = | a 0 |
| 0 b |
where a and b are the 3x3 matrixes, and by reshaping the n vectors into in single vector (called I), so that:
I = | u |
| v |
Then, the results would be:
R = M.I
Finally, the n vectors would be obtained by reshaping the R vector into n (3x1) vector.
Coming for fortran, I initially coded that by decomposing every single matrix multiplication in a for loop. I then tried to apply the above solution, but it seems a real stretch for me to do it without any loops.
I was thinking that someone already might have faced that problem.
Thanks for your help!
Vincent
|
|
|
|
| Re: Matrix operations with IDL: Avoiding for loops [message #92495 is a reply to message #92487] |
Wed, 30 December 2015 20:47  |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
On Tuesday, December 29, 2015 at 7:57:30 PM UTC-5, vince...@gmail.com wrote:
> Dear all,
>
> I was trying to improve the performance of some pieces of code that are taking forever to run.
> Basically, I'm trying to multiply a set of n matrix (3x3) by a set of n vectors (3x1) without using any for loops. The results of these operations should give me a set of n vectors (3x1).
>
> Let's take a simplified example where n=2. Therefore, I have 2 matrixes (let's call them a and b) that needs to be multiplied to 2 vector (let's call them u and v).
>
> I figured out that the operation could be done by reshaping (using rebin and reform for instance) the matrixes into a bigger array (let's call it M) whose diagonal elements are the a and b matrixes, so that:
>
> M = | a 0 |
> | 0 b |
>
> where a and b are the 3x3 matrixes, and by reshaping the n vectors into in single vector (called I), so that:
>
> I = | u |
> | v |
>
> Then, the results would be:
> R = M.I
>
> Finally, the n vectors would be obtained by reshaping the R vector into n (3x1) vector.
>
> Coming for fortran, I initially coded that by decomposing every single matrix multiplication in a for loop. I then tried to apply the above solution, but it seems a real stretch for me to do it without any loops.
>
> I was thinking that someone already might have faced that problem.
For IDL, FOR loops are not a problem as long as you do a lot of work per iteration. Here is an example, where I literally do the matrix multiplication "by hand."
;; Set up some dummy inputs
m = randomn(seed,3,3,1000) ;; M = Your 3x3xN matrices
u = randomn(seed,3,1000) ;; U = Your 3xN vectors
v = u*0 ;; V = The final result
;; Boom! Write out one row of matrix multiplication and do
;; that operation thrice.
for i = 0, 2 do v(i,*) = m(0,i,*)*u(0,*) + m(1,i,*)*u(1,*) + m(2,i,*)*u(2,*)
No FOR loops but it's so fast, who cares. Even with 100x as many matrices on my six year old laptop, it takes barely any time at all.
Craig
|
|
|
|
| Re: Matrix operations with IDL: Avoiding for loops [message #92532 is a reply to message #92495] |
Tue, 05 January 2016 13:45 |
vince33600
Messages: 3
Registered: December 2015
|
Junior Member |
|
|
Le mercredi 30 décembre 2015 22:47:12 UTC-6, Craig Markwardt a écrit :
> On Tuesday, December 29, 2015 at 7:57:30 PM UTC-5, vince...@gmail.com wrote:
>> Dear all,
>>
>> I was trying to improve the performance of some pieces of code that are taking forever to run.
>> Basically, I'm trying to multiply a set of n matrix (3x3) by a set of n vectors (3x1) without using any for loops. The results of these operations should give me a set of n vectors (3x1).
>>
>> Let's take a simplified example where n=2. Therefore, I have 2 matrixes (let's call them a and b) that needs to be multiplied to 2 vector (let's call them u and v).
>>
>> I figured out that the operation could be done by reshaping (using rebin and reform for instance) the matrixes into a bigger array (let's call it M) whose diagonal elements are the a and b matrixes, so that:
>>
>> M = | a 0 |
>> | 0 b |
>>
>> where a and b are the 3x3 matrixes, and by reshaping the n vectors into in single vector (called I), so that:
>>
>> I = | u |
>> | v |
>>
>> Then, the results would be:
>> R = M.I
>>
>> Finally, the n vectors would be obtained by reshaping the R vector into n (3x1) vector.
>>
>> Coming for fortran, I initially coded that by decomposing every single matrix multiplication in a for loop. I then tried to apply the above solution, but it seems a real stretch for me to do it without any loops.
>>
>> I was thinking that someone already might have faced that problem.
>
> For IDL, FOR loops are not a problem as long as you do a lot of work per iteration. Here is an example, where I literally do the matrix multiplication "by hand."
>
> ;; Set up some dummy inputs
> m = randomn(seed,3,3,1000) ;; M = Your 3x3xN matrices
> u = randomn(seed,3,1000) ;; U = Your 3xN vectors
> v = u*0 ;; V = The final result
>
> ;; Boom! Write out one row of matrix multiplication and do
> ;; that operation thrice.
> for i = 0, 2 do v(i,*) = m(0,i,*)*u(0,*) + m(1,i,*)*u(1,*) + m(2,i,*)*u(2,*)
>
> No FOR loops but it's so fast, who cares. Even with 100x as many matrices on my six year old laptop, it takes barely any time at all.
>
> Craig
Thanks Craig!
I guess not every for loops are evil in IDL.
Thanks for the answer anyway,
Vincent
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
