Merge pull request #9 from ivanslapnicar/ivan/morematrices

Added more matrices ...

README.md

@@ -3,13 +3,70 @@
 A [Julia](http://julialang.org) package for working with special matrix types.
 
 This package extends the `Base.LinAlg` library with support for special
-matrices which are used in linear algebra.
+matrices which are used in linear algebra. Every special matrix has its own type.
+The full matrix is accessed by the command `full(A)`.
 
 [![Build Status](https://travis-ci.org/jiahao/SpecialMatrices.jl.svg)](https://travis-ci.org/jiahao/SpecialMatrices.jl) [![Coverage Status](https://img.shields.io/coveralls/jiahao/SpecialMatrices.jl.svg)](https://coveralls.io/r/jiahao/SpecialMatrices.jl)
 
-
 ## Currently supported special matrices
 
+## [`Cauchy`](http://en.wikipedia.org/wiki/Cauchy_matrix) matrix
+
+```julia
+julia> Cauchy([1,2,3],[3,4,5])
+3x3 Cauchy{Int64}:
+ 0.25      0.2       0.166667
+ 0.2       0.166667  0.142857
+ 0.166667  0.142857  0.125
+
+julia> Cauchy([1,2,3])
+3x3 Cauchy{Int64}:
+ 0.5       0.333333  0.25
+ 0.333333  0.25      0.2
+ 0.25      0.2       0.166667
+
+julia> Cauchy(pi)
+3x3 Cauchy{Float64}:
+ 0.5       0.333333  0.25
+ 0.333333  0.25      0.2
+ 0.25      0.2       0.166667
+```
+
+## `Circulant` matrix
+
+```julia
+julia> Circulant([1:4])
+4x4 Circulant{Int64}:
+ 1  4  3  2
+ 2  1  4  3
+ 3  2  1  4
+ 4  3  2  1
+```
+
+## [`Companion`](http://en.wikipedia.org/wiki/Companion_matrix) matrix
+
+```julia
+julia> A=Companion([3,2,1])
+3x3 Companion{Int64}:
+ 0  0  -3
+ 1  0  -2
+ 0  1  -1
+```
+Also, directly from a polynomial:
+
+```julia
+julia> using Polynomials
+
+julia> P=Poly([2,3,4,5])
+Poly(2 + 3x + 4x^2 + 5x^3)
+
+julia> C=Companion(P)
+3x3 Companion{Number}:
+ 0  0  -0.4
+ 1  0  -0.6
+ 0  1  -0.8
+```
+
 ## [`Frobenius`](http://en.wikipedia.org/wiki/Frobenius_matrix) matrix
 
 Example
@@ -63,16 +120,110 @@ julia> F*[10.0:10.0:60.0]
  150.0
 ```
 
-## [`Companion`](http://en.wikipedia.org/wiki/Companion_matrix) matrix
+## [`Hankel`](http://en.wikipedia.org/wiki/Hankel_matrix) matrix
+
+Input is a vector of odd length.
 
 ```julia
-julia> A=Companion([1,2,1])
-3x3 Companion{Int64}:
- 0  0  -1
- 1  0  -2
- 0  1  -1
+julia> Hankel([-4:4])
+5x5 Hankel{Int64}:
+ -4  -3  -2  -1  0
+ -3  -2  -1   0  1
+ -2  -1   0   1  2
+ -1   0   1   2  3
+  0   1   2   3  4
+```
+
+## [`Hilbert`](http://en.wikipedia.org/wiki/Hilbert_matrix) matrix
+
+```julia
+julia> A=Hilbert(5)
+Hilbert{Rational{Int64}}(5,5)
+
+julia> full(A)
+5x5 Array{Rational{Int64},2}:
+ 1//1  1//2  1//3  1//4  1//5
+ 1//2  1//3  1//4  1//5  1//6
+ 1//3  1//4  1//5  1//6  1//7
+ 1//4  1//5  1//6  1//7  1//8
+ 1//5  1//6  1//7  1//8  1//9
+
+julia> full(Hilbert(5))
+5x5 Array{Rational{Int64},2}:
+ 1//1  1//2  1//3  1//4  1//5
+ 1//2  1//3  1//4  1//5  1//6
+ 1//3  1//4  1//5  1//6  1//7
+ 1//4  1//5  1//6  1//7  1//8
+ 1//5  1//6  1//7  1//8  1//9
+```
+Inverses are also integer matrices:
+
+```julia
+julia> inv(A)
+5x5 Array{Rational{Int64},2}:
+    25//1    -300//1     1050//1    -1400//1     630//1
+  -300//1    4800//1   -18900//1    26880//1  -12600//1
+  1050//1  -18900//1    79380//1  -117600//1   56700//1
+ -1400//1   26880//1  -117600//1   179200//1  -88200//1
+   630//1  -12600//1    56700//1   -88200//1   44100//1
 ```
 
+## [`Kahan`](http://math.nist.gov/MatrixMarket/data/MMDELI/kahan/kahan.html) matrix
+
+```julia
+julia> A=Kahan(5,5,1,35)
+5x5 Kahan{Int64,Int64}:
+ 1.0  -0.540302  -0.540302  -0.540302  -0.540302
+ 0.0   0.841471  -0.454649  -0.454649  -0.454649
+ 0.0   0.0        0.708073  -0.382574  -0.382574
+ 0.0   0.0        0.0        0.595823  -0.321925
+ 0.0   0.0        0.0        0.0        0.501368
+
+julia> A=Kahan(5,3,0.5,0)
+5x3 Kahan{Float64,Int64}:
+ 1.0  -0.877583  -0.877583
+ 0.0   0.479426  -0.420735
+ 0.0   0.0        0.229849
+ 0.0   0.0        0.0
+ 0.0   0.0        0.0
+
+julia> A=Kahan(3,5,0.5,1e-3)
+3x5 Kahan{Float64,Float64}:
+ 1.0  -0.877583  -0.877583  -0.877583  -0.877583
+ 0.0   0.479426  -0.420735  -0.420735  -0.420735
+ 0.0   0.0        0.229849  -0.201711  -0.201711
+```
+
+For more details see [N. J. Higham (1987)][Higham87].
+
+[Higham87]: http://eprints.ma.man.ac.uk/695/01/covered/MIMS_ep2007_10.pdf "N. Higham, A Survey of Condition Number Estimation for Triangular Matrices, SIMAX, Vol. 29, No. 4, pp. 588, 1987"
+
+## `Riemann` matrix
+
+Riemann matrix is defined as `A = B[2:N+1, 2:N+1]`, where
+`B[i,j] = i-1` if `i` divides `j`, and `-1` otherwise.
+[Riemann hypothesis](http://en.wikipedia.org/wiki/Riemann_hypothesis) holds
+if and only if `det(A) = O( N! N^(-1/2+epsilon))` for every `epsilon > 0`.
+
+```julia
+julia> Riemann(7)
+7x7 Riemann{Int64}:
+  1  -1   1  -1   1  -1   1
+ -1   2  -1  -1   2  -1  -1
+ -1  -1   3  -1  -1  -1   3
+ -1  -1  -1   4  -1  -1  -1
+ -1  -1  -1  -1   5  -1  -1
+ -1  -1  -1  -1  -1   6  -1
+ -1  -1  -1  -1  -1  -1   7
+```
+
+For more details see [F. Roesler (1986)][Roesler1986].
+
+[Roesler1986]: http://www.sciencedirect.com/science/article/pii/0024379586902557 "Friedrich Roesler, Riemann's hypothesis as an eigenvalue problem, Linear Algebra and its Applications, Vol. 81, (1986)"
+
+
+
+
 ## `Strang` matrix
 
 A special `SymTridiagonal` matrix named after Gilbert Strang
@@ -88,19 +239,9 @@ julia> Strang(6)
   0.0   0.0   0.0   0.0  -1.0   2.0
 ```
 
-## `Hankel` matrix
-
-```julia
-julia> Hankel([-4:4])
-5x5 Hankel{Int64}:
- -4  -3  -2  -1  0
- -3  -2  -1   0  1
- -2  -1   0   1  2
- -1   0   1   2  3
-  0   1   2   3  4
-```
+## [`Toeplitz`](http://en.wikipedia.org/wiki/Toeplitz_matrix) matrix
 
-## `Toeplitz` matrix
+Input is a vector of odd length.
 
 ```julia
 julia> Toeplitz([-4:4])
@@ -111,31 +252,8 @@ julia> Toeplitz([-4:4])
  3   2   1   0  -1
  4   3   2   1   0
 ```
+## [`Vandermonde`](http://en.wikipedia.org/wiki/Vandermonde_matrix) matrix
 
-## `Circulant` matrix
-
-```julia
-julia> Circulant([1:4])
-4x4 Circulant{Int64}:
- 1  4  3  2
- 2  1  4  3
- 3  2  1  4
- 4  3  2  1
-```
-
-## `Hilbert` matrix
-
-```julia
-julia> full(Hilbert(5))
-5x5 Array{Rational{Int64},2}:
- 1//1  1//2  1//3  1//4  1//5
- 1//2  1//3  1//4  1//5  1//6
- 1//3  1//4  1//5  1//6  1//7
- 1//4  1//5  1//6  1//7  1//8
- 1//5  1//6  1//7  1//8  1//9
-```
-
-## `Vandermonde` matrix
 
 ```julia
 julia> Vandermonde([1:5])

REQUIRE

@@ -1,2 +1,2 @@
 julia 0.2-
-
+Polynomials 0.0.3

src/SpecialMatrices.jl

@@ -2,11 +2,14 @@ module SpecialMatrices
 
 import Base: A_mul_B!, full, getindex, inv, isassigned, size, *
 
+include("cauchy.jl") #Cauchy matrix
 include("companion.jl") #Companion matrix
 include("frobenius.jl") #Frobenius matrix
-include("strang.jl") #Strang matrix
 include("hankel.jl") #Hankel matrix
 include("hilbert.jl") #Hilbert matrix
+include("kahan.jl") #Kahan matrix
+include("riemann.jl") #Riemann matrix
+include("strang.jl") #Strang matrix
 include("toeplitz.jl") #Toeplitz matrix
 include("vandermonde.jl") #Vandermonde matrix
 

src/cauchy.jl

@@ -0,0 +1,29 @@
+#### Cauchy matrix
+
+export Cauchy
+
+immutable Cauchy{T} <: AbstractMatrix{T}
+    x::Vector{T} #
+    y::Vector{T} #
+end # immutable
+
+function Cauchy(k::Number)
+         Cauchy([1:k],[1:k])
+end
+
+function Cauchy(x::Vector)
+         Cauchy(x,x)
+end
+
+# Define its size
+
+size(A::Cauchy, dim::Integer) = length(A.x)
+size(A::Cauchy)= size(A,1), size(A,1)
+
+# Index into a Cauchy
+function getindex(A::Cauchy,i::Integer,j::Integer)
+    return 1.0/(A.x[i]+A.y[j])
+end # getindex
+
+# Dense version of Cauchy
+full(A::Cauchy) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)]

src/companion.jl

@@ -3,6 +3,20 @@ export Companion
 immutable Companion{T} <: AbstractArray{T, 2}
     c :: Vector{T}
 end
+#From polynomial
+
+using Polynomials
+#Generate companion matrix from a polynomial
+
+function Companion(P::Poly)
+   n = length(P)
+   c = Array(Number,n-1)
+   d=P.a[n]
+   for i=1:n-1
+       c[i]=P.a[i]/d
+   end
+   Companion(c)   
+end
 
 #Basic property computations
 size(C::Companion, r::Int) = (r==1 || r==2) ? length(C.c) : 

src/hilbert.jl

@@ -18,3 +18,13 @@ function full{T}(H::Hilbert{T})
     end
     Hf
 end
+
+function inv{T}(H::Hilbert{T})
+    invH=zeros(T,H.m,H.m)
+    for i=1:H.m, j=1:H.m
+        invH[i,j]=(-1)^(i+j)*(i+j-1)*binomial(H.m+i-1,H.m-j)*
+                  binomial(H.m+j-1,H.m-i)*binomial(i+j-2,i-1)^2
+    end
+    invH
+end
+

src/kahan.jl

@@ -0,0 +1,33 @@
+#### Kahan matrix
+
+export Kahan
+
+immutable Kahan{T<:Number,T1<:Number} <: AbstractMatrix{T}
+    m::Int # dimension
+    n::Int # dimension
+    theta::T # angle
+    pert::T1 # perturbation is pert*eps()
+end # immutable
+
+# Define its size
+size(A::Kahan, r::Int) = r==1 ? A.m : A.n
+size(A::Kahan) = A.m, A.n
+
+# Index into a Kahan
+function getindex(A::Kahan,i::Integer,j::Integer)
+    m=minimum(size(A))
+    t=tan(A.theta)
+    c=1.0/sqrt(1.0+t^2)
+    s=t*c
+    if i>m; return 0.0
+    elseif i>j; return 0.0
+    elseif i==j; return s^(i-1)+A.pert*eps()*(m-i+1)
+    else return -c*s^(i-1)
+    end
+end # getindex
+
+# Dense version of Kahan
+full(A::Kahan) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)]
+
+
+

src/riemann.jl

@@ -0,0 +1,26 @@
+#### Riemann Matrix 
+import Base: size
+export Riemann
+
+immutable Riemann{Int} <: AbstractMatrix{Int}
+    n::Int
+end # immutable
+
+# Define its size
+
+size(A::Riemann, dim::Integer) = A.n
+size(A::Riemann)= size(A,1), size(A,1)
+
+# Index into a Riemann
+function getindex(A::Riemann,i::Integer,j::Integer)
+#    return (i+1)%(j+1)==0 ? i : -1
+    if i<=j && (j+1)%(i+1)==0
+        return i
+    else
+        return -1
+    end
+end # getindex
+
+# Dense version of Riemann
+full(A::Riemann) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)]
+