### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [April 24, 2021, 12:12pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/1 "2021-04-24T12:12:41Z")
I want to do an inplace add of a multiple of the identity to a matrix. I can get the efficiency I want with a view of the diagonal:
```julia
using LinearAlgebra
function f1!(A,s)
D = view(A, diagind(A, 0))
D .+= s;
end
```
```julia
julia> s=1e-5;
julia> A=randn(1000,1000); x=diag(A);
julia> f1!(A,s);
julia> norm(diag(A) - (x.+s))
0.0
julia> @btime f1!(A,s)
4.627 μs (3 allocations: 144 bytes)
```
For various reasons, I would prefer to use the `UniformScaling`. Is it possible to achieve the same efficiency with `UniformScaling`?
For reference:
```julia
function f2!(A,s)
A[:,:] += s*I;
end
```
```julia
julia> @btime f2!(A,s);
5.750 ms (8 allocations: 15.26 MiB)
```
As far as I understand, there are two problems with `f2!`:
- Extra memory allocation (associated parse time conversion of `+=`)
- The `+(A::Matrix,J::UniformScaling)` is a for-loop referenced below, rather than a blas level-1 `axpy`-call which is what we get in `f1!`.
>
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [April 24, 2021, 12:19pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/2 "2021-04-24T12:19:58Z")
I don’t know of a way to do this with `I`, but a loop works and is efficient (faster than your `f1!` on my machine):
```julia
function f3!(A,s)
m = min(size(A)...)
for i = 1:m; A[i,i] += s; end
return A
end
```
> [@jarl](#):
>
> The `+(A::Matrix,J::UniformScaling)` is a for-loop referenced below, rather than a blas level-1 `axpy` -call which is what we get in `f1!` .
`f1!` calls the `broadcast` machinery under the hood, not BLAS. It’s possible but a bit tricky to use a BLAS `axpy` function for this operation:
```julia
function f4!(A::Matrix{T}, s) where {T} m = min(size(A)...)
m = min(size(A)...)
incA = stride(A,1) + stride(A,2)
sr = Ref{T}(s)
GC.@preserve sr LinearAlgebra.BLAS.axpy!(m, one(T), Base.unsafe_convert(Ptr{T}, sr), 0, A, incA)
return A
end
```
but on my machine it’s only 5% faster than my looping implementation `f3!` for your 1000×1000 benchmark. And even this small performance difference goes away if I use `@inbounds` in my `f3!` loop.
---
### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [April 24, 2021, 12:57pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/3 "2021-04-24T12:57:05Z")
Thanks! Learning a lot. I thought `f1!` would be effectively be the same as `f4!`. So, the BLAS `incx=0`-trick for `vec .+ scalar` is never used in julia?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [April 24, 2021, 1:00pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/4 "2021-04-24T13:00:54Z")
> [@jarl](#):
>
> So, the BLAS `incx=0` -trick for `vec .+ scalar` is never used in julia?
Not as far as I know. (There’s little or no benefit to BLAS over a simple `for` loop for `axpy` anyway, especially for `vector .+ scalar`; compilers are good at `axpy`-like loops and there’s no possibility of fancy blocking hand optimizations like there is for BLAS-3 / matrix multiply / `gemm`, while loops are far more versatile in supporting more types etcetera.)
(As far as I know, BLAS `axpy` is not even used for `vector + vector` in Julia. There’s no point. Besides, if you care about performance you’re [much better off](https://julialang.org/blog/2017/01/moredots/#why_vectorized_code_is_not_as_fast_as_it_could_be) combining multiple vector operations into a single loop, or [using “dot fusion”](https://julialang.org/blog/2017/01/moredots/), than breaking your calculation up into a sequence of `axpy` and other elementary operations.)
In general, [performance optimization in Julia](https://docs.julialang.org/en/v1/manual/performance-tips/) doesn’t rely on “mining” the standard library in the hope of finding a “vectorized / built-in” function that does exactly what you want (unlike e.g. Matlab or Python). Properly written user code and loops are fast.
---
### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [April 24, 2021, 6:04pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/5 "2021-04-24T18:04:08Z")
Great. Thanks. That’s helpful in other places in my code.
I’ll leave this post open for a while since an efficient version involving `I` would be helpful.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 29, 2021, 3:31pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/6 "2021-04-29T15:31:15Z")
You can use `diagind` in order to do this inplace on the diagonal. It’s how we do it in OrdinaryDiffEq to make it GPU-compatible:
>
---
### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [April 30, 2021, 6:17am UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/7 "2021-04-30T06:17:55Z")
Thanks. Does that compile to the same as the use of `diagind` in `f1!` ?
Edit: Ah. Now I understand your point. It does make it “work” for `I`.
---
### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [May 4, 2021, 5:46pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/8 "2021-05-04T17:46:54Z")
Wouldn’t an extension of `mul!` be a natural place to put this functionality?
```julia
import LinearAlgebra.mul!
function mul!(X::StridedMatrix{T},a::Bool,B::UniformScaling{T},alpha::Bool,beta::Bool) where {T}
if (a & alpha & beta)
D = view(X, diagind(X)) # Or more efficient version
D .+= B.λ
else
easytoimplement()
end
return X
end
function f6!(A::StridedMatrix,s)
mul!(A,true,s*I,true,true);
end
```
```julia
julia> A=randn(1000,1000);
julia> @btime f6!(A,3.0);
3.732 μs (2 allocations: 80 bytes)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 5, 2021, 6:32am UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/9 "2021-05-05T06:32:36Z")
That’s probably the right way to add it.
---
### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [May 5, 2021, 9:54am UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/10 "2021-05-05T09:54:37Z")
Okay. I will add a PR eventually. It doesn’t really involve a multiplication, so `mul!` is not an obvious for users who don’t know how five-argument `mul!` can be used.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [May 5, 2021, 10:07am UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/11 "2021-05-05T10:07:07Z")
Yes, `mul!()` is a strange function for an addition method. I think I would look for this functionality in `axpby!()`:
[`axpby!()` in manual:](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.axpby!)
---
### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [May 5, 2021, 11:07am UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/12 "2021-05-05T11:07:30Z")
I agree. `axpby!` does have a heritage from blas (even the manual is referring to `BLAS.axpby!`) and this feature is far from blas. `mul!` is julia specific. I wonder if `axpby!` is really meant to be used as `add!` analogous to `mul!`.
---
### Author:  [@jarl](https://discourse.julialang.org/u/jarl)
#### Post date: [May 6, 2021, 11:09am UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/13 "2021-05-06T11:09:18Z")
Discussion can be continued here: [Five arg mul! for UniformScaling and improvement in exp! by jarlebring · Pull Request #40731 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/40731)
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [May 7, 2021, 1:01pm UTC](https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/14 "2021-05-07T13:01:41Z")
The elegant solution would be for this to work:
```julia
A .+= (s*I)
```
But it doesn’t.
There is something blocking this, but I’m not sure what: [Broadcasting UniformScaling Operations · Issue #23197 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/23197)