diff --git a/src/stdlib_specialmatrices_tridiagonal.fypp b/src/stdlib_specialmatrices_tridiagonal.fypp index 8fabaa5df..3e7538896 100644 --- a/src/stdlib_specialmatrices_tridiagonal.fypp +++ b/src/stdlib_specialmatrices_tridiagonal.fypp @@ -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. + A%dl = [(dl, i = 1, n-1)] + A%dv = [(dv, i = 1, n)] + A%du = [(du, i = 1, n-1)] + 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 diff --git a/test/linalg/test_linalg_specialmatrices.fypp b/test/linalg/test_linalg_specialmatrices.fypp index c445ec0da..505cdfc22 100644 --- a/test/linalg/test_linalg_specialmatrices.fypp +++ b/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,27 @@ 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 +109,7 @@ contains call check(error, state%ok(), .false.) if (allocated(error)) return end block - #:endfor + #:endfor end subroutine end module