alpha and beta now supports complex values also

Author: Mahmood-Sinan

example/specialmatrices/CMakeLists.txt

@@ -1,2 +1,3 @@
 ADD_EXAMPLE(specialmatrices_dp_spmv)
+ADD_EXAMPLE(specialmatrices_cdp_spmv)
 ADD_EXAMPLE(tridiagonal_dp_type)

example/specialmatrices/example_specialmatrices_cdp_spmv.f90

@@ -0,0 +1,30 @@
+program example_tridiagonal_matrix_cdp
+    use stdlib_linalg_constants, only: dp
+    use stdlib_specialmatrices, only: tridiagonal_cdp_type, tridiagonal, dense, spmv
+    implicit none
+
+    integer, parameter :: n = 5
+    type(tridiagonal_cdp_type) :: A
+    complex(dp) :: dl(n-1), dv(n), du(n-1)
+    complex(dp) :: x(n), y(n), y_dense(n)
+    integer :: i
+    complex(dp) :: alpha, beta
+
+    dl = [(cmplx(i,0.0_dp, dp), i=1, n - 1)]
+    dv = [(cmplx(2*i,0.0_dp, dp), i=1, n     )]
+    du = [(cmplx(3*i,0.0_dp, dp), i=1, n - 1)]
+
+    A = tridiagonal(dl, dv, du)
+
+    x = (1.0_dp, 0.0_dp)
+    y = (0.0_dp, 0.0_dp)
+    y_dense = (0.0_dp, 0.0_dp)
+    alpha = cmplx(2.0_dp, 3.0_dp)
+    beta = cmplx(-1.0_dp, 5.0_dp)
+
+    y_dense = alpha * matmul(dense(A), x) + beta * y
+    call spmv(A, x, y, alpha, beta)
+
+    print *, 'dense       :', y_dense
+    print *, 'Tridiagonal :', y
+end program example_tridiagonal_matrix_cdp

src/stdlib_extended_lapack.fypp

@@ -11,7 +11,7 @@ contains
     pure module subroutine stdlib${ii}$_glagtm_${s1}$(trans, n, nrhs, alpha, dl, d, du, x, ldx, beta, b, ldb)
         character, intent(in) :: trans
         integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
-        real(${k1}$), intent(in) :: alpha, beta
+        ${t1}$, intent(in) :: alpha, beta
         ${t1}$, intent(inout) :: b(ldb,*)
         ${t1}$, intent(in) :: d(*), dl(*), du(*), x(ldx,*)
 

src/stdlib_extended_lapack_base.fypp

@@ -12,7 +12,7 @@ module stdlib_extended_lapack_base
         pure module subroutine stdlib${ii}$_glagtm_${s1}$(trans, n, nrhs, alpha, dl, d, du, x, ldx, beta, b, ldb)
             character, intent(in) :: trans
             integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
-            real(${k1}$), intent(in) :: alpha, beta
+            ${t1}$, intent(in) :: alpha, beta
             ${t1}$, intent(inout) :: b(ldb,*)
             ${t1}$, intent(in) :: d(*), dl(*), du(*), x(ldx,*)
         end subroutine stdlib${ii}$_glagtm_${s1}$

src/stdlib_specialmatrices.fypp

@@ -146,8 +146,8 @@ module stdlib_specialmatrices
          type(tridiagonal_${s1}$_type), intent(in) :: A
             ${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$
             ${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$
-            real(${k1}$), intent(in), optional :: alpha
-            real(${k1}$), intent(in), optional :: beta
+            ${t1}$, intent(in), optional :: alpha
+            ${t1}$, intent(in), optional :: beta
             character(1), intent(in), optional :: op
         end subroutine
         #:endfor

src/stdlib_specialmatrices_tridiagonal.fypp

@@ -155,15 +155,15 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
         type(tridiagonal_${s1}$_type), intent(in) :: A
         ${t1}$, intent(in), contiguous, target :: x${ranksuffix(rank)}$
         ${t1}$, intent(inout), contiguous, target :: y${ranksuffix(rank)}$
-        real(${k1}$), intent(in), optional :: alpha
-        real(${k1}$), intent(in), optional :: beta
+        ${t1}$, intent(in), optional :: alpha
+        ${t1}$, intent(in), optional :: beta
         character(1), intent(in), optional :: op
 
         ! Internal variables.
-        real(${k1}$) :: alpha_, beta_
+        ${t1}$ :: alpha_, beta_
         integer(ilp) :: n, nrhs, ldx, ldy
         character(1) :: op_
-        logical :: alpha_special, beta_special
+        logical :: is_alpha_special, is_beta_special
 
         #:if rank == 1
         ${t1}$, pointer :: xmat(:, :), ymat(:, :)
@@ -174,8 +174,8 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
         beta_  = 0.0_${k1}$ ; if (present(beta))  beta_  = beta
         op_    = "N"        ; if (present(op))    op_    = op
 
-        alpha_special = (alpha_ ==  1.0_${k1}$  .or. alpha_ ==  0.0_${k1}$  .or. alpha_ == -1.0_${k1}$)
-        beta_special  = (beta_  ==  1.0_${k1}$  .or. beta_  ==  0.0_${k1}$  .or. beta_  == -1.0_${k1}$)
+        is_alpha_special = (alpha_ ==  1.0_${k1}$  .or. alpha_ ==  0.0_${k1}$  .or. alpha_ == -1.0_${k1}$)
+        is_beta_special  = (beta_  ==  1.0_${k1}$  .or. beta_  ==  0.0_${k1}$  .or. beta_  == -1.0_${k1}$)
 
         ! Prepare Lapack arguments.
         n = A%n ; ldx = n ; ldy = n ;
@@ -184,18 +184,26 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
         #:if rank == 1
         ! Pointer trick.
         xmat(1:n, 1:nrhs) => x ; ymat(1:n, 1:nrhs) => y
-        if(alpha_special .and. beta_special) then
+        #:if t1.startswith('complex')
+        call glagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy)
+        #:else
+        if(is_alpha_special .and. is_beta_special) then
             call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy)
         else
             call glagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, xmat, ldx, beta_, ymat, ldy)
         end if
+        #:endif
         #:else
-        if(alpha_special .and. beta_special) then
+        #:if t1.startswith('complex')
+        call glagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy)
+        #:else
+        if(is_alpha_special .and. is_beta_special) then
             call lagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy)
         else
             call glagtm(op_, n, nrhs, alpha_, A%dl, A%dv, A%du, x, ldx, beta_, y, ldy)
         end if
         #:endif
+        #:endif
     end subroutine
     #:endfor
     #:endfor