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!
