Skip to content

Commit fc61941

Browse files
loiseaujcloiseaujc
committed
Implementation of basic operations for Tridiagonal matrices is done.
1 parent 9ae0554 commit fc61941

File tree

3 files changed

+287
-78
lines changed

3 files changed

+287
-78
lines changed

‎src/CMakeLists.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,8 +52,8 @@ set(fppFiles
5252
stdlib_sparse_conversion.fypp
5353
stdlib_sparse_kinds.fypp
5454
stdlib_sparse_spmv.fypp
55-
stdlib_specialmatrices_tridiagonal.fypp
5655
stdlib_specialmatrices.fypp
56+
stdlib_specialmatrices_tridiagonal.fypp
5757
stdlib_specialfunctions_gamma.fypp
5858
stdlib_stats.fypp
5959
stdlib_stats_corr.fypp

‎src/stdlib_specialmatrices.fypp

Lines changed: 143 additions & 77 deletions
Original file line numberDiff line numberDiff line change
@@ -8,18 +8,34 @@ module stdlib_specialmatrices
88
use stdlib_linalg_constants
99
implicit none
1010
private
11-
public :: Tridiagonal, spmv, dense
12-
13-
!! version: experimental
11+
public :: Tridiagonal
12+
public :: spmv
13+
public :: dense, transpose, hermitian
14+
public :: operator(*), operator(+), operator(-)
15+
16+
!--------------------------------------
17+
!----- ------
18+
!----- TYPE DEFINITIONS ------
19+
!----- ------
20+
!--------------------------------------
21+
22+
!! Version: experimental
1423
!!
1524
!! Tridiagonal matrix
1625
#:for k1, t1, s1 in (KINDS_TYPES)
1726
type, public :: Tridiagonal_${s1}$_type
27+
private
1828
${t1}$, allocatable :: dl(:), dv(:), du(:)
19-
integer(ilp) :: n, m
29+
integer(ilp) :: n
2030
end type
2131
#:endfor
2232

33+
!--------------------------------
34+
!----- -----
35+
!----- CONSTRUCTORS -----
36+
!----- -----
37+
!--------------------------------
38+
2339
interface Tridiagonal
2440
!! This interface provides different methods to construct a
2541
!! `Tridiagonal` matrix. Only the non-zero elements of \( A \) are
@@ -61,11 +77,37 @@ module stdlib_specialmatrices
6177
!! A = Tridiagonal(a, b, c, n)
6278
!! ```
6379
#:for k1, t1, s1 in (KINDS_TYPES)
64-
module procedure initialize_tridiagonal_${s1}$
65-
module procedure initialize_constant_tridiagonal_${s1}$
80+
pure module function initialize_tridiagonal_${s1}$(dl, dv, du) result(A)
81+
!! Construct a `Tridiagonal` matrix from the rank-1 arrays
82+
!! `dl`, `dv` and `du`.
83+
${t1}$, intent(in) :: dl(:), dv(:), du(:)
84+
!! Tridiagonal matrix elements.
85+
type(Tridiagonal_${s1}$_type) :: A
86+
!! Corresponding Tridiagonal matrix.
87+
end function
88+
89+
pure module function initialize_constant_tridiagonal_${s1}$(dl, dv, du, n) result(A)
90+
!! Construct a `Tridiagonal` matrix with constant elements.
91+
${t1}$, intent(in) :: dl, dv, du
92+
!! Tridiagonal matrix elements.
93+
integer(ilp), intent(in) :: n
94+
!! Matrix dimension.
95+
type(Tridiagonal_${s1}$_type) :: A
96+
!! Corresponding Tridiagonal matrix.
97+
end function
6698
#:endfor
6799
end interface
68100

101+
!----------------------------------
102+
!----- -----
103+
!----- LINEAR ALGEBRA -----
104+
!----- -----
105+
!----------------------------------
106+
107+
!! Version: experimental
108+
!!
109+
!! Apply the matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$
110+
!! [Specifications](../page/specs/stdlib_specialmatrices.html#spmv)
69111
interface spmv
70112
#:for k1, t1, s1 in (KINDS_TYPES)
71113
#:for rank in RANKS
@@ -83,6 +125,16 @@ module stdlib_specialmatrices
83125
#:endfor
84126
end interface
85127

128+
!-------------------------------------
129+
!----- -----
130+
!----- UTILITY FUNCTIONS -----
131+
!----- -----
132+
!-------------------------------------
133+
134+
!! Version: experimental
135+
!!
136+
!! Convert a matrix of type `Tridiagonal` to its dense representation.
137+
!! [Specifications](../page/specs/stdlib_specialmatrices.html#dense)
86138
interface dense
87139
!! This interface provides methods to convert a `Tridiagonal` matrix
88140
!! to a regular rank-2 array.
@@ -93,81 +145,95 @@ module stdlib_specialmatrices
93145
!! B = dense(A)
94146
!! ```
95147
#:for k1, t1, s1 in (KINDS_TYPES)
96-
module procedure tridiagonal_to_dense_${s1}$
148+
pure module function tridiagonal_to_dense_${s1}$(A) result(B)
149+
!! Convert a `Tridiagonal` matrix to its dense representation.
150+
type(Tridiagonal_${s1}$_type), intent(in) :: A
151+
!! Input Tridiagonal matrix.
152+
${t1}$, allocatable :: B(:, :)
153+
!! Corresponding dense matrix.
154+
end function
97155
#:endfor
98156
end interface
99157

100-
contains
158+
!! Version: experimental
159+
!!
160+
!! Returns the transpose of a `Tridiagonal` matrix.
161+
!! [Specifications](../page/specs/stdlib_specialmatrices.html#transpose)
162+
interface transpose
163+
#:for k1, t1, s1 in (KINDS_TYPES)
164+
pure module function transpose_tridiagonal_${s1}$(A) result(B)
165+
type(Tridiagonal_${s1}$_type), intent(in) :: A
166+
!! Input matrix.
167+
type(Tridiagonal_${s1}$_type) :: B
168+
end function
169+
#:endfor
170+
end interface
101171

102-
!------------------------------------------------------
103-
!----- -----
104-
!----- Tridiagonal matrix implementations -----
105-
!----- -----
106-
!------------------------------------------------------
172+
!! Version: experimental
173+
!!
174+
!! Returns the Hermitian of a `Tridiagonal` matrix.
175+
!! [Specifications](../page/specs/stdlib_specialmatrices.html#hermitian)
176+
interface hermitian
177+
#:for k1, t1, s1 in (KINDS_TYPES)
178+
pure module function hermitian_tridiagonal_${s1}$(A) result(B)
179+
type(Tridiagonal_${s1}$_type), intent(in) :: A
180+
!! Input matrix.
181+
type(Tridiagonal_${s1}$_type) :: B
182+
end function
183+
#:endfor
184+
end interface
107185

108-
#:for k1, t1, s1 in (KINDS_TYPES)
109-
pure module function initialize_tridiagonal_${s1}$(dl, dv, du) result(A)
110-
!! Construct a `Tridiagonal` matrix from the rank-1 arrays
111-
!! `dl`, `dv` and `du`.
112-
${t1}$, intent(in) :: dl(:), dv(:), du(:)
113-
!! Tridiagonal matrix elements.
114-
type(Tridiagonal_${s1}$_type) :: A
115-
!! Corresponding Tridiagonal matrix.
116-
117-
! Internal variables.
118-
integer(ilp) :: n
186+
!----------------------------------------
187+
!----- -----
188+
!----- ARITHMETIC OPERATORS -----
189+
!----- -----
190+
!----------------------------------------
191+
192+
!! Version: experimental
193+
!!
194+
!! Overloads the scalar multiplication `*` for `Tridiagonal` matrices.
195+
!! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(*))
196+
interface operator(*)
197+
#:for k1, t1, s1 in (KINDS_TYPES)
198+
pure module function scalar_multiplication_tridiagonal_${s1}$(alpha, A) result(B)
199+
${t1}$, intent(in) :: alpha
200+
type(Tridiagonal_${s1}$_type), intent(in) :: A
201+
type(Tridiagonal_${s1}$_type) :: B
202+
end function
203+
pure module function scalar_multiplication_bis_tridiagonal_${s1}$(A, alpha) result(B)
204+
type(Tridiagonal_${s1}$_type), intent(in) :: A
205+
${t1}$, intent(in) :: alpha
206+
type(Tridiagonal_${s1}$_type) :: B
207+
end function
208+
#:endfor
209+
end interface
210+
211+
!! Version: experimental
212+
!!
213+
!! Overloads the addition `+` for `Tridiagonal` matrices.
214+
!! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(+))
215+
interface operator(+)
216+
#:for k1, t1, s1 in (KINDS_TYPES)
217+
pure module function matrix_add_tridiagonal_${s1}$(A, B) result(C)
218+
type(Tridiagonal_${s1}$_type), intent(in) :: A
219+
type(Tridiagonal_${s1}$_type), intent(in) :: B
220+
type(Tridiagonal_${s1}$_type) :: C
221+
end function
222+
#:endfor
223+
end interface
224+
225+
!! Version: experimental
226+
!!
227+
!! Overloads the subtraction `-` for `Tridiagonal` matrices.
228+
!! [Specifications](../page/specs/stdlib_specialmatrices.html#operator(-))
229+
interface operator(-)
230+
#:for k1, t1, s1 in (KINDS_TYPES)
231+
pure module function matrix_sub_tridiagonal_${s1}$(A, B) result(C)
232+
type(Tridiagonal_${s1}$_type), intent(in) :: A
233+
type(Tridiagonal_${s1}$_type), intent(in) :: B
234+
type(Tridiagonal_${s1}$_type) :: C
235+
end function
236+
#:endfor
237+
end interface
119238

120-
! Sanity check.
121-
n = size(dv)
122-
if (size(dl) /= n-1) error stop "Vector dl does not have the correct length."
123-
if (size(du) /= n-1) error stop "Vector du does not have the correct length."
124-
125-
! Description of the matrix.
126-
A%n = n ; A%m = n
127-
! Matrix elements.
128-
A%dl = dl ; A%dv = dv ; A%du = du
129-
end function
130-
131-
pure module function initialize_constant_tridiagonal_${s1}$(dl, dv, du, n) result(A)
132-
!! Construct a `Tridiagonal` matrix with constant elements.
133-
${t1}$, intent(in) :: dl, dv, du
134-
!! Tridiagonal matrix elements.
135-
integer(ilp), intent(in) :: n
136-
!! Matrix dimension.
137-
type(Tridiagonal_${s1}$_type) :: A
138-
!! Corresponding Tridiagonal matrix.
139-
140-
! Internal variables.
141-
integer(ilp) :: i
142-
143-
! Description of the matrix.
144-
A%n = n ; A%m = n
145-
! Matrix elements.
146-
A%dl = [(dl, i = 1, n-1)]
147-
A%dv = [(dv, i = 1, n)]
148-
A%du = [(du, i = 1, n-1)]
149-
end function
150-
151-
pure module function tridiagonal_to_dense_${s1}$(A) result(B)
152-
!! Convert a `Tridiagonal` matrix to its dense representation.
153-
type(Tridiagonal_${s1}$_type), intent(in) :: A
154-
!! Input Tridiagonal matrix.
155-
${t1}$, allocatable :: B(:, :)
156-
!! Corresponding dense matrix.
157-
158-
! Internal variables.
159-
integer(ilp) :: i
160-
161-
associate (n => A%n)
162-
allocate(B(n, n)) ; B = 0.0_${k1}$
163-
B(1, 1) = A%dv(1) ; B(1, 2) = A%du(1)
164-
do concurrent (i=2:n-1)
165-
B(i, i-1) = A%dl(i-1)
166-
B(i, i) = A%dv(i)
167-
B(i, i+1) = A%du(i)
168-
enddo
169-
B(n, n-1) = A%dl(n-1) ; B(n, n) = A%dv(n)
170-
end associate
171-
end function
172-
#:endfor
173239
end module stdlib_specialmatrices

0 commit comments

Comments
 (0)