I am trying to replicate this MATLAB code: P = (B*D) \ B.
B and D are matrices.
Using stdlib, I do:
call solve_lu( matmul(B, D), B )
I get the error:
| call solve_lu( (matmul(B, D)), B ) | 1Error: There is no specific subroutine for the generic ‘solve_lu’ at (1)
Am I using the wrong subroutine from stdlib. If so, which is the correct one to call?
Based on the error, it could be related to the args. You may want to double check the details of arguments the subroutine takes against what you’re actually passing, and pay attention to data types.
In that statement, matmul(B,D) is an expression and it cannot be modified within the subroutine. The corresponding dummy argument of all the specific subroutines has the intent(inout) attribute. The compiler cannot find a specific subroutine that has an intent(in) argument that matches your expression actual argument, so it prints that error message.
If you were calling a specific subroutine, the compiler would probably have given a more specific error message about that argument mismatch.
edit: I see that you have now switched to the function interface, solve(). That interface does have the intent(in) attribute for the matrix, so that should work alright with the expression argument.
edit2: The documentation says that intent(in) for the array is allowed for solve_lu(), so this attribute mismatch probably was not the problem after all.
My mistake. In one of my previous lines of code, I had accidentally introduced a type
o, so I was getting errors calling solve or solve_lu.
I have actually just done: result = solve( BD, B), which works fine.
It was failing before due to the typo.
@RonShepard’s solution is the correct one. Avoiding allocations was perhaps the most frequent request for the stdlib linear algebra API, so I implemented the overwrite_A (like in SciPy API) optional keyword that overwrites the matrix to spare an allocation internally. However, that mandates an intent(inout) argument also when .not.overwrite_A.
To solve this we implemented a pure interface - only for the function version with no further arguments - that you can simply call as
x = solve(matmul(z,t), b)
For the subroutine interface (or the function interface with optional arguments), you need a variable and not an expression:
x = solve(zt, b, overwrite_A=.false., err=err)
call solve(zt, b, x)
@suprit05 is working on his GSoC PR to implement an interesting matmul wrapper and I can’t wait to see what he can do on the interface side: it would be possible to do
matmul(transposed(A),B)
matmul(A,hermite(B))
...
without performing any allocation or data motion by wrapping against GEMM internally.
I got caught up in my tests and final exams, but once this semester ends (in a few days), I’ll try my best to complete the above mentioned PR and also the interface side of it 
.. regardless of my GSoC selection (the result is on May 8th).
This documentation page might need to be updated to be consistent. linalg – Fortran-lang/stdlib
edit: BTW, the only way I know to allow both expressions and variables as actual arguments is to remove the intent() attribute from the dummy argument. This shifts the burden onto the programmer to correctly specify all of the various combinations of the overwrite optional arguments with the actual arguments. The compiler can no longer detect the mismatches at compile time (as occurred in this case).
There have been previous discussions about modifying the standard to include intrinsic functions that perform operations like shallow transpose, shallow complex conjugate, and shallow hermitian conjugate. For 2D matrices, it seems like it should be possible to perform these types of operations just by modifying the array metadata, with no copying or even access to the actual memory storage required. I’m unsure if these operations can be generalized to multidimensional arrays, and they probably can only be applied to assumed shape dummy arguments, but it does seem like these operations would be useful. A side benefit is that fortran arrays could be stored in row-major order to match other language APIs. Perhaps this effort could be used to push this idea over the finish line.
The level-3 BLAS *GEMM() subroutines require each array to have one dimension with stride 1, so it only can be used for a subset of fortran assumed shape arrays, even if the above shallow operations were available. Is there any ongoing effort to generalize the level-3 BLAS interface in this respect?
For N-D, maybe there is a way, but probably the best approach would be a pure-loop implementation, in other words, let the compiler optimize as much as possible. Regarding 2D+GEMM, the idea was actually to create a thin wrapper of the operation via a derived type: (to be templated of course)
type :: matrix_state_rdp
real(dp), pointer :: A(:,:) => null()
character(1) :: Astate = 'N'
end type matrix_state_rdp
interface transposed
module procedure transposed_new_rdp
...
end interface transposed
type(matrix_state_rdp) function transposed_new_rdp(A) result(AT)
real(dp), intent(inout), target :: A(:,:)
AT%A => A
AT%Astate = 'T'
end function
Then a templated base interface for matmul would look like
! Work with normal matrices
${t1}$ function stdlib_matmul(A, B) result(C)
${t1}$, intent(in) :: A(:,:), B(:,:)
! Work with transposed/hermitian swaps
${rt}$ function stdlib_matmul(A, B) result(C)
matrix_state_${rn}$, intent(in) :: A, B
So the user writing code would have it clear:
C = stdlib_matmul(A, B)
C = stdlib_matmul(transposed(A), B)
C = stlib_matmul(A, hermitian(C))
could also be an operator:
C = strlib_matmul(A, B)
C = stdlib_matmul(.t.A, B)
C = stlib_matmul(A, .h.C)
without it triggering any actual data movement, just wrapping the requested matrix state in the derived type
I’m wondering whether the use of AT%A outside the transposed_new_rdp function is standard-conforming? I mean, nothing is guaranteed (?) about the status of the AT%A pointer outside the scope of the transposed_new_rdp, unless matrix A has target attribute at the parent scope, where the transposed and matmul are called.
That’s absolutely correct: due to the usual pointer/target restrictions, the pointer may become undefined upstream of the scope matmul is called from, even in present of intent(inout), target.
A solution to avoid this danger imho is to restrict usage of these wrappers to matmul only:
interface transposed and interface hermitianmatrix_state_${rn}$ types are private to the module, hence they can’t be used anywhere else but the matmul api