Fix(tridiagonal_matrices): spmv y vector overwrite and missing error check in impure init functions (#1054)

* Fix: spmv y-reset bug and add error guard in tridiagonal init

* minor change

* test(tridiagonal): added test cases for all alpha and beta combinations under spmv

* use allocate for matrix elements, move call spmv to separate line

Author: Mahmood-Sinan

src/stdlib_specialmatrices_tridiagonal.fypp

@@ -69,9 +69,9 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
             call linalg_error_handling(err0)
         endif
         ! Matrix elements.
-        A%dl = [(dl, i = 1, n-1)]
-        A%dv = [(dv, i = 1, n)]
-        A%du = [(du, i = 1, n-1)]
+        allocate( A%dl(n-1), source = dl )
+        allocate( A%dv(n), source= dv )
+        allocate( A%du(n-1), source = du )
     end function
 
     module function initialize_tridiagonal_impure_${s1}$(dl, dv, du, err) result(A)
@@ -103,10 +103,12 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
             call linalg_error_handling(err0, err)
         endif
 
-        ! Description of the matrix.
-        A%n = n
-        ! Matrix elements.
-        A%dl = dl ; A%dv = dv ; A%du = du
+        if(err0%ok()) then
+            ! Description of the matrix.
+            A%n = n
+            ! Matrix elements.
+            A%dl = dl ; A%dv = dv ; A%du = du
+        endif
     end function
 
     module function initialize_constant_tridiagonal_impure_${s1}$(dl, dv, du, n, err) result(A)
@@ -124,16 +126,19 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
         integer(ilp) :: i
         type(linalg_state_type) :: err0
 
-        ! Description of the matrix.
-        A%n = n
         if (n <= 0) then
             err0 = linalg_state_type(this, LINALG_VALUE_ERROR, "Matrix size needs to be positive, n = ", n, ".")
             call linalg_error_handling(err0, err)
         endif
-        ! Matrix elements.
-        A%dl = [(dl, i = 1, n-1)]
-        A%dv = [(dv, i = 1, n)]
-        A%du = [(du, i = 1, n-1)]
+
+        if(err0%ok()) then
+            ! Description of the matrix.
+            A%n = n
+            ! Matrix elements.
+            allocate( A%dl(n-1), source = dl )
+            allocate( A%dv(n), source= dv )
+            allocate( A%du(n-1), source = du )
+        endif
     end function
     #:endfor
 
@@ -168,7 +173,7 @@ submodule (stdlib_specialmatrices) tridiagonal_matrices
         op_    = "N"        ; if (present(op))    op_    = op
 
         ! Prepare Lapack arguments.
-        n = A%n ; ldx = n ; ldy = n ; y = 0.0_${k1}$
+        n = A%n ; ldx = n ; ldy = n ;
         nrhs = #{if rank==1}# 1 #{else}# size(x, dim=2, kind=ilp) #{endif}#
 
         #:if rank == 1

test/linalg/test_linalg_specialmatrices.fypp

@@ -35,6 +35,10 @@ contains
             ${t1}$, allocatable :: Amat(:,:), dl(:), dv(:), du(:)
             ${t1}$, allocatable :: x(:)
             ${t1}$, allocatable :: y1(:), y2(:)
+            ${t1}$ :: alpha, beta
+
+            integer :: i, j
+            ${t1}$, parameter :: coeffs(3) = [-1.0_wp, 0.0_wp, 1.0_wp]
 
             ! Initialize matrix.
             allocate(dl(n-1), dv(n), du(n-1))
@@ -56,13 +60,28 @@ contains
             call check(error, all_close(y1, y2), .true.)
             if (allocated(error)) return
 
-            #:if t1.startswith('complex') 
+            #:if t1.startswith('complex')
             ! Test y = A.H @ x
             y1 = 0.0_wp ; y2 = 0.0_wp
             y1 = matmul(hermitian(Amat), x) ; call spmv(A, x, y2, op="H")
             call check(error, all_close(y1, y2), .true.)
             if (allocated(error)) return
             #:endif
+
+            ! Test y = alpha * A @ x + beta * y for alpha,beta in {-1,0,1}
+            do i = 1, 3
+                do j = 1,3
+                    alpha = coeffs(i)
+                    beta = coeffs(j)
+
+                    y1 = 0.0_wp
+                    call random_number(y2)
+                    y1 = alpha * matmul(Amat, x) + beta * y2
+                    call spmv(A, x, y2, alpha=alpha, beta=beta)
+                    call check(error, all_close(y1, y2), .true.)
+                    if (allocated(error)) return
+                end do
+            end do
         end block
         #:endfor
     end subroutine
@@ -91,7 +110,7 @@ contains
             call check(error, state%ok(), .false.)
             if (allocated(error)) return
         end block
-        #:endfor       
+        #:endfor
     end subroutine
 
 end module