It looks like the qr and the eigvecs are both done on 3x3 matrices. You may
be better off just writing a special 3x3 in-place implementation of these
functions and see if that is much faster. I am guessing quite a lot of time
is spent in calling LAPACK and temporary workspace allocation.
In general, it would be nice to be able to reuse the LAPACK workspace by
allocating it once and reusing it across multiple calls. We do not have
that yet.
This may be worth filing a performance issue.
-viral
On Wednesday, November 19, 2014 6:39:49 AM UTC+5:30, Ryan Young wrote:
>
> Hello,
>
> I am new to julia and am working on some neuroimaging analysis. I have a
> function that calculates the Jacobian of the deformation field at each
> point in the deformation field. A qr decomposition is then done on the
> Jacobian. The principle eigenvector of the strain portion of the QR
> decomposition is rotated, stored in an array and returned.
>
> I was not going to parallelize this function as I think it is easier to do
> a pmap of this process over all of the samples. Currently, this is taking a
> few minutes to run on one large deformation field (225,425,225,3). All the
> time is being taken up in the last three lines of the for loop. Any advice
> is greatly appreciated.
>
> Thanks
> Ryan
>
>
> function calcJac(deformation::Array{Float32,4})
> temp=Array(Float32,3,3)
> r=Array(Float32,3,3)
> s=Array(Float32,3,3)
> ev=Array(Float32,3,3)
>
> i,j,k=size(deformation)
> final=Array(Float32,i,j,k,3)
>
> for x in 3:i-3
> for y in 3:j-3
> for z in 3:k-3
>
> #compute frist row of jacobain
>
> temp[1,1]=(-deformation[x+2,y,z,1]+8*deformation[x+1,y,z,1]-8*deformation[x-1,y,z,1]+deformation[x-2,y,z,1])/12.0
>
> temp[1,2]=(-deformation[x+2,y,z,2]+8*deformation[x+1,y,z,2]-8*deformation[x-1,y,z,2]+deformation[x-2,y,z,2])/12.0
>
> temp[1,3]=(-deformation[x+2,y,z,3]+8*deformation[x+1,y,z,3]-8*deformation[x-1,y,z,3]+deformation[x-2,y,z,3])/12.0
>
> #compute 2nd row of jacobian
>
> temp[1,1]=(-deformation[x,y+2,z,1]+8*deformation[x,y+1,z,1]-8*deformation[x,y-1,z,1]+deformation[x,y-2,z,1])/12.0
>
> temp[1,2]=(-deformation[x,y+2,z,2]+8*deformation[x,y+1,z,2]-8*deformation[x,y-1,z,2]+deformation[x,y-2,z,2])/12.0
>
> temp[1,3]=(-deformation[x,y+2,z,3]+8*deformation[x,y+1,z,3]-8*deformation[x,y-1,z,3]+deformation[x,y-2,z,3])/12.0
>
> #compute third row of jacobian
>
> temp[1,1]=(-deformation[x,y,z+2,1]+8*deformation[x,y,z+1,1]-8*deformation[x,y,z-1,1]+deformation[x,y,z-2,1])/12.0
>
> temp[1,2]=(-deformation[x,y,z+2,2]+8*deformation[x,y,z+1,2]-8*deformation[x,y,z-1,2]+deformation[x,y,z-2,2])/12.0
>
> temp[1,3]=(-deformation[x,y,z+2,3]+8*deformation[x,y,z+1,3]-8*deformation[x,y,z-1,3]+deformation[x,y,z-2,3])/12.0
>
> r,s=qr(temp)
> ev=eigvecs(s)
> final[x,y,z,:]=r*ev[:,1]
>
> end
> end
> end
>
> return final
> end
>
