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          ?    41267           
=========================================
  Hits              ?    58858           
  Misses            ?   175343           
  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?

[jalvesz]

I know that the brackets syntax looks elegant but it does 2 allocations plus an assignment. Preferer the more efficient syntax:

allocate( A%dl(n-1), source = dl )
allocate( A%dv(n), source= dv )
allocate( A%du(n-1), source = du )

[jalvesz]

Please use two lines for this. calls should always be in their own line. Same for array assignments. The use of ";" to mix in a single line multiple executions is recommended only for very basic scalar operations.

[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.