Defining Sparse Matrix from Scratch using vectors

General Usage
,
1

Hi

I want to create a custom tridiagonal matrix type whose inputs will be 3 different vectors (c, d and e).
Vecs
These 3 vectors will make the matrix in a format as follows:
TD_mat

I tried writing a struct

struct tridiag <: AbstractMatrix{<:AbstractFloat}
    c::Array{<:AbstractFloat}
    d::Array{<:AbstractFloat}
    e::Array{<:AbstractFloat}
end

but I’m getting an error

invalid subtyping in definition of tridiag

Stacktrace:
 [1] top-level scope at In[14]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

I’m also not sure how I can define a custom type for tridiagonal (sprase) matrix effectively like this. Any other ideas are welcome and I want to writing this type from scratch and not use something like Sparce matrix.

2

Will Tridiagonal do what you want? You may have to do using LinearAlgebra to get it, but it’s built in.

julia> a=[1,1,1]; b=[2,2,2,2]; c=[3,3,3];

julia> Tridiagonal(a,b,c)
4×4 Tridiagonal{Int64,Array{Int64,1}}:
 2  3  ⋅  ⋅
 1  2  3  ⋅
 ⋅  1  2  3
 ⋅  ⋅  1  2
2 Likes
3

I’m just starting out with Julia and learning by writing basic numerical methods from scratch without any libraries. I want to do something like this but want to apply from scratch. I think I can create something similar using struct. Maybe I can read the source code for this linear algebra implementation of tridiagonal. Can you suggest me something here?

4

Here are the docs Linear Algebra · The Julia Language and there’s a link to the source code in the lower right corner (you might need to select or hover for it to show up).

1 Like
6

You want the syntax:

struct Tridiag{T<:Number} <: AbstractMatrix{T}
    c::Vector{T}
    d::Vector{T}
    e::Vector{T}
end

Note that this defines a whole set of types Tridiag{Float64}, Tridiag{Int}, and so on, where for each of which the fields are concrete types. (Note also that Array{T} is abstract because you didn’t specify the dimensionality — you want Array{T,1} or equivalently Vector{T}.)

5 Likes
7

You can also see the Alan Edelman’s video on Matrix Structures for inspiration.

I have written (to learn and for fun) the structure as suggested by Steven and some two other functions, and already matrix multiplications work. This is very nice :slight_smile:

julia> A = Tridiag([1,1],[2,2,2],[3,3])
3×3 Tridiag{Int64}:
 2  3  0
 1  2  3
 0  1  2

julia> A*A
3×3 Array{Int64,2}:
 7  12   9
 4  10  12
 1   4   7

If you want to see what I did so far, here is the code:

Code 
julia> struct Tridiag{T<:Number} <: AbstractMatrix{T} 
         a :: Vector{T}
         b :: Vector{T}
         c :: Vector{T}
       end

julia> function Base.getindex(M::Tridiag{T}, i::Int, j::Int) where T
         if i == j
           M.b[i]
         elseif i == j + 1
           M.a[j]
         elseif i == j - 1
           M.c[i]
         else
           zero(T)
         end
       end

julia> Base.size(M::Tridiag) = (length(M.b),length(M.b))

julia> A = Tridiag([1,1],[2,2,2],[3,3])
3×3 Tridiag{Int64}:
 2  3  0
 1  2  3
 0  1  2

julia> A*A
3×3 Array{Int64,2}:
 7  12   9
 4  10  12
 1   4   7
4 Likes
8

@lmiq, thanks for sharing that nice video link and your code too. If one has missed the previous lessons in the series, one might not know that some of the examples shown require using several packages (ex: LinearAlgebra, SparseArrays, Statistics).

1 Like