I'm trying to create a sparse square matrix in Matlab through a mex function (written in Fortran). I want something like A = sparse(I,J,K)
. My triplets look like this, there are repetitions among the entries
femi = [1, 2, 3, 2, 2, 4, 5, 5, 4, 6, 6, 5, 5, 2]
femj = [2, 2, 1, 1, 1, 3, 3, 6, 3, 1, 1, 2, 2, 4]
femk = [2, 1, 5, 4, 2, 4, 5, 7, 2, 1, 6, 2, 1, 4]
I've written a rough piece of code, it works for small matrix dimensions, but it's much slower than the intrinsic Matlab's sparse
. Since I have almost no background in coding, I don't know what I'm doing wrong (wrong way to allocate variables? too many do loops?). Any help is appreciated. Thank you. This is the mex computational subroutine. It returns the pr, ir, jc indices array to give to the sparse matrix
subroutine new_sparse(femi, femj, femk, pr, ir, jc, n, m)
implicit none
intrinsic:: SUM, COUNT, ANY
integer :: i, j, k, n, indjc, m
real*8 :: femi(n), femj(n), femk(n)
real*8 :: pr(n)
integer :: ir(n),jc(m+1)
logical :: indices(n)
indices = .false.
k = 1
indjc = 0
jc(1) = 0
do j=1,m
do i =1,m
indices = [femi==i .and. femj==j]
if (ANY(indices .eqv. .true.)) then
ir(k) = i-1
pr(k) = SUM(femk, indices)
k = k+1
indjc = indjc + 1
end if
end do
if (indjc/=0) then
jc(j+1) = jc(j) + indjc
indjc = 0
else
jc(j+1) = jc(j)
end if
end do
return
end
Edit:
As suggested by users @jack and @veryreverie in the comments below, it's better to sort directly femi
, femj
and femk
. I guess that ranking/sorting femi
first (and sorting femj
and femk
according to femi
) and then ranking/sorting femj
(and sorting femi
and femk
according to femj
) provides the desired result. The only thing left is to deal with duplicates.
Edit #2 :
I translated line by line the serialized version of the C code by Engblom and Lukarksi . This document explains very clearly their reasoning and I think it's useful for beginners like me. However, due to my inexperience, I was unable to translate the parallelized version of the code. Maybe that prompts another question.
subroutine new_sparse(ir, jcS, pr, MatI, MatJ, MatK, n, m)
! use omp_lib
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer, intent(in) :: n, m
real(dp), intent(in) :: MatK(n), MatI(n), MatJ(n)
! integer*8, intent(out) :: nnew
integer :: i, k, col, row, c, r !, nthreads
integer :: hcol(m+1), jcS(m+1), jrS(m+1)
integer :: ixijs, irank(n), rank(n)
real*8 :: pr(*)
integer :: ir(*)
hcol = 0
jcS = 0
jrS = 0
do i = 1,n
jrS(MatI(i)+1) = jrS(MatI(i)+1)+1
end do
do r = 2,m+1
jrS(r) = jrS(r) + jrS(r-1)
end do
do i = 1,n
rank(jrS(MatI(i))+1) = i
jrS(MatI(i)) = jrS(MatI(i)) + 1
end do
k = 1
do row = 1,m
do i = k , jrS(row)
ixijs = rank(i)
col = MatJ(ixijs)
if (hcol(col) < row) then
hcol(col) = row
jcS(col+1) = jcS(col+1)+1
end if
irank(ixijs) = jcS(col+1)
k = k+1
end do
end do
do c = 2,m+1
jcS(c) = jcS(c) + jcS(c-1)
end do
do i = 1,n
irank(i) = irank(i) + jcS(MatJ(i))
end do
ir(irank) = MatI-1
do i = 1,n
pr(irank(i)) = pr(irank(i)) + MatK(i)
end do
return
end
This should work:
This writes:
The bottleneck in your algorithm comes from the instructions
indices = [femi==i .and. femj==j]
,any(indices .eqv. .true.)
andsum(femk, indices)
. These all takeO(n)
operations, and as these are within a double loop the overall cost of the subroutine isO(m^2*n)
.My algorithm works in two stages. The first stage, the
do b=jc(j),jc(j+1)-1
loop, compares each element in the input with each element in the matching column of the output, for a maximum cost ofO(mn)
operations. If the input element is found in the output, then the value is updated and nothing more needs to be done.If the input element is not found in the output, then it needs to be added to the output. This is handled by the second stage, the code after the
do b...
loop. Since this needs to move the output elements in order to make space for the new element, this stage has a maximum ofO(n'^2)
operations, wheren'
is the number of unique elements in the input, which should satisfyn'<=n
andn'<<m^2
for a sparse matrix.My algorithm should run a lot faster for large
m
andn
, but it certainly has a lot of scope for improvement. I suspect it is worth using an intermediate data structure for storingir
andpr
, so that new elements can be inserted without having to re-arrange all the elements to do so.