Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Sign up
Appearance settings

Commit 3bdcc82

Browse files
authored
Merge pull request #806 from perazz/linalg_solve
linalg: solve
2 parents 60ce18f + 5832df5 commit 3bdcc82

File tree

11 files changed

+653
-3
lines changed

11 files changed

+653
-3
lines changed

‎doc/specs/stdlib_linalg.md

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -600,6 +600,107 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l
600600
{!example/linalg/example_is_hessenberg.f90!}
601601
```
602602

603+
## `solve` - Solves a linear matrix equation or a linear system of equations.
604+
605+
### Status
606+
607+
Experimental
608+
609+
### Description
610+
611+
This function computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a square, full-rank, `real` or `complex` matrix.
612+
613+
Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned.
614+
An error is returned if the matrix is rank-deficient or singular to working precision.
615+
The solver is based on LAPACK's `*GESV` backends.
616+
617+
### Syntax
618+
619+
`Pure` interface:
620+
621+
`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b)`
622+
623+
Expert interface:
624+
625+
`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b [, overwrite_a], err)`
626+
627+
### Arguments
628+
629+
`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call.
630+
631+
`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument.
632+
633+
`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
634+
635+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. The function is not `pure` if this argument is provided.
636+
637+
### Return value
638+
639+
For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.
640+
641+
Raises `LINALG_ERROR` if the matrix is singular to working precision.
642+
Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes.
643+
If `err` is not present, exceptions trigger an `error stop`.
644+
645+
### Example
646+
647+
```fortran
648+
{!example/linalg/example_solve1.f90!}
649+
650+
{!example/linalg/example_solve2.f90!}
651+
```
652+
653+
## `solve_lu` - Solves a linear matrix equation or a linear system of equations (subroutine interface).
654+
655+
### Status
656+
657+
Experimental
658+
659+
### Description
660+
661+
This subroutine computes the solution to a linear matrix equation \( A \cdot x = b \), where \( A \) is a square, full-rank, `real` or `complex` matrix.
662+
663+
Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned.
664+
An error is returned if the matrix is rank-deficient or singular to working precision.
665+
If all optional arrays are provided by the user, no internal allocations take place.
666+
The solver is based on LAPACK's `*GESV` backends.
667+
668+
### Syntax
669+
670+
Simple (`Pure`) interface:
671+
672+
`call ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x)`
673+
674+
Expert (`Pure`) interface:
675+
676+
`call ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x [, pivot, overwrite_a, err])`
677+
678+
### Arguments
679+
680+
`a`: Shall be a rank-2 `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call.
681+
682+
`b`: Shall be a rank-1 or rank-2 array of the same kind as `a`, containing the right-hand-side vector(s). It is an `intent(in)` argument.
683+
684+
`x`: Shall be a rank-1 or rank-2 array of the same kind and size as `b`, that returns the solution(s) to the system. It is an `intent(inout)` argument, and must have the `contiguous` property.
685+
686+
`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices.
687+
688+
`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
689+
690+
### Return value
691+
692+
For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.
693+
694+
Raises `LINALG_ERROR` if the matrix is singular to working precision.
695+
Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes.
696+
If `err` is not present, exceptions trigger an `error stop`.
697+
698+
### Example
699+
700+
```fortran
701+
{!example/linalg/example_solve3.f90!}
702+
```
703+
603704
## `lstsq` - Computes the least squares solution to a linear matrix equation.
604705

605706
### Status

‎example/linalg/CMakeLists.txt

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,5 +20,8 @@ ADD_EXAMPLE(blas_gemv)
2020
ADD_EXAMPLE(lapack_getrf)
2121
ADD_EXAMPLE(lstsq1)
2222
ADD_EXAMPLE(lstsq2)
23+
ADD_EXAMPLE(solve1)
24+
ADD_EXAMPLE(solve2)
25+
ADD_EXAMPLE(solve3)
2326
ADD_EXAMPLE(determinant)
2427
ADD_EXAMPLE(determinant2)

‎example/linalg/example_solve1.f90

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
program example_solve1
2+
use stdlib_linalg_constants, only: sp
3+
use stdlib_linalg, only: solve, linalg_state_type
4+
implicit none
5+
6+
real(sp), allocatable :: A(:,:),b(:),x(:)
7+
8+
! Solve a system of 3 linear equations:
9+
! 4x + 3y + 2z = 25
10+
! -2x + 2y + 3z = -10
11+
! 3x - 5y + 2z = -4
12+
13+
! Note: Fortran is column-major! -> transpose
14+
A = transpose(reshape([ 4, 3, 2, &
15+
-2, 2, 3, &
16+
3,-5, 2], [3,3]))
17+
b = [25,-10,-4]
18+
19+
! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
20+
x = solve(A,b)
21+
22+
print *, 'solution: ',x
23+
! 5.0, 3.0, -2.0
24+
25+
end program example_solve1
26+

‎example/linalg/example_solve2.f90

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
program example_solve2
2+
use stdlib_linalg_constants, only: sp
3+
use stdlib_linalg, only: solve, linalg_state_type
4+
implicit none
5+
6+
complex(sp), allocatable :: A(:,:),b(:),x(:)
7+
8+
! Solve a system of 3 complex linear equations:
9+
! 2x + iy + 2z = (5-i)
10+
! -ix + (4-3i)y + 6z = i
11+
! 4x + 3y + z = 1
12+
13+
! Note: Fortran is column-major! -> transpose
14+
A = transpose(reshape([(2.0, 0.0),(0.0, 1.0),(2.0,0.0), &
15+
(0.0,-1.0),(4.0,-3.0),(6.0,0.0), &
16+
(4.0, 0.0),(3.0, 0.0),(1.0,0.0)] , [3,3]))
17+
b = [(5.0,-1.0),(0.0,1.0),(1.0,0.0)]
18+
19+
! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
20+
x = solve(A,b)
21+
22+
print *, 'solution: ',x
23+
! (1.0947,0.3674) (-1.519,-0.4539) (1.1784,-0.1078)
24+
25+
end program example_solve2
26+

‎example/linalg/example_solve3.f90

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
program example_solve3
2+
use stdlib_linalg_constants, only: sp,ilp
3+
use stdlib_linalg, only: solve_lu, linalg_state_type
4+
implicit none
5+
6+
integer(ilp) :: test
7+
integer(ilp), allocatable :: pivot(:)
8+
complex(sp), allocatable :: A(:,:),b(:),x(:)
9+
10+
! Solve a system of 3 complex linear equations:
11+
! 2x + iy + 2z = (5-i)
12+
! -ix + (4-3i)y + 6z = i
13+
! 4x + 3y + z = 1
14+
15+
! Note: Fortran is column-major! -> transpose
16+
A = transpose(reshape([(2.0, 0.0),(0.0, 1.0),(2.0,0.0), &
17+
(0.0,-1.0),(4.0,-3.0),(6.0,0.0), &
18+
(4.0, 0.0),(3.0, 0.0),(1.0,0.0)] , [3,3]))
19+
20+
! Pre-allocate x
21+
allocate(b(size(A,2)),pivot(size(A,2)))
22+
allocate(x,mold=b)
23+
24+
! Call system many times avoiding reallocation
25+
do test=1,100
26+
b = test*[(5.0,-1.0),(0.0,1.0),(1.0,0.0)]
27+
call solve_lu(A,b,x,pivot)
28+
print "(i3,'-th solution: ',*(1x,f12.6))", test,x
29+
end do
30+
31+
end program example_solve3
32+

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ set(fppFiles
2525
stdlib_linalg_outer_product.fypp
2626
stdlib_linalg_kronecker.fypp
2727
stdlib_linalg_cross_product.fypp
28+
stdlib_linalg_solve.fypp
2829
stdlib_linalg_determinant.fypp
2930
stdlib_linalg_state.fypp
3031
stdlib_optval.fypp

‎src/stdlib_linalg.fypp

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,8 @@ module stdlib_linalg
2222
public :: eye
2323
public :: lstsq
2424
public :: lstsq_space
25+
public :: solve
26+
public :: solve_lu
2527
public :: solve_lstsq
2628
public :: trace
2729
public :: outer_product
@@ -228,6 +230,100 @@ module stdlib_linalg
228230
#:endfor
229231
end interface is_hessenberg
230232

233+
! Solve linear system system Ax=b.
234+
interface solve
235+
!! version: experimental
236+
!!
237+
!! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \).
238+
!! ([Specification](../page/specs/stdlib_linalg.html#solve-solves-a-linear-matrix-equation-or-a-linear-system-of-equations))
239+
!!
240+
!!### Summary
241+
!! Interface for solving a linear system arising from a general matrix.
242+
!!
243+
!!### Description
244+
!!
245+
!! This interface provides methods for computing the solution of a linear matrix system.
246+
!! Supported data types include `real` and `complex`. No assumption is made on the matrix
247+
!! structure.
248+
!! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`)
249+
!! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
250+
!!
251+
!!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
252+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
253+
!!
254+
#:for nd,ndsuf,nde in ALL_RHS
255+
#:for rk,rt,ri in RC_KINDS_TYPES
256+
#:if rk!="xdp"
257+
module function stdlib_linalg_${ri}$_solve_${ndsuf}$(a,b,overwrite_a,err) result(x)
258+
!> Input matrix a[n,n]
259+
${rt},ドル intent(inout), target :: a(:,:)
260+
!> Right hand side vector or array, b[n] or b[n,nrhs]
261+
${rt},ドル intent(in) :: b${nd}$
262+
!> [optional] Can A data be overwritten and destroyed?
263+
logical(lk), optional, intent(in) :: overwrite_a
264+
!> [optional] state return flag. On error if not requested, the code will stop
265+
type(linalg_state_type), intent(out) :: err
266+
!> Result array/matrix x[n] or x[n,nrhs]
267+
${rt},ドル allocatable, target :: x${nd}$
268+
end function stdlib_linalg_${ri}$_solve_${ndsuf}$
269+
pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x)
270+
!> Input matrix a[n,n]
271+
${rt},ドル intent(in) :: a(:,:)
272+
!> Right hand side vector or array, b[n] or b[n,nrhs]
273+
${rt},ドル intent(in) :: b${nd}$
274+
!> Result array/matrix x[n] or x[n,nrhs]
275+
${rt},ドル allocatable, target :: x${nd}$
276+
end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$
277+
#:endif
278+
#:endfor
279+
#:endfor
280+
end interface solve
281+
282+
! Solve linear system Ax = b using LU decomposition (subroutine interface).
283+
interface solve_lu
284+
!! version: experimental
285+
!!
286+
!! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \).
287+
!! ([Specification](../page/specs/stdlib_linalg.html#solve-lu-solves-a-linear-matrix-equation-or-a-linear-system-of-equations-subroutine-interface))
288+
!!
289+
!!### Summary
290+
!! Subroutine interface for solving a linear system using LU decomposition.
291+
!!
292+
!!### Description
293+
!!
294+
!! This interface provides methods for computing the solution of a linear matrix system using
295+
!! a subroutine. Supported data types include `real` and `complex`. No assumption is made on the matrix
296+
!! structure. Preallocated space for the solution vector `x` is user-provided, and it may be provided
297+
!! for the array of pivot indices, `pivot`. If all pre-allocated work spaces are provided, no internal
298+
!! memory allocations take place when using this interface.
299+
!! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`)
300+
!! or several (from a 2-d right-hand-side vector `b(:,:)`) systems.
301+
!!
302+
!!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`.
303+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
304+
!!
305+
#:for nd,ndsuf,nde in ALL_RHS
306+
#:for rk,rt,ri in RC_KINDS_TYPES
307+
#:if rk!="xdp"
308+
pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err)
309+
!> Input matrix a[n,n]
310+
${rt},ドル intent(inout), target :: a(:,:)
311+
!> Right hand side vector or array, b[n] or b[n,nrhs]
312+
${rt},ドル intent(in) :: b${nd}$
313+
!> Result array/matrix x[n] or x[n,nrhs]
314+
${rt},ドル intent(inout), contiguous, target :: x${nd}$
315+
!> [optional] Storage array for the diagonal pivot indices
316+
integer(ilp), optional, intent(inout), target :: pivot(:)
317+
!> [optional] Can A data be overwritten and destroyed?
318+
logical(lk), optional, intent(in) :: overwrite_a
319+
!> [optional] state return flag. On error if not requested, the code will stop
320+
type(linalg_state_type), optional, intent(out) :: err
321+
end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$
322+
#:endif
323+
#:endfor
324+
#:endfor
325+
end interface solve_lu
326+
231327
! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized.
232328
interface lstsq
233329
!! version: experimental

0 commit comments

Comments
(0)

AltStyle によって変換されたページ (->オリジナル) /