Skip to content

Commit 27906b3

Browse files
Mahmood-SinanMahmood-Sinan
committed
extended lagtm implemented, but only for real alpha and beta values
1 parent 34723c7 commit 27906b3

File tree

5 files changed

+125
-2
lines changed

5 files changed

+125
-2
lines changed

‎src/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@ set(fppFiles
99
stdlib_codata_type.fypp
1010
stdlib_constants.fypp
1111
stdlib_error.fypp
12+
stdlib_extended_lapack_base.fypp
13+
stdlib_extended_lapack.fypp
1214
stdlib_hash_32bit.fypp
1315
stdlib_hash_32bit_fnv.fypp
1416
stdlib_hash_32bit_nm.fypp

‎src/stdlib_extended_lapack.fypp‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4+
#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
5+
6+
submodule(stdlib_extended_lapack_base) stdlib_extended_lapack
7+
implicit none
8+
contains
9+
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
10+
#:for k1,t1,s1 in KINDS_TYPES
11+
pure module subroutine stdlib${ii}$_glagtm_${s1}$(trans, n, nrhs, alpha, dl, d, du, x, ldx, beta, b, ldb)
12+
character, intent(in) :: trans
13+
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
14+
real(${k1}$), intent(in) :: alpha, beta
15+
${t1}$, intent(inout) :: b(ldb,*)
16+
${t1}$, intent(in) :: d(*), dl(*), du(*), x(ldx,*)
17+
18+
! Internal variables.
19+
integer(${ik}$) :: i, j
20+
${t1}$ :: temp
21+
if(n == 0) then
22+
return
23+
endif
24+
if(beta == 0.0_${k1}$) then
25+
b(1:n, 1:nrhs) = 0.0_${k1}$
26+
else
27+
b(1:n, 1:nrhs) = beta * b(1:n, 1:nrhs)
28+
end if
29+
30+
if(trans == 'N') then
31+
do j = 1, nrhs
32+
if(n == 1_${ik}$) then
33+
temp = d(1_${ik}$) * x(1_${ik}$, j)
34+
b(1_${ik}$, j) = b(1_${ik}$, j) + alpha * temp
35+
else
36+
temp = d(1_${ik}$) * x(1_${ik}$, j) + du(1_${ik}$) * x(2_${ik}$, j)
37+
b(1_${ik}$, j) = b(1_${ik}$, j) + alpha * temp
38+
do i = 2, n - 1
39+
temp = dl(i - 1) * x(i - 1, j) + d(i) * x(i, j) + du(i) * x(i + 1, j)
40+
b(i, j) = b (i, j) + alpha * temp
41+
end do
42+
temp = dl(n - 1) * x(n - 1, j) + d(n) * x(n, j)
43+
b(n, j) = b(n, j) + alpha * temp
44+
end if
45+
end do
46+
#:if t1.startswith('complex')
47+
else if(trans == 'C') then
48+
do j = 1, nrhs
49+
if(n == 1_${ik}$) then
50+
temp = conjg(d(1_${ik}$)) * x(1_${ik}$, j)
51+
b(1_${ik}$, j) = b(1_${ik}$, j) + alpha * temp
52+
else
53+
temp = conjg(d(1_${ik}$)) * x(1_${ik}$, j) + conjg(dl(1_${ik}$)) * x(2_${ik}$, j)
54+
b(1_${ik}$, j) = b(1_${ik}$, j) + alpha * temp
55+
do i = 2, n - 1
56+
temp = conjg(du(i - 1)) * x(i - 1, j) + conjg(d(i)) * x(i, j) + conjg(dl(i)) * x(i + 1, j)
57+
b(i, j) = b (i, j) + alpha * temp
58+
end do
59+
temp = conjg(du(n - 1)) * x(n - 1, j) + conjg(d(n)) * x(n, j)
60+
b(n, j) = b(n, j) + alpha * temp
61+
end if
62+
end do
63+
#:endif
64+
else
65+
do j = 1, nrhs
66+
if(n == 1_${ik}$) then
67+
temp = d(1_${ik}$) * x(1_${ik}$, j)
68+
b(1_${ik}$, j) = b(1_${ik}$, j) + alpha * temp
69+
else
70+
temp = d(1_${ik}$) * x(1_${ik}$, j) + dl(1_${ik}$) * x(2_${ik}$, j)
71+
b(1_${ik}$, j) = b(1_${ik}$, j) + alpha * temp
72+
do i = 2, n - 1
73+
temp = du(i - 1) * x(i - 1, j) + d(i) * x(i, j) + dl(i) * x(i + 1, j)
74+
b(i, j) = b (i, j) + alpha * temp
75+
end do
76+
temp = du(n - 1) * x(n - 1, j) + d(n) * x(n, j)
77+
b(n, j) = b(n, j) + alpha * temp
78+
end if
79+
end do
80+
end if
81+
end subroutine stdlib${ii}$_glagtm_${s1}$
82+
#:endfor
83+
#:endfor
84+
85+
end submodule

‎src/stdlib_extended_lapack_base.fypp‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4+
#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
5+
module stdlib_extended_lapack_base
6+
use stdlib_linalg_constants
7+
implicit none
8+
9+
interface glagtm
10+
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
11+
#:for k1,t1,s1 in KINDS_TYPES
12+
pure module subroutine stdlib${ii}$_glagtm_${s1}$(trans, n, nrhs, alpha, dl, d, du, x, ldx, beta, b, ldb)
13+
character, intent(in) :: trans
14+
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
15+
real(${k1}$), intent(in) :: alpha, beta
16+
${t1}$, intent(inout) :: b(ldb,*)
17+
${t1}$, intent(in) :: d(*), dl(*), du(*), x(ldx,*)
18+
end subroutine stdlib${ii}$_glagtm_${s1}$
19+
#:endfor
20+
#:endfor
21+
end interface
22+
end module

‎src/stdlib_specialmatrices.fypp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ module stdlib_specialmatrices
1212
use stdlib_constants
1313
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
1414
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
15+
use stdlib_extended_lapack_base
1516
implicit none
1617
private
1718
public :: tridiagonal

‎src/stdlib_specialmatrices_tridiagonal.fypp‎

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -163,6 +163,8 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
163163
real(${k1}$) :: alpha_, beta_
164164
integer(ilp) :: n, nrhs, ldx, ldy
165165
character(1) :: op_
166+
logical :: alpha_special, beta_special
167+
166168
#:if rank == 1
167169
${t1}$, pointer :: xmat(:, :), ymat(:, :)
168170
#:endif
@@ -172,16 +174,27 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
172174
beta_ = 0.0_${k1}$ ; if (present(beta)) beta_ = beta
173175
op_ = "N" ; if (present(op)) op_ = op
174176

177+
alpha_special = (alpha_ == 1.0_${k1}$ .or. alpha_ == 0.0_${k1}$ .or. alpha_ == -1.0_${k1}$)
178+
beta_special = (beta_ == 1.0_${k1}$ .or. beta_ == 0.0_${k1}$ .or. beta_ == -1.0_${k1}$)
179+
175180
! Prepare Lapack arguments.
176181
n = A%n ; ldx = n ; ldy = n ;
177182
nrhs = #{if rank==1}# 1 #{else}# size(x, dim=2, kind=ilp) #{endif}#
178183

179184
#:if rank == 1
180185
! Pointer trick.
181186
xmat(1:n, 1:nrhs) => x ; ymat(1:n, 1:nrhs) => y
182-
call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy)
187+
if(alpha_special .and. beta_special) then
188+
call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy)
189+
else
190+
call glagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy)
191+
end if
183192
#:else
184-
call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy)
193+
if(alpha_special .and. beta_special) then
194+
call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy)
195+
else
196+
call glagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy)
197+
end if
185198
#:endif
186199
end subroutine
187200
#:endfor

0 commit comments

Comments
 (0)