Remove bounds-checking in Bidiagonal rdiv! (#1281)

This improves performance
```julia
julia> B = Bidiagonal(fill(10.0,400), fill(2.0,399), :L);

julia> x = UpperTriangular(fill(400.0, size(B)));

julia> @Btime $x/$B;
  193.799 μs (3 allocations: 1.22 MiB) # nightly
  163.389 μs (3 allocations: 1.22 MiB) # this PR
```
Also, seems to reduce TTFX slightly, but at the expense of allocations:
```julia
julia> T = Bidiagonal(ones(4), ones(3), :U); TM = Matrix(T); x = UpperTriangular(TM);

julia> @time x/T;
  0.125422 seconds (104.07 k allocations: 5.261 MiB, 99.95% compilation time) # nightly
  0.099386 seconds (113.82 k allocations: 5.748 MiB, 99.94% compilation time) # this PR
```

Author: jishnub

src/bidiag.jl

@@ -1334,27 +1334,27 @@ function _rdiv!(C::AbstractMatrix, A::AbstractMatrix, B::Bidiagonal)
     isnothing(zi) || throw(SingularException(zi))
 
     if B.uplo == 'L'
-        diagB = B.dv[n]
-        for i in 1:m
-            C[i,n] = A[i,n] / diagB
+        diagB = @inbounds B.dv[n]
+        for i in axes(A,1)
+            @inbounds C[i,n] = A[i,n] / diagB
         end
-        for j in n-1:-1:1
-            diagB = B.dv[j]
-            offdiagB = B.ev[j]
-            for i in 1:m
-                C[i,j] = (A[i,j] - C[i,j+1]*offdiagB)/diagB
+        for j in reverse(axes(A,2)[1:end-1]) # n-1:-1:1
+            diagB = @inbounds B.dv[j]
+            offdiagB = @inbounds B.ev[j]
+            for i in axes(A,1)
+                @inbounds C[i,j] = (A[i,j] - C[i,j+1]*offdiagB)/diagB
             end
         end
     else
-        diagB = B.dv[1]
-        for i in 1:m
-            C[i,1] = A[i,1] / diagB
+        diagB = @inbounds B.dv[1]
+        for i in axes(A,1)
+            @inbounds C[i,1] = A[i,1] / diagB
         end
-        for j in 2:n
-            diagB = B.dv[j]
-            offdiagB = B.ev[j-1]
-            for i = 1:m
-                C[i,j] = (A[i,j] - C[i,j-1]*offdiagB)/diagB
+        for j in axes(A,2)[2:end]
+            diagB = @inbounds B.dv[j]
+            offdiagB = @inbounds B.ev[j-1]
+            for i in axes(A,1)
+                @inbounds C[i,j] = (A[i,j] - C[i,j-1]*offdiagB)/diagB
             end
         end
     end