diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index a081291d0..c414c071a 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -600,6 +600,107 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l {!example/linalg/example_is_hessenberg.f90!} ``` +## `solve` - Solves a linear matrix equation or a linear system of equations. + +### Status + +Experimental + +### Description + +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. + +Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. +An error is returned if the matrix is rank-deficient or singular to working precision. +The solver is based on LAPACK's `*GESV` backends. + +### Syntax + +`Pure` interface: + +`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b)` + +Expert interface: + +`x = ` [[stdlib_linalg(module):solve(interface)]] `(a, b [, overwrite_a], err)` + +### Arguments + +`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. + +`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. + +`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. + +`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. + +### Return value + +For a full-rank matrix, returns an array value that represents the solution to the linear system of equations. + +Raises `LINALG_ERROR` if the matrix is singular to working precision. +Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_solve1.f90!} + +{!example/linalg/example_solve2.f90!} +``` + +## `solve_lu` - Solves a linear matrix equation or a linear system of equations (subroutine interface). + +### Status + +Experimental + +### Description + +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. + +Result vector or array `x` returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. +An error is returned if the matrix is rank-deficient or singular to working precision. +If all optional arrays are provided by the user, no internal allocations take place. +The solver is based on LAPACK's `*GESV` backends. + +### Syntax + +Simple (`Pure`) interface: + +`call ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x)` + +Expert (`Pure`) interface: + +`call ` [[stdlib_linalg(module):solve_lu(interface)]] `(a, b, x [, pivot, overwrite_a, err])` + +### Arguments + +`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. + +`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. + +`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. + +`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. + +`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. + +### Return value + +For a full-rank matrix, returns an array value that represents the solution to the linear system of equations. + +Raises `LINALG_ERROR` if the matrix is singular to working precision. +Raises `LINALG_VALUE_ERROR` if the matrix and rhs vectors have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_solve3.f90!} +``` + ## `lstsq` - Computes the least squares solution to a linear matrix equation. ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 729f117d3..2d80f067c 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -20,5 +20,8 @@ ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) +ADD_EXAMPLE(solve1) +ADD_EXAMPLE(solve2) +ADD_EXAMPLE(solve3) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) diff --git a/example/linalg/example_solve1.f90 b/example/linalg/example_solve1.f90 new file mode 100644 index 000000000..0b5c1bbf1 --- /dev/null +++ b/example/linalg/example_solve1.f90 @@ -0,0 +1,26 @@ +program example_solve1 + use stdlib_linalg_constants, only: sp + use stdlib_linalg, only: solve, linalg_state_type + implicit none + + real(sp), allocatable :: A(:,:),b(:),x(:) + + ! Solve a system of 3 linear equations: + ! 4x + 3y + 2z = 25 + ! -2x + 2y + 3z = -10 + ! 3x - 5y + 2z = -4 + + ! Note: Fortran is column-major! -> transpose + A = transpose(reshape([ 4, 3, 2, & + -2, 2, 3, & + 3,-5, 2], [3,3])) + b = [25,-10,-4] + + ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3) + x = solve(A,b) + + print *, 'solution: ',x + ! 5.0, 3.0, -2.0 + +end program example_solve1 + diff --git a/example/linalg/example_solve2.f90 b/example/linalg/example_solve2.f90 new file mode 100644 index 000000000..a761ae7bd --- /dev/null +++ b/example/linalg/example_solve2.f90 @@ -0,0 +1,26 @@ +program example_solve2 + use stdlib_linalg_constants, only: sp + use stdlib_linalg, only: solve, linalg_state_type + implicit none + + complex(sp), allocatable :: A(:,:),b(:),x(:) + + ! Solve a system of 3 complex linear equations: + ! 2x + iy + 2z = (5-i) + ! -ix + (4-3i)y + 6z = i + ! 4x + 3y + z = 1 + + ! Note: Fortran is column-major! -> transpose + A = transpose(reshape([(2.0, 0.0),(0.0, 1.0),(2.0,0.0), & + (0.0,-1.0),(4.0,-3.0),(6.0,0.0), & + (4.0, 0.0),(3.0, 0.0),(1.0,0.0)] , [3,3])) + b = [(5.0,-1.0),(0.0,1.0),(1.0,0.0)] + + ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3) + x = solve(A,b) + + print *, 'solution: ',x + ! (1.0947,0.3674) (-1.519,-0.4539) (1.1784,-0.1078) + +end program example_solve2 + diff --git a/example/linalg/example_solve3.f90 b/example/linalg/example_solve3.f90 new file mode 100644 index 000000000..ae92ea5da --- /dev/null +++ b/example/linalg/example_solve3.f90 @@ -0,0 +1,32 @@ +program example_solve3 + use stdlib_linalg_constants, only: sp,ilp + use stdlib_linalg, only: solve_lu, linalg_state_type + implicit none + + integer(ilp) :: test + integer(ilp), allocatable :: pivot(:) + complex(sp), allocatable :: A(:,:),b(:),x(:) + + ! Solve a system of 3 complex linear equations: + ! 2x + iy + 2z = (5-i) + ! -ix + (4-3i)y + 6z = i + ! 4x + 3y + z = 1 + + ! Note: Fortran is column-major! -> transpose + A = transpose(reshape([(2.0, 0.0),(0.0, 1.0),(2.0,0.0), & + (0.0,-1.0),(4.0,-3.0),(6.0,0.0), & + (4.0, 0.0),(3.0, 0.0),(1.0,0.0)] , [3,3])) + + ! Pre-allocate x + allocate(b(size(A,2)),pivot(size(A,2))) + allocate(x,mold=b) + + ! Call system many times avoiding reallocation + do test=1,100 + b = test*[(5.0,-1.0),(0.0,1.0),(1.0,0.0)] + call solve_lu(A,b,x,pivot) + print "(i3,'-th solution: ',*(1x,f12.6))", test,x + end do + +end program example_solve3 + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b58ba0b9c..a031ab823 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -25,6 +25,7 @@ set(fppFiles stdlib_linalg_outer_product.fypp stdlib_linalg_kronecker.fypp stdlib_linalg_cross_product.fypp + stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp stdlib_linalg_state.fypp stdlib_optval.fypp diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 460aa0593..134bf7af5 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -22,6 +22,8 @@ module stdlib_linalg public :: eye public :: lstsq public :: lstsq_space + public :: solve + public :: solve_lu public :: solve_lstsq public :: trace public :: outer_product @@ -228,6 +230,100 @@ module stdlib_linalg #:endfor end interface is_hessenberg + ! Solve linear system system Ax=b. + interface solve + !! version: experimental + !! + !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#solve-solves-a-linear-matrix-equation-or-a-linear-system-of-equations)) + !! + !!### Summary + !! Interface for solving a linear system arising from a general matrix. + !! + !!### Description + !! + !! This interface provides methods for computing the solution of a linear matrix system. + !! Supported data types include `real` and `complex`. No assumption is made on the matrix + !! structure. + !! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`) + !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems. + !! + !!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module function stdlib_linalg_${ri}$_solve_${ndsuf}$(a,b,overwrite_a,err) result(x) + !> Input matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out) :: err + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, allocatable, target :: x${nd}$ + end function stdlib_linalg_${ri}$_solve_${ndsuf}$ + pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x) + !> Input matrix a[n,n] + ${rt}$, intent(in) :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, allocatable, target :: x${nd}$ + end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface solve + + ! Solve linear system Ax = b using LU decomposition (subroutine interface). + interface solve_lu + !! version: experimental + !! + !! Solves the linear system \( A \cdot x = b \) for the unknown vector \( x \) from a square matrix \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#solve-lu-solves-a-linear-matrix-equation-or-a-linear-system-of-equations-subroutine-interface)) + !! + !!### Summary + !! Subroutine interface for solving a linear system using LU decomposition. + !! + !!### Description + !! + !! This interface provides methods for computing the solution of a linear matrix system using + !! a subroutine. Supported data types include `real` and `complex`. No assumption is made on the matrix + !! structure. Preallocated space for the solution vector `x` is user-provided, and it may be provided + !! for the array of pivot indices, `pivot`. If all pre-allocated work spaces are provided, no internal + !! memory allocations take place when using this interface. + !! The function can solve simultaneously either one (from a 1-d right-hand-side vector `b(:)`) + !! or several (from a 2-d right-hand-side vector `b(:,:)`) systems. + !! + !!@note The solution is based on LAPACK's generic LU decomposition based solvers `*GESV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err) + !> Input matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] Storage array for the diagonal pivot indices + integer(ilp), optional, intent(inout), target :: pivot(:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$ + #:endif + #:endfor + #:endfor + end interface solve_lu + ! Least squares solution to system Ax=b, i.e. such that the 2-norm abs(b-Ax) is minimized. interface lstsq !! version: experimental diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp new file mode 100644 index 000000000..8612d5052 --- /dev/null +++ b/src/stdlib_linalg_solve.fypp @@ -0,0 +1,170 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set RHS_SUFFIX = ["one","many"] +#:set RHS_SYMBOL = [ranksuffix(r) for r in [1,2]] +#:set RHS_EMPTY = [emptyranksuffix(r) for r in [1,2]] +#:set ALL_RHS = list(zip(RHS_SYMBOL,RHS_SUFFIX,RHS_EMPTY)) +submodule (stdlib_linalg) stdlib_linalg_solve +!! Solve linear system Ax=b + use stdlib_linalg_constants + use stdlib_linalg_lapack, only: gesv + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + implicit none(type,external) + + character(*), parameter :: this = 'solve' + + contains + + elemental subroutine handle_gesv_info(info,lda,n,nrhs,err) + integer(ilp), intent(in) :: info,lda,n,nrhs + type(linalg_state_type), intent(out) :: err + + ! Process output + select case (info) + case (0) + ! Success + case (-1) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid problem size n=',n) + case (-2) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid rhs size n=',nrhs) + case (-4) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size a=',[lda,n]) + case (-7) + err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n]) + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'singular matrix') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select + + end subroutine handle_gesv_info + + #:for nd,ndsuf,nde in ALL_RHS + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + ! Compute the solution to a real system of linear equations A * X = B + module function stdlib_linalg_${ri}$_solve_${ndsuf}$(a,b,overwrite_a,err) result(x) + !> Input matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), intent(out) :: err + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, allocatable, target :: x${nd}$ + + ! Initialize solution shape from the rhs array + allocate(x,mold=b) + + call stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,overwrite_a=overwrite_a,err=err) + + end function stdlib_linalg_${ri}$_solve_${ndsuf}$ + + !> Compute the solution to a real system of linear equations A * X = B (pure interface) + pure module function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$(a,b) result(x) + !> Input matrix a[n,n] + ${rt}$, intent(in) :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, allocatable, target :: x${nd}$ + + ! Local variables + ${rt}$, allocatable :: amat(:,:) + + ! Copy `a` so it can be intent(in) + allocate(amat,source=a) + + ! Initialize solution shape from the rhs array + allocate(x,mold=b) + + call stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(amat,b,x,overwrite_a=.true.) + + end function stdlib_linalg_${ri}$_pure_solve_${ndsuf}$ + + !> Compute the solution to a real system of linear equations A * X = B (pure interface) + pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err) + !> Input matrix a[n,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Right hand side vector or array, b[n] or b[n,nrhs] + ${rt}$, intent(in) :: b${nd}$ + !> Result array/matrix x[n] or x[n,nrhs] + ${rt}$, intent(inout), contiguous, target :: x${nd}$ + !> [optional] Storage array for the diagonal pivot indices + integer(ilp), optional, intent(inout), target :: pivot(:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + ! Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: lda,n,ldb,ldx,nrhsx,nrhs,info,npiv + integer(ilp), pointer :: ipiv(:) + logical(lk) :: copy_a + ${rt}$, pointer :: xmat(:,:),amat(:,:) + + ! Problem sizes + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + ldb = size(b,1,kind=ilp) + nrhs = size(b ,kind=ilp)/ldb + ldx = size(x,1,kind=ilp) + nrhsx = size(x ,kind=ilp)/ldx + + ! Has a pre-allocated pivots storage array been provided? + if (present(pivot)) then + ipiv => pivot + else + allocate(ipiv(n)) + endif + npiv = size(ipiv,kind=ilp) + + ! Can A be overwritten? By default, do not overwrite + if (present(overwrite_a)) then + copy_a = .not.overwrite_a + else + copy_a = .true._lk + endif + + if (any([lda,n,ldb]<1) .or. any([lda,ldb,ldx]/=n) .or. nrhsx/=nrhs .or. npiv/=n) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid sizes: a=',[lda,n], & + 'b=',[ldb,nrhs],' x=',[ldx,nrhsx], & + 'pivot=',n) + call linalg_error_handling(err0,err) + return + end if + + ! Initialize a matrix temporary + if (copy_a) then + allocate(amat(lda,n),source=a) + else + amat => a + endif + + ! Initialize solution with the rhs + x = b + xmat(1:n,1:nrhs) => x + + ! Solve system + call gesv(n,nrhs,amat,lda,ipiv,xmat,ldb,info) + + ! Process output + call handle_gesv_info(info,lda,n,nrhs,err0) + + if (copy_a) deallocate(amat) + if (.not.present(pivot)) deallocate(ipiv) + + ! Process output and return + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$ + + #:endif + #:endfor + #:endfor + +end submodule stdlib_linalg_solve diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 5a07b5bc4..9ec242213 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -2,6 +2,7 @@ set( fppFiles "test_linalg.fypp" "test_blas_lapack.fypp" + "test_linalg_solve.fypp" "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" "test_linalg_matrix_property_checks.fypp" @@ -11,5 +12,6 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_determinant) ADDTEST(linalg_matrix_property_checks) +ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(blas_lapack) diff --git a/test/linalg/test_linalg_lstsq.fypp b/test/linalg/test_linalg_lstsq.fypp index 29cd07519..9aba94871 100644 --- a/test/linalg/test_linalg_lstsq.fypp +++ b/test/linalg/test_linalg_lstsq.fypp @@ -70,14 +70,17 @@ module test_linalg_least_squares type(linalg_state_type) :: state integer(ilp), parameter :: n = 12, m = 3 + real :: Arnd(n,m),xrnd(m) ${rt}$ :: xsol(m),x(m),y(n),A(n,m) ! Random coefficient matrix and solution - call random_number(A) - call random_number(xsol) + call random_number(Arnd) + call random_number(xrnd) ! Compute rhs - y = matmul(A,xsol) + A = real(Arnd,${rk}$) + xsol = real(xrnd,${rk}$) + y = matmul(A,xsol) ! Find polynomial x = lstsq(A,y,err=state) diff --git a/test/linalg/test_linalg_solve.fypp b/test/linalg/test_linalg_solve.fypp new file mode 100644 index 000000000..25234a40b --- /dev/null +++ b/test/linalg/test_linalg_solve.fypp @@ -0,0 +1,190 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test linear system solver +module test_linalg_solve + use stdlib_linalg_constants + use stdlib_linalg_state + use stdlib_linalg, only: solve + use testdrive, only: error_type, check, new_unittest, unittest_type + + implicit none (type,external) + private + + public :: test_linear_systems + + contains + + !> Solve real and complex linear systems + subroutine test_linear_systems(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + tests = [tests,new_unittest("solve_${ri}$",test_${ri}$_solve), & + new_unittest("solve_${ri}$_multiple",test_${ri}$_solve_multiple)] + #:endif + #:endfor + + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + tests = [tests,new_unittest("solve_complex_${ci}$",test_${ci}$_solve), & + new_unittest("solve_2x2_complex_${ci}$",test_2x2_${ci}$_solve)] + #:endif + #:endfor + + end subroutine test_linear_systems + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + !> Simple linear system + subroutine test_${ri}$_solve(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$ :: A(3,3) = transpose(reshape([${rt}$ :: 1, 3, 3, & + 1, 3, 4, & + 1, 4, 3], [3,3])) + ${rt}$ :: b (3) = [${rt}$ :: 1, 4, -1] + ${rt}$ :: res(3) = [${rt}$ :: -2, -2, 3] + ${rt}$ :: x(3) + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) Simple linear system with multiple right hand sides + subroutine test_${ri}$_solve_multiple(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$ :: A(3,3) = transpose(reshape([${rt}$ :: 1,-1, 2, & + 0, 1, 1, & + 1,-1, 3], [3,3])) + ${rt}$ :: b(3,3) = transpose(reshape([${rt}$ :: 0, 1, 2, & + 1,-2,-1, & + 2, 3,-1], [3,3])) + ${rt}$ :: res(3,3) = transpose(reshape([${rt}$ ::-5,-7,10, & + -1,-4, 2, & + 2, 2,-3], [3,3])) + ${rt}$ :: x(3,3) + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) Complex linear system + !> Militaru, Popa, "On the numerical solving of complex linear systems", + !> Int J Pure Appl Math 76(1), 113-122, 2012. + subroutine test_${ri}$_solve(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$ :: A(5,5), b(5), res(5), x(5) + integer(ilp) :: i + + ! Fill in linear system + A = (0.0_${rk}$,0.0_${rk}$) + + A(1:2,1) = [(19.73_${rk}$,0.0_${rk}$),(0.0_${rk}$,-0.51_${rk}$)] + A(1:3,2) = [(12.11_${rk}$,-1.0_${rk}$),(32.3_${rk}$,7.0_${rk}$),(0.0_${rk}$,-0.51_${rk}$)] + A(1:4,3) = [(0.0_${rk}$,5.0_${rk}$),(23.07_${rk}$,0.0_${rk}$),(70.0_${rk}$,7.3_${rk}$),(1.0_${rk}$,1.1_${rk}$)] + A(2:5,4) = [(0.0_${rk}$,1.0_${rk}$),(3.95_${rk}$,0.0_${rk}$),(50.17_${rk}$,0.0_${rk}$),(0.0_${rk}$,-9.351_${rk}$)] + A(3:5,5) = [(19.0_${rk}$,31.83_${rk}$),(45.51_${rk}$,0.0_${rk}$),(55.0_${rk}$,0.0_${rk}$)] + + b = [(77.38_${rk}$,8.82_${rk}$),(157.48_${rk}$,19.8_${rk}$),(1175.62_${rk}$,20.69_${rk}$),(912.12_${rk}$,-801.75_${rk}$),(550.0_${rk}$,-1060.4_${rk}$)] + + ! Exact result + res = [(3.3_${rk}$,-1.0_${rk}$),(1.0_${rk}$,0.17_${rk}$),(5.5_${rk}$,0.0_${rk}$),(9.0_${rk}$,0.0_${rk}$),(10.0_${rk}$,-17.75_${rk}$)] + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) 2x2 Complex linear system + subroutine test_2x2_${ri}$_solve(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + ${rt}$, parameter :: i = (0.0_${rk}$,1.0_${rk}$) + + ${rt}$ :: A(2,2), b(2), res(2), x(2) + + ! Fill in linear system + A(1,:) = [ 1+2*i, 2-i] + A(2,:) = [ 2+i , i] + + b = [1,-1] + + ! Exact result + res = [(-0.28_${rk}$,-0.04_${rk}$),(0.36_${rk}$,0.48_${rk}$)] + + x = solve(a,b,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(x-res) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_solve +