Skip to content

Commit d64ebef

Browse files
loiseaujcloiseaujc
committed
SymTridiagonal is extended from Tridiagonal instead of being its own type.
1 parent 9a37e13 commit d64ebef

File tree

3 files changed

+106
-243
lines changed

3 files changed

+106
-243
lines changed

‎src/stdlib_specialmatrices.fypp

Lines changed: 19 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -37,12 +37,12 @@ module stdlib_specialmatrices
3737

3838
!--> Symmetric Tridiagonal matrices
3939
#:for k1, t1, s1 in (KINDS_TYPES)
40-
type, public :: symtridiagonal_${s1}$_type
40+
type, extends(tridiagonal_${s1}$_type), public :: symtridiagonal_${s1}$_type
4141
!! Base type to define a `symtridiagonal` matrix.
4242
private
43-
${t1}$, allocatable :: dv(:), ev(:)
44-
integer(ilp) :: n
43+
#:if t1.startswith('real')
4544
logical(lk) :: is_posdef
45+
#:endif
4646
end type
4747
#:endfor
4848

@@ -238,16 +238,7 @@ module stdlib_specialmatrices
238238
#:for k1, t1, s1 in (KINDS_TYPES)
239239
#:for rank in RANKS
240240
module subroutine spmv_tridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op)
241-
type(tridiagonal_${s1}$_type), intent(in) :: A
242-
${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$
243-
${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$
244-
real(${k1}$), intent(in), optional :: alpha
245-
real(${k1}$), intent(in), optional :: beta
246-
character(1), intent(in), optional :: op
247-
end subroutine
248-
249-
module subroutine spmv_symtridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op)
250-
type(symtridiagonal_${s1}$_type), intent(in) :: A
241+
class(tridiagonal_${s1}$_type), intent(in) :: A
251242
${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$
252243
${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$
253244
real(${k1}$), intent(in), optional :: alpha
@@ -271,19 +262,11 @@ module stdlib_specialmatrices
271262
#:for k1, t1, s1 in (KINDS_TYPES)
272263
pure module function tridiagonal_to_dense_${s1}$(A) result(B)
273264
!! Convert a `tridiagonal` matrix to its dense representation.
274-
type(tridiagonal_${s1}$_type), intent(in) :: A
265+
class(tridiagonal_${s1}$_type), intent(in) :: A
275266
!! Input Tridiagonal matrix.
276267
${t1}$, allocatable :: B(:, :)
277268
!! Corresponding dense matrix.
278269
end function
279-
280-
pure module function symtridiagonal_to_dense_${s1}$(A) result(B)
281-
!! Convert a `symtridiagonal` matrix to its dense representation.
282-
type(symtridiagonal_${s1}$_type), intent(in) :: A
283-
!! Input SymTridiagonal matrix.
284-
${t1}$, allocatable :: B(:, :)
285-
!! Corresponding dense matrix.
286-
end function
287270
#:endfor
288271
end interface
289272

@@ -293,15 +276,9 @@ module stdlib_specialmatrices
293276
!! [Specifications](../page/specs/stdlib_specialmatrices.html#transpose)
294277
#:for k1, t1, s1 in (KINDS_TYPES)
295278
pure module function transpose_tridiagonal_${s1}$(A) result(B)
296-
type(tridiagonal_${s1}$_type), intent(in) :: A
279+
class(tridiagonal_${s1}$_type), intent(in) :: A
297280
!! Input matrix.
298-
type(tridiagonal_${s1}$_type) :: B
299-
end function
300-
301-
pure module function transpose_symtridiagonal_${s1}$(A) result(B)
302-
type(symtridiagonal_${s1}$_type), intent(in) :: A
303-
!! Input matrix.
304-
type(symtridiagonal_${s1}$_type) :: B
281+
class(tridiagonal_${s1}$_type), allocatable :: B
305282
end function
306283
#:endfor
307284
end interface
@@ -313,15 +290,9 @@ module stdlib_specialmatrices
313290
!! [Specifications](../page/specs/stdlib_specialmatrices.html#hermitian)
314291
#:for k1, t1, s1 in (KINDS_TYPES)
315292
pure module function hermitian_tridiagonal_${s1}$(A) result(B)
316-
type(tridiagonal_${s1}$_type), intent(in) :: A
317-
!! Input matrix.
318-
type(tridiagonal_${s1}$_type) :: B
319-
end function
320-
321-
pure module function hermitian_symtridiagonal_${s1}$(A) result(B)
322-
type(symtridiagonal_${s1}$_type), intent(in) :: A
293+
class(tridiagonal_${s1}$_type), intent(in) :: A
323294
!! Input matrix.
324-
type(symtridiagonal_${s1}$_type) :: B
295+
class(tridiagonal_${s1}$_type), allocatable :: B
325296
end function
326297
#:endfor
327298
end interface
@@ -339,24 +310,13 @@ module stdlib_specialmatrices
339310
#:for k1, t1, s1 in (KINDS_TYPES)
340311
pure module function scalar_multiplication_tridiagonal_${s1}$(alpha, A) result(B)
341312
${t1}$, intent(in) :: alpha
342-
type(tridiagonal_${s1}$_type), intent(in) :: A
343-
type(tridiagonal_${s1}$_type) :: B
313+
class(tridiagonal_${s1}$_type), intent(in) :: A
314+
class(tridiagonal_${s1}$_type), allocatable :: B
344315
end function
345316
pure module function scalar_multiplication_bis_tridiagonal_${s1}$(A, alpha) result(B)
346-
type(tridiagonal_${s1}$_type), intent(in) :: A
347-
${t1}$, intent(in) :: alpha
348-
type(tridiagonal_${s1}$_type) :: B
349-
end function
350-
351-
pure module function scalar_multiplication_symtridiagonal_${s1}$(alpha, A) result(B)
317+
class(tridiagonal_${s1}$_type), intent(in) :: A
352318
${t1}$, intent(in) :: alpha
353-
type(symtridiagonal_${s1}$_type), intent(in) :: A
354-
type(symtridiagonal_${s1}$_type) :: B
355-
end function
356-
pure module function scalar_multiplication_bis_symtridiagonal_${s1}$(A, alpha) result(B)
357-
type(symtridiagonal_${s1}$_type), intent(in) :: A
358-
${t1}$, intent(in) :: alpha
359-
type(symtridiagonal_${s1}$_type) :: B
319+
class(tridiagonal_${s1}$_type), allocatable :: B
360320
end function
361321
#:endfor
362322
end interface
@@ -367,15 +327,9 @@ module stdlib_specialmatrices
367327
!! [Specifications](../page/specs/stdlib_specialmatrices.html#operators)
368328
#:for k1, t1, s1 in (KINDS_TYPES)
369329
pure module function matrix_add_tridiagonal_${s1}$(A, B) result(C)
370-
type(tridiagonal_${s1}$_type), intent(in) :: A
371-
type(tridiagonal_${s1}$_type), intent(in) :: B
372-
type(tridiagonal_${s1}$_type) :: C
373-
end function
374-
375-
pure module function matrix_add_symtridiagonal_${s1}$(A, B) result(C)
376-
type(symtridiagonal_${s1}$_type), intent(in) :: A
377-
type(symtridiagonal_${s1}$_type), intent(in) :: B
378-
type(symtridiagonal_${s1}$_type) :: C
330+
class(tridiagonal_${s1}$_type), intent(in) :: A
331+
class(tridiagonal_${s1}$_type), intent(in) :: B
332+
class(tridiagonal_${s1}$_type), allocatable :: C
379333
end function
380334
#:endfor
381335
end interface
@@ -386,15 +340,9 @@ module stdlib_specialmatrices
386340
!! [Specifications](../page/specs/stdlib_specialmatrices.html#operators)
387341
#:for k1, t1, s1 in (KINDS_TYPES)
388342
pure module function matrix_sub_tridiagonal_${s1}$(A, B) result(C)
389-
type(tridiagonal_${s1}$_type), intent(in) :: A
390-
type(tridiagonal_${s1}$_type), intent(in) :: B
391-
type(tridiagonal_${s1}$_type) :: C
392-
end function
393-
394-
pure module function matrix_sub_symtridiagonal_${s1}$(A, B) result(C)
395-
type(symtridiagonal_${s1}$_type), intent(in) :: A
396-
type(symtridiagonal_${s1}$_type), intent(in) :: B
397-
type(symtridiagonal_${s1}$_type) :: C
343+
class(tridiagonal_${s1}$_type), intent(in) :: A
344+
class(tridiagonal_${s1}$_type), intent(in) :: B
345+
class(tridiagonal_${s1}$_type), allocatable :: C
398346
end function
399347
#:endfor
400348
end interface

‎src/stdlib_specialmatrices_symtridiagonal.fypp

Lines changed: 6 additions & 142 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ submodule (stdlib_specialmatrices) symtridiagonal_matrices
4242
! Description of the matrix.
4343
A%n = n
4444
! Matrix elements.
45-
A%dv = dv ; A%ev = ev
45+
A%dl = ev ; A%dv = dv ; A%du = ev
4646
end function
4747

4848
pure module function initialize_constant_symtridiagonal_pure_${s1}$(dv, ev, n) result(A)
@@ -65,8 +65,9 @@ submodule (stdlib_specialmatrices) symtridiagonal_matrices
6565
call linalg_error_handling(err0)
6666
endif
6767
! Matrix elements.
68+
A%dl = [(ev, i = 1, n-1)]
6869
A%dv = [(dv, i = 1, n-1)]
69-
A%ev = [(ev, i = 1, n-1)]
70+
A%du = [(ev, i = 1, n-1)]
7071
end function
7172

7273
module function initialize_symtridiagonal_impure_${s1}$(dv, ev, err) result(A)
@@ -97,7 +98,7 @@ submodule (stdlib_specialmatrices) symtridiagonal_matrices
9798
! Description of the matrix.
9899
A%n = n
99100
! Matrix elements.
100-
A%dv = dv ; A%ev = ev
101+
A%dl = ev ; A%dv = dv ; A%du = ev
101102
end function
102103

103104
module function initialize_constant_symtridiagonal_impure_${s1}$(dv, ev, n, err) result(A)
@@ -122,146 +123,9 @@ submodule (stdlib_specialmatrices) symtridiagonal_matrices
122123
call linalg_error_handling(err0, err)
123124
endif
124125
! Matrix elements.
126+
A%dl = [(ev, i = 1, n)]
125127
A%dv = [(dv, i = 1, n-1)]
126-
A%ev = [(ev, i = 1, n)]
127-
end function
128-
#:endfor
129-
130-
!-----------------------------------------
131-
!----- -----
132-
!----- MATRIX-VECTOR PRODUCT -----
133-
!----- -----
134-
!-----------------------------------------
135-
136-
!! spmv_tridiag
137-
#:for k1, t1, s1 in (KINDS_TYPES)
138-
#:for rank in RANKS
139-
module subroutine spmv_symtridiag_${rank}$d_${s1}$(A, x, y, alpha, beta, op)
140-
type(symtridiagonal_${s1}$_type), intent(in) :: A
141-
${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$
142-
${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$
143-
real(${k1}$), intent(in), optional :: alpha
144-
real(${k1}$), intent(in), optional :: beta
145-
character(1), intent(in), optional :: op
146-
147-
! Internal variables.
148-
real(${k1}$) :: alpha_, beta_
149-
integer(ilp) :: n, nrhs, ldx, ldy
150-
character(1) :: op_
151-
#:if rank == 1
152-
${t1}$, pointer :: xmat(:, :), ymat(:, :)
153-
#:endif
154-
155-
! Deal with optional arguments.
156-
alpha_ = 1.0_${k1}$ ; if (present(alpha)) alpha_ = alpha
157-
beta_ = 0.0_${k1}$ ; if (present(beta)) beta_ = beta
158-
op_ = "N" ; if (present(op)) op_ = op
159-
if (op_ == "H") op_ = "C"
160-
161-
! Prepare Lapack arguments.
162-
n = A%n ; ldx = n ; ldy = n ; y = 0.0_${k1}$
163-
nrhs = #{if rank==1}# 1 #{else}# size(x, dim=2, kind=ilp) #{endif}#
164-
165-
#:if rank == 1
166-
! Pointer trick.
167-
xmat(1:n, 1:nrhs) => x ; ymat(1:n, 1:nrhs) => y
168-
call lagtm(op_, n, nrhs, alpha_, A%ev, A%dv, A%ev, xmat, ldx, beta_, ymat, ldy)
169-
#:else
170-
call lagtm(op_, n, nrhs, alpha_, A%ev, A%dv, A%ev, x, ldx, beta_, y, ldy)
171-
#:endif
172-
end subroutine
173-
#:endfor
174-
#:endfor
175-
176-
!-------------------------------------
177-
!----- -----
178-
!----- UTILITY FUNCTIONS -----
179-
!----- -----
180-
!-------------------------------------
181-
182-
#:for k1, t1, s1 in (KINDS_TYPES)
183-
pure module function symtridiagonal_to_dense_${s1}$(A) result(B)
184-
!! Convert a `symtridiagonal` matrix to its dense representation.
185-
type(symtridiagonal_${s1}$_type), intent(in) :: A
186-
!! Input symtridiagonal matrix.
187-
${t1}$, allocatable :: B(:, :)
188-
!! Corresponding dense matrix.
189-
190-
! Internal variables.
191-
integer(ilp) :: i
192-
193-
associate (n => A%n)
194-
#:if t1.startswith('complex')
195-
allocate(B(n, n), source=zero_c${k1}$)
196-
#:else
197-
allocate(B(n, n), source=zero_${k1}$)
198-
#:endif
199-
B(1, 1) = A%dv(1) ; B(1, 2) = A%ev(1)
200-
do concurrent (i=2:n-1)
201-
B(i, i-1) = A%ev(i-1)
202-
B(i, i) = A%dv(i)
203-
B(i, i+1) = A%ev(i)
204-
enddo
205-
B(n, n-1) = A%ev(n-1) ; B(n, n) = A%dv(n)
206-
end associate
207-
end function
208-
#:endfor
209-
210-
#:for k1, t1, s1 in (KINDS_TYPES)
211-
pure module function transpose_symtridiagonal_${s1}$(A) result(B)
212-
type(symtridiagonal_${s1}$_type), intent(in) :: A
213-
!! Input matrix.
214-
type(symtridiagonal_${s1}$_type) :: B
215-
B = A
216-
end function
217-
#:endfor
218-
219-
#:for k1, t1, s1 in (KINDS_TYPES)
220-
pure module function hermitian_symtridiagonal_${s1}$(A) result(B)
221-
type(symtridiagonal_${s1}$_type), intent(in) :: A
222-
!! Input matrix.
223-
type(symtridiagonal_${s1}$_type) :: B
224-
#:if t1.startswith("complex")
225-
B = symtridiagonal(A%dv, conjg(A%ev))
226-
#:else
227-
B = symtridiagonal(A%dv, A%ev)
228-
#:endif
229-
end function
230-
#:endfor
231-
232-
#:for k1, t1, s1 in (KINDS_TYPES)
233-
pure module function scalar_multiplication_symtridiagonal_${s1}$(alpha, A) result(B)
234-
${t1}$, intent(in) :: alpha
235-
type(symtridiagonal_${s1}$_type), intent(in) :: A
236-
type(symtridiagonal_${s1}$_type) :: B
237-
B = symtridiagonal(A%dv, A%ev)
238-
B%dv = alpha*B%dv; B%ev = alpha*B%ev
239-
end function
240-
241-
pure module function scalar_multiplication_bis_symtridiagonal_${s1}$(A, alpha) result(B)
242-
type(symtridiagonal_${s1}$_type), intent(in) :: A
243-
${t1}$, intent(in) :: alpha
244-
type(symtridiagonal_${s1}$_type) :: B
245-
B = symtridiagonal(A%dv, A%ev)
246-
B%dv = alpha*B%dv; B%ev = alpha*B%ev
247-
end function
248-
#:endfor
249-
250-
#:for k1, t1, s1 in (KINDS_TYPES)
251-
pure module function matrix_add_symtridiagonal_${s1}$(A, B) result(C)
252-
type(symtridiagonal_${s1}$_type), intent(in) :: A
253-
type(symtridiagonal_${s1}$_type), intent(in) :: B
254-
type(symtridiagonal_${s1}$_type) :: C
255-
C = symtridiagonal(A%dv, A%ev)
256-
C%dv = C%dv + B%dv; C%ev = C%ev + B%ev
257-
end function
258-
259-
pure module function matrix_sub_symtridiagonal_${s1}$(A, B) result(C)
260-
type(symtridiagonal_${s1}$_type), intent(in) :: A
261-
type(symtridiagonal_${s1}$_type), intent(in) :: B
262-
type(symtridiagonal_${s1}$_type) :: C
263-
C = symtridiagonal(A%dv, A%ev)
264-
C%dv = C%dv - B%dv; C%ev = C%ev - B%ev
128+
A%du = [(ev, i = 1, n)]
265129
end function
266130
#:endfor
267131

0 commit comments

Comments
 (0)