Just curious--in this case how slow? Anyway, you should check out the FastAnonymous.jl package--it should help in your case.
Cheers, Kevin On Wednesday, February 11, 2015, Andrei Berceanu <[email protected]> wrote: > I have written the following function for generating a sparse matrix: > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > *function genspmat(ω0::Float64; N=pm["N"], α=pm["α"], γ=pm["γ"], > κ=pm["κ"]) # Determine memory usage nz = countnonzeros(; N=N) # > Preallocate I = Array(Int64,nz) J = Array(Int64,nz) V = > Array(Complex{Float64},nz) function setnzelem(i,n,m; pos="self") > if pos=="left" k += 1 J[k] = i-N; I[k] = i; V[k] = > 1 elseif pos=="right" k += 1 J[k] = i+N; I[k] > = i; V[k] = 1 elseif pos=="up" k += 1 J[k] = > i-1; I[k] = i; V[k] = exp(-im*2π*α*m) elseif pos=="down" > k += 1 J[k] = i+1; I[k] = i; V[k] = exp(im*2π*α*m) elseif > pos=="self" k += 1 J[k] = i; I[k] = i; V[k] = ω0 + > im*γ - 1/2*κ*(m^2+n^2) end end # maximum value of > m or n indices maxm = div(N-1,2) k = 0 for i in 1:N^2 m = > getm(i; N=N) n = getn(i; N=N) #self interaction is always > present setnzelem(i,n,m) #corners #top left if > n==maxm && m==-maxm setnzelem(i,n,m; pos="right") > setnzelem(i,n,m; pos="down") #top right elseif n==maxm && > m==maxm setnzelem(i,n,m; pos="left") setnzelem(i,n,m; > pos="down") #bottom right elseif n==-maxm && m==maxm > setnzelem(i,n,m; pos="left") setnzelem(i,n,m; > pos="up") #bottom left elseif n==-maxm && m==-maxm > setnzelem(i,n,m; pos="right") setnzelem(i,n,m; > pos="up") #edges #top elseif n == maxm > setnzelem(i,n,m; pos="right") setnzelem(i,n,m; > pos="left") setnzelem(i,n,m; pos="down") #right > elseif m == maxm setnzelem(i,n,m; pos="left") > setnzelem(i,n,m; pos="up") setnzelem(i,n,m; pos="down") > #bottom elseif n == -maxm setnzelem(i,n,m; > pos="left") setnzelem(i,n,m; pos="up") > setnzelem(i,n,m; pos="right") #left elseif m == > -maxm setnzelem(i,n,m; pos="down") setnzelem(i,n,m; > pos="up") setnzelem(i,n,m; pos="right") else > #bulk setnzelem(i,n,m; pos="down") setnzelem(i,n,m; > pos="up") setnzelem(i,n,m; pos="right") > setnzelem(i,n,m; pos="left") end end return sparse(I,J,V)end* > Notice that inside this function I have defined *setnzelem(i,n,m; > pos="self"), *which is where all the action really takes place :) > Now I would like to generalize the *genspmat* to values V[k] which can be > arbitrary functions of m and n. One way of doing so is to define > > > *function > genspmat(l::Function,r::Function,u::Function,d::Function,s::Function, > N::Int,nz::Int)*and then inside setnzelem do > > function setnzelem(i::Int,n::Int,m::Int; pos="self") > if pos=="left" > k += 1 > J[k] = i-N; I[k] = i; V[k] = l(m,n) > elseif pos=="right" > k += 1 > J[k] = i+N; I[k] = i; V[k] = r(m,n) > elseif pos=="up" > k += 1 > J[k] = i-1; I[k] = i; V[k] = u(m,n) > elseif pos=="down" > k += 1 > J[k] = i+1; I[k] = i; V[k] = d(m,n) > elseif pos=="self" > k += 1 > J[k] = i; I[k] = i; V[k] = s(m,n) > end > end > > The problem with this is that passing functions as arguments tends to be > slow in Julia. > So what is your advice on accomplishing what I want? > > Tnx! > > >
