Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

linalg: determinant #798

Merged
merged 29 commits into from
Apr 25, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
4a4b899
import files
perazz Apr 13, 2024
d54cfa3
add `pure` interface
perazz Apr 13, 2024
c987fa1
add operator interface
perazz Apr 13, 2024
b9a3b0e
add tests
perazz Apr 13, 2024
4a29219
make `pure`
perazz Apr 13, 2024
3d05cc3
remove `xdp`
perazz Apr 13, 2024
4beab1f
Documentation
perazz Apr 13, 2024
15023e1
`submodule version`
perazz Apr 14, 2024
5923a54
Update src/stdlib_linalg.fypp
perazz Apr 15, 2024
31cdc84
Update src/stdlib_linalg.fypp
perazz Apr 15, 2024
45a606f
add `det` example
perazz Apr 15, 2024
cacb585
Update src/stdlib_linalg_determinant.fypp
perazz Apr 15, 2024
5c16ff8
Update src/stdlib_linalg_determinant.fypp
perazz Apr 15, 2024
ab030c5
Update src/stdlib_linalg_determinant.fypp
perazz Apr 15, 2024
13bd98a
warn about `xdp`
perazz Apr 15, 2024
7bf7141
Merge branch 'determinant' of github.com:perazz/stdlib into determinant
perazz Apr 15, 2024
504d90d
cleanup xdp notes
perazz Apr 15, 2024
e80b508
spacing
perazz Apr 15, 2024
c6076ea
relax error thresholds
perazz Apr 15, 2024
eaebe5c
restore `pure` attribute
perazz Apr 15, 2024
5d52d48
Update test/linalg/test_linalg_determinant.fypp
perazz Apr 15, 2024
cff995d
add docs
perazz Apr 21, 2024
2fe2428
link to specs
perazz Apr 21, 2024
45447a8
Update doc/specs/stdlib_linalg.md
perazz Apr 21, 2024
3ee20ad
Update doc/specs/stdlib_linalg.md
perazz Apr 21, 2024
1e50115
Update doc/specs/stdlib_linalg.md
perazz Apr 21, 2024
fe933ad
Update doc/specs/stdlib_linalg.md
perazz Apr 21, 2024
59dae89
Update doc/specs/stdlib_linalg.md
perazz Apr 21, 2024
3c72b06
Update doc/specs/stdlib_linalg.md
perazz Apr 21, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
74 changes: 74 additions & 0 deletions doc/specs/stdlib_linalg.md
Original file line number Diff line number Diff line change
Expand Up @@ -595,3 +595,77 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l
```fortran
{!example/linalg/example_is_hessenberg.f90!}
```

## `det` - Computes the determinant of a square matrix

### Status

Experimental

### Description

This function computes the determinant of a `real` or `complex` square matrix.

This interface comes with a `pure` version `det(a)`, and a non-pure version `det(a,overwrite_a,err)` that
allows for more expert control.

### Syntax

`c = ` [[stdlib_linalg(module):det(interface)]] `(a [, overwrite_a, err])`

### Arguments

`a`: Shall be a rank-2 square array

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

### Return value

Returns a `real` scalar value of the same kind of `a` that represents the determinant of the matrix.

Raises `LINALG_ERROR` if the matrix is singular.
Raises `LINALG_VALUE_ERROR` if the matrix is non-square.
Exceptions are returned to the `err` argument if provided; an `error stop` is triggered otherwise.

### Example

```fortran
{!example/linalg/example_determinant.f90!}
```

## `.det.` - Determinant operator of a square matrix

### Status

Experimental

### Description

This operator returns the determinant of a real square matrix.

This interface is equivalent to the `pure` version of determinant [[stdlib_linalg(module):det(interface)]].

### Syntax

`c = ` [[stdlib_linalg(module):operator(.det.)(interface)]] `(a)`

### Arguments

`a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument.

### Return value

Returns a real scalar value that represents the determinnt of the matrix.

Raises `LINALG_ERROR` if the matrix is singular.
Raises `LINALG_VALUE_ERROR` if the matrix is non-square.
Exceptions trigger an `error stop`.

### Example

```fortran
{!example/linalg/example_determinant2.f90!}
```
3 changes: 3 additions & 0 deletions example/linalg/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -18,3 +18,6 @@ ADD_EXAMPLE(state1)
ADD_EXAMPLE(state2)
ADD_EXAMPLE(blas_gemv)
ADD_EXAMPLE(lapack_getrf)
ADD_EXAMPLE(determinant)
ADD_EXAMPLE(determinant2)

14 changes: 14 additions & 0 deletions example/linalg/example_determinant.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
program example_determinant
use stdlib_kinds, only: dp
use stdlib_linalg, only: det, linalg_state_type
implicit none
type(linalg_state_type) :: err

real(dp) :: d

! Compute determinate of a real matrix
d = det(reshape([real(dp)::1,2,3,4],[2,2]))

print *, d ! a*d-b*c = -2.0

end program example_determinant
13 changes: 13 additions & 0 deletions example/linalg/example_determinant2.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
program example_determinant2
use stdlib_kinds, only: dp
use stdlib_linalg, only: operator(.det.)
implicit none

real(dp) :: d

! Compute determinate of a real matrix
d = .det.reshape([real(dp)::1,2,3,4],[2,2])

print *, d ! a*d-b*c = -2.0

end program example_determinant2
1 change: 1 addition & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ set(fppFiles
stdlib_linalg_outer_product.fypp
stdlib_linalg_kronecker.fypp
stdlib_linalg_cross_product.fypp
stdlib_linalg_determinant.fypp
stdlib_linalg_state.fypp
stdlib_optval.fypp
stdlib_selection.fypp
Expand Down
114 changes: 112 additions & 2 deletions src/stdlib_linalg.fypp
Original file line number Diff line number Diff line change
@@ -1,15 +1,19 @@
#:include "common.fypp"
#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
#:set RCI_KINDS_TYPES = RC_KINDS_TYPES + INT_KINDS_TYPES
module stdlib_linalg
!!Provides a support for various linear algebra procedures
!! ([Specification](../page/specs/stdlib_linalg.html))
use stdlib_kinds, only: sp, dp, xdp, qp, &
use stdlib_kinds, only: sp, dp, xdp, qp, lk, &
int8, int16, int32, int64
use stdlib_error, only: error_stop
use stdlib_optval, only: optval
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling
implicit none
private

public :: det
public :: operator(.det.)
public :: diag
public :: eye
public :: trace
Expand All @@ -23,6 +27,9 @@ module stdlib_linalg
public :: is_hermitian
public :: is_triangular
public :: is_hessenberg

! Export linalg error handling
public :: linalg_state_type, linalg_error_handling
Comment on lines +30 to +32
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This let me think that we should review the structure of the specs. Currently we have 1 page per module. However, if these become available through stdlib_linalg, the specs should be modified accordingly (probably in another PR; I can start to work on a proposition when this is ready).


interface diag
!! version: experimental
Expand Down Expand Up @@ -214,6 +221,109 @@ module stdlib_linalg
#:endfor
end interface is_hessenberg


interface det
!! version: experimental
!!
!! Computes the determinant of a square matrix
!! ([Specification](../page/specs/stdlib_linalg.html#det-computes-the-determinant-of-a-square-matrix))
!!
!!### Summary
!! Interface for computing matrix determinant.
!!
!!### Description
!!
!! This interface provides methods for computing the determinant of a matrix.
!! Supported data types include `real` and `complex`.
!!
!!@note The provided functions are intended for square matrices only.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
!!### Example
!!
!!```fortran
!!
!! real(sp) :: a(3,3), d
!! type(linalg_state_type) :: state
!! a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
!!
!! ! ...
!! d = det(a,err=state)
!! if (state%ok()) then
!! print *, 'Success! det=',d
!! else
!! print *, state%print()
!! endif
!! ! ...
!!```
!!
#:for rk,rt in RC_KINDS_TYPES
#:if rk!="xdp"
module procedure stdlib_linalg_${rt[0]}$${rk}$determinant
module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant
#:endif
#:endfor
end interface det

interface operator(.det.)
!! version: experimental
!!
!! Determinant operator of a square matrix
!! ([Specification](../page/specs/stdlib_linalg.html#det-determinant-operator-of-a-square-matrix))
!!
!!### Summary
!! Pure operator interface for computing matrix determinant.
!!
!!### Description
!!
!! This pure operator interface provides a convenient way to compute the determinant of a matrix.
!! Supported data types include real and complex.
!!
!!@note The provided functions are intended for square matrices.
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
!!
!!### Example
!!
!!```fortran
!!
!! ! ...
!! real(sp) :: matrix(3,3), d
!! matrix = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
!! d = .det.matrix
!! ! ...
!!
!!```
!
#:for rk,rt in RC_KINDS_TYPES
#:if rk!="xdp"
module procedure stdlib_linalg_pure_${rt[0]}$${rk}$determinant
#:endif
#:endfor
end interface operator(.det.)

interface
#:for rk,rt in RC_KINDS_TYPES
#:if rk!="xdp"
module function stdlib_linalg_${rt[0]}$${rk}$determinant(a,overwrite_a,err) result(det)
!> Input matrix a[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> State return flag.
type(linalg_state_type), intent(out) :: err
!> Matrix determinant
${rt}$ :: det
end function stdlib_linalg_${rt[0]}$${rk}$determinant
pure module function stdlib_linalg_pure_${rt[0]}$${rk}$determinant(a) result(det)
!> Input matrix a[m,n]
${rt}$, intent(in) :: a(:,:)
!> Matrix determinant
${rt}$ :: det
end function stdlib_linalg_pure_${rt[0]}$${rk}$determinant
#:endif
#:endfor
end interface

contains


Expand Down
Loading
Loading