forked from JuliaLang/julia
-
Notifications
You must be signed in to change notification settings - Fork 0
/
arpack.jl
106 lines (89 loc) · 3.59 KB
/
arpack.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
begin
srand(1234)
local n,a,asym,b,bsym,d,v
n = 10
areal = sprandn(n,n,0.4)
breal = sprandn(n,n,0.4)
acmplx = complex(sprandn(n,n,0.4), sprandn(n,n,0.4))
bcmplx = complex(sprandn(n,n,0.4), sprandn(n,n,0.4))
testtol = 1e-6
for elty in (Float64, Complex128)
if elty == Complex64 || elty == Complex128
a = acmplx
b = bcmplx
else
a = areal
b = breal
end
a = convert(SparseMatrixCSC{elty}, a)
asym = a' + a # symmetric indefinite
apd = a'*a # symmetric positive-definite
b = convert(SparseMatrixCSC{elty}, b)
bsym = b' + b
bpd = b'*b
(d,v) = eigs(a, nev=3)
@test_approx_eq a*v[:,2] d[2]*v[:,2]
@test norm(v) > testtol # eigenvectors cannot be null vectors
# (d,v) = eigs(a, b, nev=3, tol=1e-8) not handled yet
# @test_approx_eq_eps a*v[:,2] d[2]*b*v[:,2] testtol
# @test norm(v) > testtol # eigenvectors cannot be null vectors
(d,v) = eigs(asym, nev=3)
@test_approx_eq asym*v[:,1] d[1]*v[:,1]
@test_approx_eq eigs(asym; nev=1, sigma=d[3])[1][1] d[3]
@test norm(v) > testtol # eigenvectors cannot be null vectors
(d,v) = eigs(apd, nev=3)
@test_approx_eq apd*v[:,3] d[3]*v[:,3]
@test_approx_eq eigs(apd; nev=1, sigma=d[3])[1][1] d[3]
(d,v) = eigs(apd, bpd, nev=3, tol=1e-8)
@test_approx_eq_eps apd*v[:,2] d[2]*bpd*v[:,2] testtol
@test norm(v) > testtol # eigenvectors cannot be null vectors
# test (shift-and-)invert mode
(d,v) = eigs(apd, nev=3, sigma=0)
@test_approx_eq apd*v[:,3] d[3]*v[:,3]
@test norm(v) > testtol # eigenvectors cannot be null vectors
(d,v) = eigs(apd, bpd, nev=3, sigma=0, tol=1e-8)
@test_approx_eq_eps apd*v[:,1] d[1]*bpd*v[:,1] testtol
@test norm(v) > testtol # eigenvectors cannot be null vectors
end
end
# Example from Quantum Information Theory
import Base: size, issym, ishermitian
type CPM{T<:Base.LinAlg.BlasFloat}<:AbstractMatrix{T} # completely positive map
kraus::Array{T,3} # kraus operator representation
end
size(Phi::CPM)=(size(Phi.kraus,1)^2,size(Phi.kraus,3)^2)
issym(Phi::CPM)=false
ishermitian(Phi::CPM)=false
function *{T<:Base.LinAlg.BlasFloat}(Phi::CPM{T},rho::Vector{T})
rho=reshape(rho,(size(Phi.kraus,3),size(Phi.kraus,3)))
rho2=zeros(T,(size(Phi.kraus,1),size(Phi.kraus,1)))
for s=1:size(Phi.kraus,2)
As=slice(Phi.kraus,:,s,:)
rho2+=As*rho*As'
end
return reshape(rho2,(size(Phi.kraus,1)^2,))
end
# Generate random isometry
(Q,R)=qr(randn(100,50))
Q=reshape(Q,(50,2,50))
# Construct trace-preserving completely positive map from this
Phi=CPM(Q)
(d,v,nconv,numiter,numop,resid) = eigs(Phi,nev=1,which="LM")
# Properties: largest eigenvalue should be 1, largest eigenvector, when reshaped as matrix
# should be a Hermitian positive definite matrix (up to an arbitrary phase)
@test_approx_eq d[1] 1. # largest eigenvalue should be 1.
v=reshape(v,(50,50)) # reshape to matrix
v/=trace(v) # factor out arbitrary phase
@test isapprox(vecnorm(imag(v)),0.) # it should be real
v=real(v)
# @test isapprox(vecnorm(v-v')/2,0.) # it should be Hermitian
# Since this fails sometimes (numerical precision error),this test is commented out
v=(v+v')/2
@test isposdef(v)
# Repeat with starting vector
(d2,v2,nconv2,numiter2,numop2,resid2) = eigs(Phi,nev=1,which="LM",v0=reshape(v,(2500,)))
v2=reshape(v2,(50,50))
v2/=trace(v2)
@test numiter2<numiter
@test_approx_eq v v2
@test_approx_eq eigs(speye(50), nev=10)[1] ones(10) #Issue 4246