← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to test/linalg.jl

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-02-11 03:51:26 UTC
  • mfrom: (1.1.4)
  • Revision ID: package-import@ubuntu.com-20130211035126-hap464pbhd97wjbl
Tags: 0.1~20130211.git86fbe98d-1
* New upstream snapshot.
* debian/control:
   + add git to Recommends (for Julia package manager)
   + remove dependencies on libglpk-dev (it moved to its own package)
   + add explicit dependency on libgmp10 (there is no more a wrapper)
* fix-clean-rules.patch: remove patch, applied upstream
* gsvddense_blasint.patch: new patch
* Refresh other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/linalg.jl
4
4
b     = rand(n)
5
5
for elty in (Float32, Float64, Complex64, Complex128)
6
6
        a     = convert(Matrix{elty}, a)
7
 
        asym  = a' + a                          # symmetric indefinite
8
 
        apd   = a'*a                            # symmetric positive-definite
 
7
        asym  = a' + a                  # symmetric indefinite
 
8
        apd   = a'*a                    # symmetric positive-definite
9
9
        b     = convert(Vector{elty}, b)
10
10
        
11
 
        capd  = chold(apd)                      # upper Cholesky factor
 
11
        capd  = cholfact(apd)              # upper Cholesky factor
12
12
        r     = factors(capd)
13
13
        @assert_approx_eq r'*r apd
14
14
        @assert_approx_eq b apd * (capd\b)
15
15
        @assert_approx_eq apd * inv(capd) eye(elty, n)
16
 
        @assert_approx_eq a*(capd\(a'*b)) b          # least squares soln for square a
 
16
        @assert_approx_eq a*(capd\(a'*b)) b # least squares soln for square a
17
17
        @assert_approx_eq det(capd) det(apd)
18
18
 
19
 
        l     = factors(chold(apd, 'L'))      # lower Cholesky factor
 
19
        l     = factors(cholfact(apd, 'L')) # lower Cholesky factor
20
20
        @assert_approx_eq l*l' apd
21
21
 
22
 
        cpapd = cholpd(apd)                     # pivoted Choleksy decomposition
 
22
        cpapd = cholpfact(apd)           # pivoted Choleksy decomposition
23
23
        @test rank(cpapd) == n
24
 
        @test all(diff(diag(real(cpapd.LR))).<=0.)  # diagonal show be non-increasing
 
24
        @test all(diff(diag(real(cpapd.LR))).<=0.) # diagonal should be non-increasing
25
25
        @assert_approx_eq b apd * (cpapd\b)
26
26
        @assert_approx_eq apd * inv(cpapd) eye(elty, n)
27
27
 
28
 
        bc1   = BunchKaufman(asym)              # Bunch-Kaufman factor of indefinite matrix
 
28
        bc1   = BunchKaufman(asym) # Bunch-Kaufman factor of indefinite matrix
29
29
        @assert_approx_eq inv(bc1) * asym eye(elty, n)
30
30
        @assert_approx_eq asym * (bc1\b) b
31
 
        bc2   = BunchKaufman(apd)               # Bunch-Kaufman factors of a pos-def matrix
 
31
        bc2   = BunchKaufman(apd) # Bunch-Kaufman factors of a pos-def matrix
32
32
        @assert_approx_eq inv(bc2) * apd eye(elty, n)
33
33
        @assert_approx_eq apd * (bc2\b) b
34
34
 
35
 
        lua   = lud(a)                          # LU decomposition
 
35
        lua   = lufact(a)                  # LU decomposition
36
36
        l,u,p = lu(a)
37
37
        L,U,P = factors(lua)
38
38
        @test l == L && u == U && p == P
41
41
        @assert_approx_eq a * inv(lua) eye(elty, n)
42
42
        @assert_approx_eq a*(lua\b) b
43
43
 
44
 
        qra   = qrd(a)                          # QR decomposition
 
44
        qra   = qrfact(a)                  # QR decomposition
45
45
        q,r   = factors(qra)
46
46
        @assert_approx_eq q'*q eye(elty, n)
47
47
        @assert_approx_eq q*q' eye(elty, n)
48
48
        Q,R   = qr(a)
49
49
        @test q == Q && r == R
50
50
        @assert_approx_eq q*r a
51
 
        @assert_approx_eq qra*b Q*b
52
 
        @assert_approx_eq qra'*b Q'*b
 
51
        @assert_approx_eq qmulQR(qra,b) Q*b
 
52
        @assert_approx_eq qTmulQR(qra,b) Q'*b
53
53
        @assert_approx_eq a*(qra\b) b
54
54
 
55
 
        qrpa  = qrpd(a)                         # pivoted QR decomposition
 
55
        qrpa  = qrpfact(a)                 # pivoted QR decomposition
56
56
        q,r,p = factors(qrpa)
57
57
        @assert_approx_eq q'*q eye(elty, n)
58
58
        @assert_approx_eq q*q' eye(elty, n)
62
62
        @assert_approx_eq q*r[:,invperm(p)] a
63
63
        @assert_approx_eq a*(qrpa\b) b
64
64
 
65
 
        d,v   = eig(asym)                       # symmetric eigen-decomposition
 
65
        d,v   = eig(asym)              # symmetric eigen-decomposition
66
66
        @assert_approx_eq asym*v[:,1] d[1]*v[:,1]
67
67
        @assert_approx_eq v*diagmm(d,v') asym
68
68
 
69
 
        d,v   = eig(a)                          # non-symmetric eigen decomposition
 
69
        d,v   = eig(a)                 # non-symmetric eigen decomposition
70
70
        for i in 1:size(a,2) @assert_approx_eq a*v[:,i] d[i]*v[:,i] end
71
71
    
72
 
        u, q, v = schur(a)                         # Schur
 
72
        u, q, v = schur(a)             # Schur
73
73
        @assert_approx_eq q*u*q' a
74
74
        @assert_approx_eq sort(real(v)) sort(real(d))
75
75
        @assert_approx_eq sort(imag(v)) sort(imag(d))
76
76
        @test istriu(u) || isreal(a)
77
77
 
78
 
        u,s,vt = svdt(a)                        # singular value decomposition
 
78
        u,s,vt = svdt(a)                # singular value decomposition
79
79
        @assert_approx_eq u*diagmm(s,vt) a
80
80
    
 
81
        gsvd = svd(a,a[1:5,:])         # Generalized svd
 
82
        @assert_approx_eq gsvd[1]*gsvd[4]*gsvd[6]*gsvd[3]' a
 
83
        @assert_approx_eq gsvd[2]*gsvd[5]*gsvd[6]*gsvd[3]' a[1:5,:]
 
84
 
81
85
        x = a \ b
82
86
        @assert_approx_eq a*x b
83
87
    
87
91
        x = tril(a)\b
88
92
        @assert_approx_eq tril(a)*x b
89
93
 
90
 
        # Test null
 
94
                                        # Test null
91
95
        a15null = null(a[:,1:5]')
92
96
        @assert_approx_eq_eps norm(a[:,1:5]'a15null) zero(elty) n*eps(one(elty))
93
97
        @assert_approx_eq_eps norm(a15null'a[:,1:5]) zero(elty) n*eps(one(elty))
94
98
        @test size(null(b), 2) == 0
95
99
 
96
 
        # Test pinv
 
100
                                        # Test pinv
97
101
        pinva15 = pinv(a[:,1:5])
98
102
        @assert_approx_eq a[:,1:5]*pinva15*a[:,1:5] a[:,1:5]
99
103
        @assert_approx_eq pinva15*a[:,1:5]*pinva15 pinva15
100
104
    
101
 
        # Complex complex rhs real lhs
 
105
                                        # Complex vector rhs
102
106
        x = a\complex(b)
103
107
        @assert_approx_eq a*x complex(b)
104
108
 
105
 
        # Test cond
 
109
                                        # Test cond
106
110
        @assert_approx_eq_eps cond(a, 1) 4.837320054554436e+02 0.01
107
111
        @assert_approx_eq_eps cond(a, 2) 1.960057871514615e+02 0.01
108
112
        @assert_approx_eq_eps cond(a, Inf) 3.757017682707787e+02 0.01
109
113
        @assert_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01
110
114
 
111
 
        # Matrix square root
 
115
                                        # Matrix square root
112
116
        asq = sqrtm(a)
113
117
        @assert_approx_eq asq*asq a
114
118
end
 
119
 
 
120
## Least squares solutions
115
121
a = [ones(20) 1:20 1:20]
116
122
b = reshape(eye(8, 5), 20, 2)
117
123
for elty in (Float32, Float64, Complex64, Complex128)
118
124
        a = convert(Matrix{elty}, a)
119
125
        b = convert(Matrix{elty}, b)
120
126
 
121
 
        # Matrix and vector multiplication
122
 
        @assert_approx_eq b'b convert(Matrix{elty}, [3 0; 0 2])
123
 
        @assert_approx_eq b'b[:,1] convert(Vector{elty}, [3, 0])
124
 
        @assert_approx_eq dot(b[:,1], b[:,1]) convert(elty, 3.0)
125
 
 
126
 
        # Least squares        
127
 
        x = a[:,1:2]\b[:,1]                             # Vector rhs
 
127
        x = a[:,1:2]\b[:,1]             # Vector rhs
128
128
        @assert_approx_eq ((a[:,1:2]*x-b[:,1])'*(a[:,1:2]*x-b[:,1]))[1] convert(elty, 2.546616541353384)
129
129
    
130
 
        x = a[:,1:2]\b                                  # Matrix rhs
 
130
        x = a[:,1:2]\b                  # Matrix rhs
131
131
        @assert_approx_eq det((a[:,1:2]*x-b)'*(a[:,1:2]*x-b)) convert(elty, 4.437969924812031)
132
132
    
133
 
        x = a\b                            # Rank deficient
 
133
        x = a\b                         # Rank deficient
134
134
        @assert_approx_eq det((a*x-b)'*(a*x-b)) convert(elty, 4.437969924812031)
135
135
 
136
 
        x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1]) # Underdetermined minimum norm
 
136
                                        # Underdetermined minimum norm
 
137
        x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1]) 
137
138
        @assert_approx_eq x convert(Vector{elty}, [1, 0.5, -0.5])
138
139
 
139
 
        # symmetric, positive definite
 
140
                                        # symmetric, positive definite
140
141
        @assert_approx_eq inv(convert(Matrix{elty}, [6. 2; 2 1])) convert(Matrix{elty}, [0.5 -1; -1 3])
141
 
        # symmetric, negative definite
 
142
                                        # symmetric, negative definite
142
143
        @assert_approx_eq inv(convert(Matrix{elty}, [1. 2; 2 1])) convert(Matrix{elty}, [-1. 2; 2 -1]/3)               
143
144
end
144
145
 
145
146
## Test Julia fallbacks to BLAS routines
146
 
# matrices with zero dimensions
 
147
                                        # matrices with zero dimensions
147
148
@test ones(0,5)*ones(5,3) == zeros(0,3)
148
149
@test ones(3,5)*ones(5,0) == zeros(3,0)
149
150
@test ones(3,0)*ones(0,4) == zeros(3,4)
151
152
@test ones(0,0)*ones(0,4) == zeros(0,4)
152
153
@test ones(3,0)*ones(0,0) == zeros(3,0)
153
154
@test ones(0,0)*ones(0,0) == zeros(0,0)
154
 
# 2x2
 
155
                                        # 2x2
155
156
A = [1 2; 3 4]
156
157
B = [5 6; 7 8]
157
158
@test A*B == [19 22; 43 50]
164
165
@test Ac_mul_B(Ai, Bi) == [68.5-12im 57.5-28im; 88-3im 76.5-25im]
165
166
@test A_mul_Bc(Ai, Bi) == [64.5+5.5im 43+31.5im; 104-18.5im 80.5+31.5im]
166
167
@test Ac_mul_Bc(Ai, Bi) == [-28.25-66im 9.75-58im; -26-89im 21-73im]
167
 
# 3x3
 
168
                                        # 3x3
168
169
A = [1 2 3; 4 5 6; 7 8 9]-5
169
170
B = [1 0 5; 6 -10 3; 2 -4 -1]
170
171
@test A*B == [-26 38 -27; 1 -4 -6; 28 -46 15]
177
178
@test Ac_mul_B(Ai, Bi) == [-21+2im -1.75+49im -51.25+19.5im; 25.5+56.5im -7-35.5im 22+35.5im; -3+12im -32.25+43im -34.75-2.5im]
178
179
@test A_mul_Bc(Ai, Bi) == [-20.25+15.5im -28.75-54.5im 22.25+68.5im; -12.25+13im -15.5+75im -23+27im; 18.25+im 1.5+94.5im -27-54.5im]
179
180
@test Ac_mul_Bc(Ai, Bi) == [1+2im 20.75+9im -44.75+42im; 19.5+17.5im -54-36.5im 51-14.5im; 13+7.5im 11.25+31.5im -43.25-14.5im]
180
 
# Generic integer matrix multiplication
 
181
                                        # Generic integer matrix multiplication
181
182
A = [1 2 3; 4 5 6] - 3
182
183
B = [2 -2; 3 -5; -4 7]
183
184
@test A*B == [-7 9; -4 9]
189
190
B = rand(1:20, 5, 5) - 10
190
191
@test At_mul_B(A, B) == A'*B
191
192
@test A_mul_Bt(A, B) == A*B'
192
 
# Preallocated
 
193
                                        # Preallocated
193
194
C = Array(Int, size(A, 1), size(B, 2))
194
195
@test A_mul_B(C, A, B) == A*B
195
196
@test At_mul_B(C, A, B) == A'*B
196
197
@test A_mul_Bt(C, A, B) == A*B'
197
198
@test At_mul_Bt(C, A, B) == A'*B'
198
 
# matrix algebra with subarrays of floats (stride != 1)
 
199
                                        # matrix algebra with subarrays of floats (stride != 1)
199
200
A = reshape(float64(1:20),5,4)
200
201
Aref = A[1:2:end,1:2:end]
201
202
Asub = sub(A, 1:2:5, 1:2:4)
208
209
Asub = sub(Ai, 1:2:5, 1:2:4)
209
210
@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref)
210
211
@test A_mul_Bc(Asub, Asub) == A_mul_Bc(Aref, Aref)
211
 
# syrk & herk
 
212
                                        # syrk & herk
212
213
A = reshape(1:1503, 501, 3)-750.0
213
214
res = float64([135228751 9979252 -115270247; 9979252 10481254 10983256; -115270247 10983256 137236759])
214
215
@test At_mul_B(A, A) == res
223
224
Aref = Ai[1:2:2*cutoff, 1:3]
224
225
@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref)
225
226
 
226
 
# Matrix exponential and Hessenberg
 
227
                                        # Matrix exponential
227
228
for elty in (Float32, Float64, Complex64, Complex128)
228
229
        A1  = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4])
229
230
        eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181  127.781085523182;
247
248
        0.135335281175235 0.406005843524598 0.541341126763207]')
248
249
        @assert_approx_eq expm(A3) eA3
249
250
 
250
 
        # Hessenberg
 
251
                                        # Hessenberg
251
252
        @assert_approx_eq hess(A1) convert(Matrix{elty}, 
252
253
                        [4.000000000000000  -1.414213562373094  -1.414213562373095
253
254
                        -1.414213562373095   4.999999999999996  -0.000000000000000
254
255
                                         0  -0.000000000000002   3.000000000000000])
255
256
end
256
257
 
257
 
# matmul for types w/o sizeof (issue #1282)
 
258
                                        # matmul for types w/o sizeof (issue #1282)
258
259
A = Array(ComplexPair{Int},10,10)
259
260
A[:] = complex(1,1)
260
261
A2 = A^2
261
262
@test A2[1,1] == 20im
262
263
 
263
 
# basic tridiagonal operations
 
264
                                        # basic tridiagonal operations
264
265
n = 5
265
266
d = 1 + rand(n)
266
267
dl = -rand(n-1)
267
268
du = -rand(n-1)
268
269
v = randn(n)
269
270
B = randn(n,2)
270
 
# Woodbury
 
271
                                        # Woodbury
271
272
U = randn(n,2)
272
273
V = randn(2,n)
273
274
C = randn(2,2)
285
286
            F[i+1,i] = dl[i]
286
287
        end
287
288
        @test full(T) == F
288
 
 
289
 
        # tridiagonal linear algebra
 
289
                                        # tridiagonal linear algebra
290
290
        v = convert(Vector{elty}, v)
291
291
        @assert_approx_eq T*v F*v
292
292
        invFv = F\v
294
294
        @assert_approx_eq solve(T,v) invFv
295
295
        B = convert(Matrix{elty}, B)
296
296
        @assert_approx_eq solve(T, B) F\B
297
 
        Tlu = lud(T)
 
297
        Tlu = lufact(T)
298
298
        x = Tlu\v
299
299
        @assert_approx_eq x invFv
300
300
        @assert_approx_eq det(T) det(F)
301
 
 
302
 
        # symmetric tridiagonal
 
301
                                        # symmetric tridiagonal
303
302
        Ts = SymTridiagonal(d, dl)
304
303
        Fs = full(Ts)
305
304
        invFsv = Fs\v
306
305
        Tldlt = ldltd(Ts)
307
306
        x = Tldlt\v
308
307
        @assert_approx_eq x invFsv
309
 
 
310
 
        # eigenvalues/eigenvectors of symmetric tridiagonal
 
308
                                        # eigenvalues/eigenvectors of symmetric tridiagonal
311
309
        if elty === Float32 || elty === Float64
312
310
            DT, VT = eig(Ts)
313
311
            D, Vecs = eig(Fs)
314
312
            @assert_approx_eq DT D
315
313
            @assert_approx_eq abs(VT'Vecs) eye(elty, n)
316
314
        end
317
 
        
318
 
        # Woodbury
 
315
                                        # Woodbury
319
316
        U = convert(Matrix{elty}, U)
320
317
        V = convert(Matrix{elty}, V)
321
318
        C = convert(Matrix{elty}, C)
355
352
        @assert_approx_eq_eps det(ones(elty, 3,3)) zero(elty) 3*eps(one(elty))
356
353
end
357
354
 
358
 
# LAPACK tests
 
355
                                        # LAPACK tests
359
356
Ainit = randn(5,5)
360
357
for elty in (Float32, Float64, Complex64, Complex128)
361
 
        # syevr!
 
358
                                        # syevr!
362
359
        A = convert(Array{elty, 2}, Ainit)
363
360
        Asym = A'A
364
361
        Z = Array(elty, 5, 5)