In [1]:
include("./ModuleB.jl")
using .ModuleB
You can see the availabe images at http://timholy.github.io/TestImages.jl/.
After installing the package with Pkg.add("TestImages"), the images are located in the directory .julia/packages/TestImages/..../images.
In [4]:
# import Pkg; Pkg.add("TestImages")
In [5]:
using TestImages, Images
In [6]:
varinfo(TestImages)
Out[6]:
In [10]:
img=testimage("lighthouse")
Out[10]:
In [11]:
# For example
img[1,1]
Out[11]:
In [12]:
show(img[1,1])
In [13]:
# Separate the image into R, G and B components
channels=channelview(img)
Out[13]:
In [14]:
R=channels[1,:,:]
G=channels[2,:,:]
B=channels[3,:,:]
Out[14]:
In [16]:
# Let us take a look at the red component
colorview(Gray,G)
Out[16]:
In [17]:
using LinearAlgebra
In [18]:
# Compute the SVD of each component
UR,σR,VR=svd(R)
UG,σG,VG=svd(G)
UB,σB,VB=svd(B);
In [19]:
# Plot the singular values
using Plots
x=collect(1:minimum(size(img)))
plot(x,[σR σG σB])
Out[19]:
In [24]:
# Watch the compression quality
k=150
Rc=UR[:,1:k]*Diagonal(σR[1:k])*VR[:,1:k]'
Gc=UG[:,1:k]*Diagonal(σG[1:k])*VG[:,1:k]'
Bc=UB[:,1:k]*Diagonal(σB[1:k])*VB[:,1:k]'
colorview(RGB, Rc, Gc, Bc)
Out[24]:
In [25]:
# Part of the preamble of lapack.jl
const liblapack = Base.liblapack_name
import LinearAlgebra.BLAS.@blasfunc
# import ..LinAlg: BlasFloat, Char, BlasInt, LAPACKException,
# DimensionMismatch, SingularException, PosDefException, chkstride1, chksquare
import LinearAlgebra.BlasInt
function chklapackerror(ret::BlasInt)
if ret == 0
return
elseif ret < 0
throw(ArgumentError("invalid argument #$(-ret) to LAPACK call"))
else # ret > 0
throw(LAPACKException(ret))
end
end
Out[25]:
In [29]:
for (gesvj, elty) in
((:dgesvj_,:Float64),
(:sgesvj_,:Float32))
@eval begin
function gesvj!(A::Array{$elty})
joba='G'
jobu='U'
jobv='V'
m,n = size(A)
lda=m
sva=Array{$elty}(undef,n)
mv=1
V=Array{$elty}(undef,n,n)
ldv=n
lwork=max(6,m+n)
work = Array{$elty}(undef,lwork)
info = Array{BlasInt}(undef,1)
ccall((@blasfunc($gesvj), liblapack), Cvoid,
(Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
Ref{BlasInt}, Ref{BlasInt}),
joba, jobu, jobv,
m, n, A,
lda, sva, mv,
V, ldv, work,
lwork, info)
chklapackerror(info[])
A, sva, V
end
end
end
In [30]:
# Small, strongly scaled matrix
m=20
n=15
using Random
Random.seed!(421)
B=rand(m,n)
D=exp.(80*(rand(n).-0.5))
A=B*Diagonal(D);
In [31]:
U,σ,V=ModuleB.myJacobiR(A)
U₁,σ₁,V₁=gesvj!(copy(A))
(sort(σ,rev=true)-σ₁)./σ₁
Out[31]:
In [33]:
# Timing for large matrices
n=1000
A=rand(n,n)
@time svd(A);
@time gesvj!(A);
In [ ]: