How to output selective elements from Matrix?

I have the below configuration. I am trying to output selective elements from matrix L.sl based on values of L.Seg. More specifically, I want to write syntax such as L.sl[L.Seg] to give me the same results as S. Any simple way to do that?

using SparseArrays, LinearAlgebra

Base.@kwdef mutable struct Ls
   tol::Vector{Float64} = [0,0,0]
   Seg::Vector{Int8} = [1,2,3]
   sl = []
end


L = Ls[];
push!(L, Ls());
L[1].sl = [1 2 3; 4 5 6; 7 8 9];
push!(L, Ls());
L[2].sl = [1 2 3; 4 5 6; 7 8 9];

@eval L = Ls((reduce(vcat, getfield.(L, pn)) for pn in fieldnames(Ls))...);

S = zeros(length(L.tol));
for j in eachindex(L.tol)
   S[j] = L.sl[j,L.Seg[j]];
end
#########################################

julia> L.sl
6Ă—3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9
 1  2  3
 4  5  6
 7  8  9

julia> S
6-element Vector{Float64}:
 1.0
 5.0
 9.0
 1.0
 5.0
 9.0

julia> L.sl[L.Seg]
6-element Vector{Int64}:
 1
 4
 7
 1
 4
 7

Any comment here please?

I’m not sure I completely understand what you’re after. But it sounds like you might be looking for CartesianIndex? For example,

inds = CartesianIndex.(eachindex(L.Seg), L.Seg)
S = L.sl[inds]

You could also accomplish this with broadcasted getindex, although it’s a bit less typical and can be trickier to understand

S = getindex.(Ref(L.sl), eachindex(L.Seg), L.Seg)

Thank you very much. Exactly this is what I need.
I have another issue please. Below is an MWE, in which I am trying to update the sparse matrix S. However, non of the two methods are working. Could you please suggest a solution?

using SparseArrays, LinearAlgebra

function dump(R)
  R.nzval .= [7,8,9];
end
S = sparse([1 0 0 0 0; 0 0 0 0 0; 0 0 1 0 0; 0 0 0 2 0; 0 0 0 0 5])
# Method-1
dump(S[3:end,3:end])
5Ă—5 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  â‹…  â‹…  â‹…  â‹…
 â‹…  â‹…  â‹…  â‹…  â‹…
 â‹…  â‹…  1  â‹…  â‹…
 â‹…  â‹…  â‹…  2  â‹…
 â‹…  â‹…  â‹…  â‹…  5
# Method-2
dump( @view( S[3:end,3:end]) )
ERROR: type SubArray has no field nzval

You should think very hard before mucking with the internal fields of objects, like you’re doing with R.nzval .= [7,8,9]. It’s generally a bad idea that can break things or at least be difficult to reason about and maintain.

In the dump(S[3:end,3:end]) version, S[3:end,3:end] copies the indexed part of the array, so changing it’s nzval doesn’t do anything to the original S. You need to return R if you want this modified array to be accessible outside of dump (although note that the first two rows/columns of S were never a part of R and will not be returned).

The @views version puts a wrapper object around the original array. The wrapper does not have a nzvals field to access, hence your error.

Thank you very much. It is clear now.
For performance purposes, what is the best way to mutate S, as below?

function dump(R)
  R.nzval .= [7,8,9];
return R
end
S = sparse([1 0 0 0 0; 0 0 0 0 0; 0 0 1 0 0; 0 0 0 2 0; 0 0 0 0 5])
# Method-1
S[3:end,3:end] .= dump(S[3:end,3:end])

Nothing wrong with loops in Julia:

using SparseArrays
k = 0
newvals = [7,8,9]
ii,jj,vv = findnz(S)
for (i,j) in zip(ii,jj)
  if i >= 3 && j >= 3
    S[i,j] = newvals[k+=1]
  end
end
S

If you’re just trying to mutate the diagonal of the matrix:

newvals = [7,8,9]
for (i,v) in zip(3:5,newvals)
  S[i,i] = v
end
1 Like

Thank you for your note.
Actually, my main issue is how to integrate this mutation inside the function dump (which exist in my main code)?

using SparseArrays
function dump(R,startat)
	k = 0
	newvals = 7:9
	ii,jj,vv = findnz(R)
	for (i,j) in zip(ii,jj)
		if i >= startat && j >= startat
			R[i,j] = newvals[k+=1]
		end
	end
	return R
end

# S will be mutated since dump mutates the first argument
dump(S,3)

Since dump mutates its argument, we’d typically name it dump! instead.

1 Like

Perfect. There are two things:
1- Since dump(S,3) is not assigned to any parameters, so what will return R do?
2- I deleted the return R and the results are not effected.

The return R allows you to write X = dump(S,3) and now get X === S. Often times this can be useful, but it’s certainly not necessary. Without this, dump returns the value of it’s last line. Since the last “line” here is a loop body, it will return nothing if no explicit return is written.

1 Like

Thank you very much for your details!