Proposal for linalg

[certik]

eig, eigh, inv, solve, det, svd, ... . A possible implementation is here: https://github.com/certik/fortran-utils/blob/b43bd24cd421509a5bc6d3b9c3eeae8ce856ed88/src/linalg.f90.

All these functions will be implemented in stdlib_linalg module, and they would probably just call Lapack. The general idea of these routines is to be general routines that will just work, with a simple intuitive interface, and the highest performance given the simple API. One can always achieve higher performance with more specialized routines for a particular problem (and more complicated API), but that is not the point here. Rather we would like a Matlab / NumPy style routines to do linear algebra.

In particular, let's start with eig, for an eigenvalue and eigenvectors of a general (non-symmetric) matrix. Both NumPy and Matlab have a very similar interface called eig, that I propose we use:

• NumPy: https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html
• Matlab: https://www.mathworks.com/help/matlab/ref/eig.html

Julia seems to have more of an "object oriented" interface called eigen: https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/index.html#LinearAlgebra.eigen, which uses some Julia language features to emulate the Matlab style vals, vecs = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0]).

[certik]

Let's discuss the API for eig. Fortran cannot return more than one result in a function, and eig must return both eigenvalues and eigenvectors. So the only solution is a subroutine. The first argument should be the input matrix, to be consistent with NumPy and Matlab. The second argument can be an optional B matrix for the generalized eigenvector problem. The next argument can be eigenvalues and the last argument eigenvectors (NumPy, Julia style) or they can be switched (Matlab style). For efficiency reasons the eigenvalues should be just a vector (NumPy, Julia style) not a diagonal matrix (Matlab style: the default is a matrix, but the eigvalOption input option can specify vector to return a vector as well).

That reasoning leads to the following API:

  interface eig
     module procedure seig
     module procedure seig_generalized
     module procedure deig
     module procedure deig_generalized
     module procedure zeig
     module procedure zeig_generalized
  end interface eig

Here s means single precision, d double, z complex (double). Here is an example of a double precision interface:

    subroutine deig(A, lam, c)
    ! Solves a standard eigenproblem A*c = lam*c
    real(dp), intent(in) :: A(:, :)  ! matrix A
    complex(dp), intent(out) :: lam(:)  ! eigenvalues: A c = lam c
    complex(dp), intent(out) :: c(:, :)  ! eigenvectors: A c = lam c; c(i,j) = ith component of jth vec.
    ...
    end subroutine

    subroutine deig_generalized(A, B, lam, c)
    ! Solves a generalized eigenproblem A*c = lam*B*c
    real(dp), intent(in) :: A(:, :)  ! matrix A
    real(dp), intent(in) :: B(:, :)  ! matrix B
    complex(dp), intent(out) :: lam(:)  ! eigenvalues: A c = lam c
    complex(dp), intent(out) :: c(:, :)  ! eigenvectors: A c = lam c; c(i,j) = ith component of jth vec.
    ...
    end subroutine

[certik]

Here are some common points that should be discussed:

• Should the return arrays be allocatable, intent(out)? I would say no, for efficiency reasons.

• Should there be some optional options such as algorithm in the Matlab's API? I think there could.

• We need to implement all combinations of types. The above proposal assumes real input matrices. We should also allow a complex input matrix.

• Modify inplace? Lapack sometimes returns the result by overwriting the input matrix, although in the case of eig that does not seem to be the case (it is the case for cholesky), but the input matrix is destroyed, so internally there is an extra copy happening. The simple solution is the above proposal, that requires an extra copy, but the input matrix is never modified. Another solution is to have two kinds of eig subroutines, the one above, and another one called eig_inplace which will modify the matrix A and return the eigenvalues there.

[certik]

@cmacmackin, @milancurcic, @septcolor what do you think of this proposal? I chose it so that it's something we might be able to agree quickly, and then we can figure out the proper workflow.

[milancurcic]

I'm not experienced with linalg stuff, but all looks reasonable. Procedure names are short and intuitive.

I agree that allocatable for intent(out) dummy arguments should be avoided unless necessary.

I don't like APIs that modify variables in-place, so I suggest we make an extra copy in default implementation, but if desired by the community we can also have a variant with an intent(in out) argument as you suggested.

Otherwise, full steam ahead as far as I'm concerned.

[ivan-pi]

From my experience with using LAPACK routines in some fluid-solvers (immersed boundary method, lattice boltzmann equation, etc) and PDE routines (orthogonal collocation), the right hand-side is rarely required after the solution of the system. If it is needed the user can always make a copy beforehand. Would the A matrix be returned factorize or would there be another copy?

real, allocatable(:) :: x, b
real, allocatable(:) :: A
A = reshape([1,3,2,4],[2,2])
b = [5, 6]
x = b  ! user creates the copy
call solve(A,x) ! solve A*x=b

If you follow a more functional style you could write

x = solve(A,b)

The annoying thing is that a subroutine and a function cannot be in the same generic interface block (perhaps something to discuss in the J3 Fortran proposal), meaning we would need two routines with different names.

While I like the idea to have convenience routines, I think it is also worth considering a more object-oriented interface, similar to the one in the Eigen library: http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html. A major difference compared to Eigen is that in Fortran we would like to use native real arrays and not some abstract matrix data types, meaning the solvers would likely become the abstract data type, e.g.:

type(linalg_solver) :: solver
! ... create A and b ...
call solver%init(A,factorization="LU") 
x = solver%solve(b)

[certik]

@milancurcic, @ivan-pi I think you just demonstrated that we need both "copy" and "inplace" API. Julia does that, they append ! after the name of functions that modify arguments inplace, e.g. cholesky vs cholesky!. Fortran does not allow ! as part of the name, so we need some other convention. Any ideas?

Finally, @ivan-pi also has a good point that besides the direct "copy" and "inplace" API, some people might also prefer an object oriented API. So I think we might need 3 APIs.

@milancurcic so I can submit a PR for eig. But I would like new code to land in some kind of "experimental" namespace in stdlib. So that we can start using it and trying it out, but to reserve a right to change the API, if we discover issues. Only after some time and real usage in real codes, we should consider moving a functionality from experimental to the actual stdlib, after which we must support the API in a backwards compatible way.

[certik]

@jacobwilliams, @marshallward do you have any thoughts on this workflow? Let's just concentrate on eig for now. We'll tackle solve and others later.

[cmacmackin]

Personally, my preference would be for using derived types to store the matrix and/or its factorisation. With LAPACK, the factorisation requires more information than can be stored in the original matrix (e.g., pivot information), requiring additional arrays to be passed in. There are use-cases that reuse the factorisation, so we don't want that information to be lost. An example of an OOP approach can be found here: https://github.com/cmacmackin/OOP-Fortran-Examples/blob/master/04-lapack-wrapper.f90.

In terms of returning multiple results, another approach, aside from using a subroutine, would be to define a transparent derived type with these results as components and just to return that.

[milancurcic]

@milancurcic, @ivan-pi I think you just demonstrated that we need both "copy" and "inplace" API.

At best, I'd say we may need both. Let's focus on one, and consider another later.

I'd caution against involving derived types or OOP here from the get-go. Obviously there are implementations with procedures only. I suggest taking an incremental approach here:

• Implement the most basic building block, say, eig that doesn't modify in-place (makes a copy);
• Evaluate and hear from community if in-place variant is really needed. It's likely to be more computationally efficient, albeit slightly less elegant. Does the community want it? Great! Let's do it;
• Evaluate the basic 2 implementations and discuss if higher-level abstractions (derived types) are needed.

In summary, let's make a banana first, and jungle later. Gorilla will come for the banana sooner or later. Let's worry about it then.

[milancurcic]

Fortran does not allow ! as part of the name, so we need some other convention. Any ideas?

In case of eig, it seems that both variants would need to be subroutines. The number of arguments will be different for the "copy" and "in-place" specific procedures. Can we then use a generic interface for both?

[certik]

I agree with your plan to start with eig, then possibly an inplace variant, and then possibly an OO interface.