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

[Mahmood-Sinan]

Summary:
This pr fixes two issues within tridiagonal_matrices module.

1. Inside the impure functions for initializing tridiagonal_matrices, the matrix fields were assigned even if there was an error(for example, when the sizes of du and dl vectors differed from dv by more than one).
2. Inside spmv subroutine, the y vector was forcefully being zeroed instead of using the passed argument.

Changes:

1. Added an if (err0%ok()) check in the impure initialization functions.
2. Removed the unintended y = 0.0_${k1}$ in the spmv subroutine.

[loiseaujc]

Thanks for spotting some things I completely missed. Could you also extend the tests to ensure that the cases with $(\alpha, \beta) \neq (1, 0)$ work correctly.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
⚠️ Please upload report for BASE (master@21c65ab). Learn more about missing BASE report.

Additional details and impacted files

@@            Coverage Diff            @@
##             master    #1054   +/-   ##
=========================================
  Coverage          ?   25.13%           
=========================================
  Files             ?      570           
  Lines             ?   234201           
  Branches          ?    41277           
=========================================
  Hits              ?    58869           
  Misses            ?   175332           
  Partials          ?        0           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[loiseaujc]

As it is, we do not perform any check in spmv to ensure $(\alpha, \beta) \in (-1, 0, 1)^2$ nor that op = N, T or C. I was wondering whether we should include such checks to dumb-proof the implementation, in particular the $(\alpha, \beta)$ one since this is a quirk from the lapack implementation of *LAGTM differing from *GEMM.

@perazz, @jvdp1 and @jalvesz : What do you think? Should we modify the interface to be something like

subroutine spmv(A, x, y, alpha, beta, op, err)

where err is a linalg_state_type and have the subroutine potentially raise an error if the alpha, beta or op provided by the user are not admissible? Note that *GEMM does not have an info return flag. I'm not entirely sure what it does if transa or transb are not admissible, neither does *LAGTM. Also, the spmv kernels by @jalvesz for sparse matrices do not ensure that op is admissible.

[jalvesz]

I guess we could add an optional error output argument to check for the admissibility of op. Regarding alpha and beta, they could in principle take on any real (or complex) value, so I don't think it is necessary to be so strict as to check on their value.

[loiseaujc]

Well the thing is the underlying spmv implementation relies on *LAGTM. For reasons I don't quite understand, the docs reads

!>
!> ZLAGTM performs a matrix-matrix product of the form
!>
!>    B := alpha * A * X + beta * B
!>
!> where A is a tridiagonal matrix of order N, B and X are N by NRHS
!> matrices, and alpha and beta are real scalars, each of which may be
!> 0., 1., or -1.
!> 

Hence, compared to the gemm or the sparse spmv you implemented, lagtm has this oddity with alpha and beta. There probably is a historical/practical reason for that but I don't get it.

[jalvesz]

I saw the reason ... *LAGTM uses the coefficients within if\else switches instead of as actual coefficients, thus the limitation. But this does looks like a legacy oddity, I think this is a deliberate choice of the time for performance reason, not necessarily to build a general purpose routine but actually an specialized one ...

Have you tried to call these routines from OpenBLAS/MKL with values different than -1,0,1 to check how the libraries behave?

[loiseaujc]

I've just tried with openblas 0.3.30 and it follows the lapack documentation. Any value of alpha not in $(-1, 0, 1)$ is treated as 0, and any value of beta not in the same set is treated as 1.

[Mahmood-Sinan]

Thank you for the suggestions. I have made the changes. Please let me know if anything else needs adjustment.

[jalvesz]

Thanks @Mahmood-Sinan for these fixes, LGTM!

@loiseaujc I would propose that this PR is kept as is in terms of scope and then open 2 different ones:
PR1. Add an optional err output to the spmv kernels for sparse and special matrices checking only on the admissibility of the op argument.
PR2. As is, these higher level spmv are supposed to be general with respect to alpha and beta, but indeed blas\lapack only provide specialized ones. I would propose to included within the tridiagonal module an extention to *LAGTM to cover for the general case and that within the spmv kernel a check is carried out to switch between the lapack backed *LAGTM for the -1,0,1 cases, else use the a stdlib internal generic version:

pure subroutine stdlib_glagtm_sp( trans, n, nrhs, alpha, dl, d, du, x, ldx, beta,b, ldb )
     !! GLAGTM performs a matrix-vector product of the form
     !! B := alpha * A * X + beta * B
     !! where A is a tridiagonal matrix of order N, B and X are N by NRHS
     !! matrices, and alpha and beta are real scalars
           use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           integer, parameter :: ilp = int32
           integer, parameter :: wp = sp
           character, intent(in) :: trans
           
           integer(ilp), intent(in) :: ldb, ldx, n, nrhs
           real(wp), intent(in) :: alpha, beta
           ! Array Arguments 
           real(wp), intent(inout) :: b(ldb,*)
           real(wp), intent(in) :: d(*), dl(*), du(*), x(ldx,*)
        ! =====================================================================
           
           ! Local Scalars 
           integer(ilp) :: i, j
           real(wp) :: temp
           ! Executable Statements 
           if( n==0 )return
           ! multiply b by beta if beta/=0
           if( beta==zero ) then
            b(1:n,1:nrhs) = zero
           else
            b(1:n,1:nrhs) = beta * b(1:n,1:nrhs)
           end if

           if( trans(1:1) == 'N' ) then
              ! compute b := b + alpha * a*x
              do j = 1, nrhs
                 if( n==1_ilp ) then
                    temp = d( 1_ilp )*x( 1_ilp, j )
                    b( 1_ilp, j ) = b( 1_ilp, j ) + alpha*temp
                 else
                    temp = d( 1_ilp )*x( 1_ilp, j ) + du( 1_ilp )*x( 2_ilp, j )
                    b( 1_ilp, j ) = b( 1_ilp, j ) + alpha*temp
                    do i = 2, n - 1
                       temp = dl( i-1 )*x( i-1, j ) + d( i )*x( i, j )  + du( i )*x( i+1, j )
                       b( i, j ) = b( i, j ) + alpha*temp 
                    end do
                    temp = dl( n-1 )*x( n-1, j ) + d( n )*x( n, j )
                    b( n, j ) = b( n, j ) + alpha*temp
                 end if
              end do
           else
              ! compute b := b + alpha * a**t*x
              do j = 1, nrhs
                 if( n==1_ilp ) then
                    temp = d( 1_ilp )*x( 1_ilp, j )
                    b( 1_ilp, j ) = b( 1_ilp, j ) + alpha*temp
                 else
                    temp = d( 1_ilp )*x( 1_ilp, j ) + dl( 1_ilp )*x( 2_ilp, j )
                    b( 1_ilp, j ) = b( 1_ilp, j ) + alpha*temp
                    do i = 2, n - 1
                       temp = du( i-1 )*x( i-1, j ) + d( i )*x( i, j ) + dl( i )*x( i+1, j )
                       b( i, j ) = b( i, j ) + alpha*temp
                    end do
                    temp = du( n-1 )*x( n-1, j ) + d( n )*x( n, j )
                    b( n, j ) = b( n, j ) + alpha*temp
                 end if
              end do
           end if
           return
end subroutine stdlib_glagtm_sp

[loiseaujc]

All good for me. @Mahmood-Sinan: do you want to take the lead with the generalized lagtm ?

[Mahmood-Sinan]

Thanks for approving. And @loiseaujc Yes, I’d be happy to work on the generalized LAGTM extension. I’ll likely need some help or pointers as I go.