Use BLAS and LAPACK

[ivan-pi]

Currently, BFPack uses a few linear algebra primitives including:

• FINDinv (inverse of square matrix)
• spofa (factorization of a real symmetric positive definite matrix)

These can be replaced with procedures from libraries that are already shipped with R (libRblas.so, libRlapack.so).

?pofa can be replaced with ?potrf, simply as

call dpofa(a,lda,n,info) ! <- replace this
call dpotrf('U',n,a,lda,info)

The routine FINDinv can be implemented with the LAPACK procedures, dgetrf and dgetri, as follows:

! Subroutine to find the inverse of a square matrix
SUBROUTINE FINDinv(matrix, inverse, n, errorflag)

    implicit none

    !Declarations
    INTEGER(rint), INTENT(IN) :: n
    REAL(rdp), INTENT(IN) :: matrix(n,n)  !Input matrix
    INTEGER(rint), INTENT(OUT) :: errorflag  !Return error status. -1 for error, 0 for normal
    REAL(rdp), INTENT(OUT) :: inverse(n,n) !Inverted matrix

    integer :: ipiv(n), info
    real(rdp) :: work(n), lwork

    external :: dgetrf, dgetri

    errorflag = 0

    inverse = matrix
    call dgetrf(n,n,inverse,n,ipiv,info)
    if (info > 0) then
        inverse = 0
        errorflag = -1
        return
    end if

    lwork = n
    call dgetri(n,inverse,n,ipiv,work,lwork,info)
    if (info > 0) then
        inverse = 0
        errorflag = -1
        return
    end if

END SUBROUTINE FINDinv

[jomulder]

Hi Ivan, thanks a lot for checking this and for these suggestions. I will change this and will also check whether other primitive routines could be replaced. Best, Joris

[jomulder]

When using this FINDinv subroutine, I do get errors from BLAS or LAPACK on examples that did run fine previously

** On entry to DGETRI, parameter number 6 had an illegal value

I will google (and debug) a bit to see what's causing this.

[ivan-pi]

Ooops, that is my mistake. The 6th parameter is supposed to be an integer:

    integer :: ipiv(n), info, lwork
    real(rdp) :: work(n)

(Classic example of why external procedures are to be avoided, because they don't support argument type-checking.)

[jomulder]

Thanks again Ivan!

Just for my understanding, this statement

external :: dgetrf, dgetri

is required to 'tell' that these subroutines are external? Should I also add this for the other recommended subroutine dpotrf? It does seem to work fine without such a statement for dpotrf. And about your comment that external procedures should be avoided, these BLAS/LAPACK procedures are external (right?) but still advisable then. So is there an alternative to these external procedures? Just for me that I understand the rational.

[ivan-pi]

Yes, I'd also recommend this for dpotrf. As you say it works without it, because a Fortran processor/compiler implicitly assumes that it's calling an external procedure. In Fortran 2018, the implicit none statement was amended, and there is now a version,

implicit none (type, external)

If added, the compiler will complain about any procedures that don't have an explicit interface (e.g. they come from a module, have an interface block, or are an internal procedure) or aren't marked as external.

Usually you will see,

implicit none    ! same as implicit none (type)

that only complains about variables without a type declaration.

Ideally LAPACK procedures would be provided in a module, but historically this wasn't the case (originally it was written in F77 which didn't have modules). Today this would cause a huge portability problem, as most LAPACK libraries (libRlapack.so, Intel MKL, OpenBLAS, ...) assume the procedures are external (i.e. they don't exist in a module).

The best practice is to provide an interface block manually. In the scopes where the LAPACK procedures are used, you would insert:

interface
  subroutine dpotrf(uplo,n,a,lda,info)
     character :: uplo
     integer :: n, lda, info
     double precision :: a(lda,*)
  end subroutine

   subroutine dgetrf(m,n,a,lda,ipiv,info)
      integer :: m, n, lda, info
      double precision :: a(lda,*)
      integer :: ipiv(*)
   end subroutine

   subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
      integer :: n, lda, lwork, info
      double precision :: a(lda,*), work(*)
      integer :: ipiv(*)
   end subroutine
end interface

When providing this interface block you can "cheat" a little by adding intent attributes (they do not impact the linker). You could also place these interface blocks in a module or an include file and reuse them.

[ivan-pi]

I'm glad this seems to work at least. Would you like to try and move the linear algebra routines and interfaces into a separate module, shared by the other procedures?

[jomulder]

Well, in the end I only need the main subroutine in bct_mixedordinal (bct_continuous is a special case), and the others at simple tests. So therefore I don't think it is necessary to move the LA routines and interfaces to a module. But I understand now that this would be preferred when having more subroutines.