A few more points:
- Simon's suggestion of checking the types is spot on; if one is a matrix of
Any, for example, you're doomed to a slower code path. If the types are
Array{Float64,2} and Array{Float64,1}, that's not the problem.
- It should be slightly faster if you change your parentheses,
shrt*(diagm(expr)*shrt'). Jutho's scale suggestion would be even better.
- Try the profiler (see the docs). If you're running julia 0.4 and see it's
spending a lot of time in generic_matmatmul, do a git pull and rebuild---it
should fix the problem.
--Tim
On Wednesday, January 21, 2015 01:05:06 AM Simon Byrne wrote:
> As Jutho said, this shouldn't happen, but is difficult to diagnose without
> further information. What are the types of shrt and expr? (these can be
> found using the typeof function).
>
> Simon
>
> On Wednesday, 21 January 2015 07:59:34 UTC, Jutho wrote:
> > Not sure what is causing the slowness, but you could avoid creating a
> > diagonal matrix and then doing the matrix multiplication with diagm(expr)
> > which will be treated as a full matrix.
> > Instead of shrt*diagm(expr) which is interpreted as the multiplication of
> > two full matrices, try scale(shrt,expr) .
> >
> > Op woensdag 21 januari 2015 07:56:19 UTC+1 schreef Micah McClimans:
> >> I'm running into trouble with a line of matrix multiplication going very
> >> slowly in one of my programs. The line I'm looking into is:
> >> objectivematrix=shrt*diagm(expr)*(shrt')
> >> where shrt is 12,000x600 and expr is 600 long. This line takes several
> >> HOURS to run, on a computer that can run
> >>
> >> k=rand(12000,12000)
> >> k3=k*k*k
> >>
> >> in under a minute. I've tried devectorizing the line into the following
> >> loop (shrt is block-diagonal with each block ONevecs and -ONevecs
> >> respectively, so I split the loop in half)
> >>
> >> objectivematrix=zeros(2*size(ONevecs,1),2*size(ONevecs,1))
> >> for i in 1:size(ONevecs,1)
> >>
> >> print(i)
> >> for j in 1:size(ONevecs,1)
> >>
> >> for k in 1:size(ONevecs,2)
> >> objectivematrix[i,j]+=ONevecs[i,k]*ONevecs[j,k]*expr[k]
> >> end
> >>
> >> end
> >>
> >> end
> >> for i in 1:size(ONevecs,1)
> >>
> >> print(i)
> >> for j in 1:size(ONevecs,1)
> >>
> >> for k in 1:size(ONevecs,2)
> >>
> >> objectivematrix[i+size(ONevecs,1),j+size(ONevecs,1)]+=ONevecs[i,k]*ONevec
> >> s[j,k]*expr[k+size(ONevecs,2)]>>
> >> end
> >>
> >> end
> >>
> >> end
> >>
> >> and this give a print out every couple seconds- it's faster than the
> >> matrix multiplication version, but not enough. Why is this taking so
> >> long?
> >> This should not be a hard operation.
