On Tuesday, September 15, 2015 at 11:02:26 PM UTC-4, Jack Poulson wrote:
>
> I believe that Tony is suggesting manually applying the sparse operator
> rather than explicitly constructing it and then applying it. This is a
> common (and significant) performance optimization when a sparse operator is
> only used once or twice, as the construction of the sparse matrix is often
> *more* expensive than applying it once.
>
> The answer to your question (since you are applying the sparse matrix from
> the right to a row vector) would be a loop of the form: forall nonzero
> triplets (i,j,value), y[j] += x[i] * value.
>
> On Tuesday, September 15, 2015 at 8:52:46 AM UTC-7, Frank Kampas wrote:
>>
>> Could you post a link to the part of the documentation that describes how
>> to do that?
>>
>> On Tuesday, September 15, 2015 at 3:53:11 AM UTC-4, Tony Kelman wrote:
>>>
>>> Instead of constructing a sparse matrix in the inner loop it would be
>>> more efficient to write an in place stencil kernel function to perform the
>>> equivalent operation.
>>
>>
This is what I'm doing now. It's a little faster.
for i = 1:n
for j = i:n
if i == j
lhs[rctr,i] = -2 * start[i]
lhs[rctr,n+i] = -2 * start[n+i]
lhs[rctr,2n+1] = -2 * (radii[i] - start[2n+1])
rhs[rctr] = -start[i]^2 - start[n+i]^2 + start[2n+1]^2-
radii[i]^2
else
lhs[rctr,i] = 2 * (start[i] - start[j])
lhs[rctr,j] = 2 * (start[j] - start[i])
lhs[rctr,n+i] = 2 * (start[n+i]-start[n+j])
lhs[rctr,n+j] = 2 *(start[n+j] - start[n+i])
rhs[rctr] = (start[i]-start[j])^2 +
(start[n+i]-start[n+j])^2 + (radii[i]+radii[j])^2
end
rctr += 1
end
end
