# Stdlib: Solving a system of equations with a matrix right hand side

**URL:** <https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660>\
**Category:** stdlib\
**Created:** [May 3, 2025, 1:27pm UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660 "2025-05-03T13:27:58Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![general\_rishkin](https://avatars.discourse-cdn.com/v4/letter/g/a9a28c/32.png) [@general\_rishkin](https://fortran-lang.discourse.group/u/general_rishkin)\
**Post date:** [May 3, 2025, 1:27pm UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/1 "2025-05-03T13:27:58Z")

</div>

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 )
> | 1
> 
> ```
> 
> Error: 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?

---

<div class="post-metadata">

**Author:** ![fxm](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/fxm/32/5726_2.png) [@fxm](https://fortran-lang.discourse.group/u/fxm)\
**Post date:** [May 3, 2025, 1:35pm UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/2 "2025-05-03T13:35:18Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [May 3, 2025, 3:09pm UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/3 "2025-05-03T15:09:28Z")

</div>

> [@general\_rishkin](#):
>
> `call solve_lu( (matmul(B, D)), B )`

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.

---

<div class="post-metadata">

**Author:** ![general\_rishkin](https://avatars.discourse-cdn.com/v4/letter/g/a9a28c/32.png) [@general\_rishkin](https://fortran-lang.discourse.group/u/general_rishkin)\
**Post date:** [May 3, 2025, 3:11pm UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/4 "2025-05-03T15:11:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\
**Post date:** [May 5, 2025, 7:08am UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/5 "2025-05-05T07:08:25Z")

</div>

@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

```Fortran
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:

```Fortran
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](https://github.com/fortran-lang/stdlib/pull/951) and I can’t wait to see what he can do on the interface side: it would be possible to do

```auto
matmul(transposed(A),B)
matmul(A,hermite(B))
...

```

without performing any allocation or data motion by wrapping against `GEMM` internally.

---

<div class="post-metadata">

**Author:** ![suprit05](https://avatars.discourse-cdn.com/v4/letter/s/e5b9ba/32.png) [@suprit05](https://fortran-lang.discourse.group/u/suprit05)\
**Post date:** [May 5, 2025, 10:10am UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/6 "2025-05-05T10:10:50Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [May 5, 2025, 3:22pm UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/7 "2025-05-05T15:22:28Z")

</div>

> [@FedericoPerini](#):
>
> However, that mandates an `intent(inout)` argument also when `.not.overwrite_A`.

This documentation page might need to be updated to be consistent. [linalg – Fortran-lang/stdlib](https://stdlib.fortran-lang.org/page/specs/stdlib_linalg.html#arguments_15)

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

> [@FedericoPerini](#):
>
> ```auto
> matmul(transposed(A),B)
> matmul(A,hermite(B))
> 
> ```

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.

> [@FedericoPerini](#):
>
> without performing any allocation or data motion by wrapping against `GEMM` internally.

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?

---

<div class="post-metadata">

**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\
**Post date:** [May 6, 2025, 6:44am UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/8 "2025-05-06T06:44:34Z")

</div>

> [@RonShepard](#):
>
> 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.

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)

```fortran
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

```fortran
! 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:

```fortran
C = stdlib_matmul(A, B)
C = stdlib_matmul(transposed(A), B)
C = stlib_matmul(A, hermitian(C))

```

could also be an operator:

```fortran
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

---

<div class="post-metadata">

**Author:** ![rsci](https://avatars.discourse-cdn.com/v4/letter/r/e480ec/32.png) [@rsci](https://fortran-lang.discourse.group/u/rsci)\
**Post date:** [May 6, 2025, 11:43am UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/9 "2025-05-06T11:43:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\
**Post date:** [May 6, 2025, 12:15pm UTC](https://fortran-lang.discourse.group/t/stdlib-solving-a-system-of-equations-with-a-matrix-right-hand-side/9660/10 "2025-05-06T12:15:08Z")

</div>

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:

- module exports the generation interfaces: `interface transposed` and `interface hermitian`
- the `matrix_state_${rn}$` types are private to the module, hence they can’t be used anywhere else but the `matmul` api