Upper or lower triangular part of an array

[ivan-pi]

triu and tril

tril - lower triangular part of an array
triu - upper triangular part of an array

Return a copy of the lower/upper triangular part of a rank-2 array. The elements below/above the k-th diagonal are replaced with zeroes (default k=0)

Useful to recover the lower or upper part of a matrix factorization.

Interface

interface tril
    module function tril_rsp(A) result(L)
        real(sp), intent(in) :: A(:,:)
        real(sp) :: L(size(A,1),size(A,2))
    end function
    module function tril_k_rsp(A,k) result(L)
        real(sp), intent(in) :: A(:,:)
        integer, intent(in) :: k
        real(sp) :: L(size(A,1),size(A,2))
    end function
    ! .. repeat for all real and integer kinds ..
end interface

Analogous interface for triu. The interfaces would go to the file stdlib_experimental_linalg.f90. Implementations would go to the submodule stdlib_experimental_linalg_trilu.f90.

Point for discussion: two separate functions (without or with diagonal) as shown above or only a single interface using the present intrinsic?

Other languages

Julia

• tril: https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.tril
• triu: https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.triu

MATLAB

• tril: https://de.mathworks.com/help/matlab/ref/tril.html
• triu: https://de.mathworks.com/help/matlab/ref/triu.html

Python (NumPy)

• tril: https://numpy.org/doc/stable/reference/generated/numpy.tril.html
• triu: https://numpy.org/doc/stable/reference/generated/numpy.triu.html

C++

• (Eigen) Triangular view: https://eigen.tuxfamily.org/dox/group__QuickRefPage.html#title14 (since matrices are classes, this gives a triangular view of a dense matrix with specialized operations)
• (Armadillo) trimatu/trimatl: http://arma.sourceforge.net/docs.html#trimat

Other

• Would we like a subroutine version which works in-place? In Julia they use the tril!(M) and triu!(M) syntax for this purpose.

[ivan-pi]

I was working on some Kalman filter stuff today, and I realized this would be useful to check I am calling the right sequence of BLAS and LAPACK routines.

[jvdp1]

Nice proposition. It is something I often use in Octave.

Point for discussion: two separate functions (without or with diagonal) as shown above or only a single interface using the present intrinsic?

Why would you like 2 separate functions? I would suggest to implement only one function and use the optval function for the optional integer.

• Would we like a subroutine version which works in-place? In Julia they use the tril!(M) and triu!(M) syntax for this purpose.

Yes, I would. However, I don't really like the sign "!" used in Julia.

[certik]

Great idea, thanks.

Regarding in-place, see #177 for a discussion about the syntax. We can't use ! in Fortran (plus I don't really like it anyway).