Merge branch 'master' into qr

Author: perazz

CHANGELOG.md

@@ -1,3 +1,43 @@
+# Version 0.7.0
+
+Full release notes available at [v0.7.0] tag.
+
+[v0.7.0]: https://github.com/fortran-lang/stdlib/releases/tag/v0.7.0
+
+- new module `stdlib_constants`
+  [#800](https://github.com/fortran-lang/stdlib/pull/800)
+  - Many mathematical constants and most common physical ([codata](https://codata.org)) constants
+
+Changes to existing scripts and modules
+  - changes in CI
+    - Use of `fortran-setup` for GCC, Intel LLVM and Intel Classic
+      [#834](https://github.com/fortran-lang/stdlib/pull/834)
+  - change in module `stdlib_hashmaps`
+    - Support of hash map key generic interfaces
+      [#827](https://github.com/fortran-lang/stdlib/pull/827)
+  - changes in module `stdlib_io`
+    - Addition of a Fortran format specifier in `loadtxt`
+      [#805](https://github.com/fortran-lang/stdlib/pull/805)
+  - changes in module `stdlib_linalg`
+    - Support of extended and quad precision checking 
+      [#821](https://github.com/fortran-lang/stdlib/pull/821)
+    - Several fixes
+      [#815](https://github.com/fortran-lang/stdlib/pull/815)
+      [#818](https://github.com/fortran-lang/stdlib/pull/818)
+      [#826](https://github.com/fortran-lang/stdlib/pull/826)
+      [#830](https://github.com/fortran-lang/stdlib/pull/830)
+      [#836](https://github.com/fortran-lang/stdlib/pull/836)
+    - New procedures for Eigenvalues and Eigenvectors computation: `eig`, `eigh`, `eigvals`, `eigvalsh`
+      [#816](https://github.com/fortran-lang/stdlib/pull/816)
+    - New procedures for Singular Value Decomposition: `svd`, `svdvals`
+      [#808](https://github.com/fortran-lang/stdlib/pull/808)
+  - changes in module `stdlib_sorting`
+    - Renamed variable from `int_size` to `int_index`
+      [#824](https://github.com/fortran-lang/stdlib/pull/824)
+    - Support of `int32` `index` array in `sort_index`
+      [#829](https://github.com/fortran-lang/stdlib/pull/829)
+
+
 # Version 0.6.1
 
 Full release notes available at [v0.6.1] tag.

VERSION

@@ -1 +1 @@
-0.6.1
+0.7.0

doc/specs/stdlib_linalg.md

@@ -687,6 +687,8 @@ Expert (`Pure`) interface:
 
 `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
 
 For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.
@@ -975,6 +977,179 @@ This subroutine computes the internal working space requirements for the QR fact
 {!example/linalg/example_qr_space.f90!}
 ```
 
+## `eig` - Eigenvalues and Eigenvectors of a Square Matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the solution to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank, `real` or `complex` matrix.
+
+Result array `lambda` returns the eigenvalues of \( A \). The user can request eigenvectors to be returned: if provided, on output `left` will contain the left eigenvectors, `right` the right eigenvectors of \( A \).
+Both `left` and `right` are rank-2 arrays, where eigenvectors are stored as columns.
+The solver is based on LAPACK's `*GEEV` backends.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):eig(interface)]] `(a, lambda [, right] [,left] [,overwrite_a] [,err])`
+
+### Arguments
+
+`a` : `real` or `complex` square array containing the coefficient matrix. If `overwrite_a=.false.`, it is an `intent(in)` argument. Otherwise, it is an `intent(inout)` argument and is destroyed by the call. 
+
+`lambda`: Shall be a `complex` or `real` rank-1 array of the same kind as `a`, containing the eigenvalues, or their `real` component only. It is an `intent(out)` argument.
+
+`right` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the right eigenvectors of `a`. It is an `intent(out)` argument.
+
+`left` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the left eigenvectors of `a`. It is an `intent(out)` 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.
+
+### Return value
+
+Raises `LINALG_ERROR` if the calculation did not converge.
+Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
+Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_eig.f90!}
+```
+
+## `eigh` - Eigenvalues and Eigenvectors of a Real symmetric or Complex Hermitian Square Matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the solution to the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \),
+where \( A \) is a square, full-rank, `real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix.
+
+Result array `lambda` returns the `real` eigenvalues of \( A \). The user can request the orthogonal eigenvectors 
+to be returned: on output `vectors` may contain the matrix of eigenvectors, returned as a column.
+
+Normally, only the lower triangular part of \( A \) is accessed. On input, `logical` flag `upper_a` 
+allows the user to request what triangular part of the matrix should be used.
+
+The solver is based on LAPACK's `*SYEV` and `*HEEV` backends.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):eigh(interface)]] `(a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])`
+
+### Arguments
+
+`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. 
+
+`lambda`: Shall be a `complex` rank-1 array of the same precision as `a`, containing the eigenvalues. It is an `intent(out)` argument. 
+
+`vectors` (optional): Shall be a rank-2 array of the same type, size and kind as `a`, containing the eigenvectors of `a`. It is an `intent(out)` argument.
+
+`upper_a` (optional): Shall be an input `logical` flag. If `.true.`, the upper triangular part of `a` will be accessed. Otherwise, the lower triangular part will be accessed. 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.
+
+### Return value
+
+Raises `LINALG_ERROR` if the calculation did not converge.
+Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_eigh.f90!}
+```
+
+## `eigvals` - Eigenvalues of a Square Matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This function returns the eigenvalues to matrix \( A \): a square, full-rank, `real` or `complex` matrix.
+The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \).
+
+Result array `lambda` is `complex`, and returns the eigenvalues of \( A \). 
+The solver is based on LAPACK's `*GEEV` backends.
+
+### Syntax
+
+`lambda = ` [[stdlib_linalg(module):eigvals(interface)]] `(a, [,err])`
+
+### Arguments
+
+`a` : `real` or `complex` square array containing the coefficient matrix. It 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 `complex` array containing the eigenvalues of `a`. 
+
+Raises `LINALG_ERROR` if the calculation did not converge.
+Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_eigvals.f90!}
+```
+
+## `eigvalsh` - Eigenvalues of a Real Symmetric or Complex Hermitian Square Matrix
+
+### Status
+
+Experimental
+
+### Description
+
+This function returns the eigenvalues to matrix \( A \): a where \( A \) is a square, full-rank, 
+`real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix.
+The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \).
+
+Result array `lambda` is `real`, and returns the eigenvalues of \( A \). 
+The solver is based on LAPACK's `*SYEV` and `*HEEV` backends.
+
+### Syntax
+
+`lambda = ` [[stdlib_linalg(module):eigvalsh(interface)]] `(a, [, upper_a] [,err])`
+
+### Arguments
+
+`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument.
+
+`upper_a` (optional): Shall be an input logical flag. If `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed (default). It 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` array containing the eigenvalues of `a`. 
+
+Raises `LINALG_ERROR` if the calculation did not converge.
+Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_eigvalsh.f90!}
+```
+
 ## `svd` - Compute the singular value decomposition of a rank-2 array (matrix).
 
 ### Status
@@ -1066,4 +1241,3 @@ Exceptions trigger an `error stop`, unless argument `err` is present.
 ```fortran
 {!example/linalg/example_svdvals.f90!}
 ```
-

example/linalg/CMakeLists.txt

@@ -13,6 +13,10 @@ ADD_EXAMPLE(is_square)
 ADD_EXAMPLE(is_symmetric)
 ADD_EXAMPLE(is_triangular)
 ADD_EXAMPLE(outer_product)
+ADD_EXAMPLE(eig)
+ADD_EXAMPLE(eigh)
+ADD_EXAMPLE(eigvals)
+ADD_EXAMPLE(eigvalsh)
 ADD_EXAMPLE(trace)
 ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)

example/linalg/example_eig.f90

@@ -0,0 +1,26 @@
+! Eigendecomposition of a real square matrix 
+program example_eig
+  use stdlib_linalg, only: eig
+  implicit none
+
+  integer :: i
+  real, allocatable :: A(:,:)
+  complex, allocatable :: lambda(:),vectors(:,:)
+
+  ! Decomposition of this square matrix
+  ! NB Fortran is column-major -> transpose input
+  A = transpose(reshape( [ [2, 2, 4], &
+                           [1, 3, 5], &
+                           [2, 3, 4] ], [3,3] )) 
+
+  ! Get eigenvalues and right eigenvectors
+  allocate(lambda(3),vectors(3,3))
+  
+  call eig(A, lambda, right=vectors)
+  
+  do i=1,3
+     print *, 'eigenvalue  ',i,': ',lambda(i)
+     print *, 'eigenvector ',i,': ',vectors(:,i)
+  end do
+  
+end program example_eig

example/linalg/example_eigh.f90

@@ -0,0 +1,38 @@
+! Eigendecomposition of a real symmetric matrix 
+program example_eigh
+  use stdlib_linalg, only: eigh
+  implicit none
+
+  integer :: i
+  real, allocatable :: A(:,:),lambda(:),vectors(:,:)
+  complex, allocatable :: cA(:,:),cvectors(:,:)
+
+  ! Decomposition of this symmetric matrix
+  ! NB Fortran is column-major -> transpose input
+  A = transpose(reshape( [ [2, 1, 4], &
+                           [1, 3, 5], &
+                           [4, 5, 4] ], [3,3] )) 
+
+  ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors
+  allocate(lambda(3),vectors(3,3))  
+  call eigh(A, lambda, vectors=vectors)
+  
+  print *, 'Real matrix'
+  do i=1,3
+     print *, 'eigenvalue  ',i,': ',lambda(i)
+     print *, 'eigenvector ',i,': ',vectors(:,i)
+  end do
+  
+  ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors
+  cA = A
+  
+  allocate(cvectors(3,3))  
+  call eigh(cA, lambda, vectors=cvectors)
+  
+  print *, 'Complex matrix'
+  do i=1,3
+     print *, 'eigenvalue  ',i,': ',lambda(i)
+     print *, 'eigenvector ',i,': ',cvectors(:,i)
+  end do  
+  
+end program example_eigh

example/linalg/example_eigvals.f90

@@ -0,0 +1,24 @@
+! Eigenvalues of a general real / complex matrix 
+program example_eigvals
+  use stdlib_linalg, only: eigvals
+  implicit none
+
+  integer :: i
+  real, allocatable :: A(:,:),lambda(:)
+  complex, allocatable :: cA(:,:),clambda(:)
+
+  ! NB Fortran is column-major -> transpose input
+  A = transpose(reshape( [ [2, 8, 4], &
+                           [1, 3, 5], &
+                           [9, 5,-2] ], [3,3] )) 
+
+  ! Note: real symmetric matrix
+  lambda = eigvals(A)
+  print *, 'Real    matrix eigenvalues: ',lambda
+  
+  ! Complex general matrix
+  cA = cmplx(A, -2*A)
+  clambda = eigvals(cA)
+  print *, 'Complex matrix eigenvalues: ',clambda
+  
+end program example_eigvals

example/linalg/example_eigvalsh.f90

@@ -0,0 +1,25 @@
+! Eigenvalues of a real symmetric / complex hermitian matrix 
+program example_eigvalsh
+  use stdlib_linalg, only: eigvalsh
+  implicit none
+
+  integer :: i
+  real, allocatable :: A(:,:),lambda(:)
+  complex, allocatable :: cA(:,:)
+
+  ! Decomposition of this symmetric matrix
+  ! NB Fortran is column-major -> transpose input
+  A = transpose(reshape( [ [2, 1, 4], &
+                           [1, 3, 5], &
+                           [4, 5, 4] ], [3,3] )) 
+
+  ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors
+  lambda = eigvalsh(A)
+  print *, 'Symmetric matrix eigenvalues: ',lambda
+  
+  ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors
+  cA = A
+  lambda = eigvalsh(cA)
+  print *, 'Hermitian matrix eigenvalues: ',lambda
+  
+end program example_eigvalsh

include/common.fypp

@@ -60,6 +60,19 @@
 #:set CMPLX_INIT = CMPLX_INIT + ["w"]
 #:endif
 
+#! BLAS/LAPACK complex->real kind initial conversion
+#! Converts a BLAS/LAPACK complex kind initial to a real kind initial
+#!
+#! Args:
+#!     ci (character): Complex kind initial in ["c","z","y","w"]
+#!
+#! Returns:
+#!     Real kind initial in ["s","d","x","q"] or an empty string on invalid input
+#!
+#:def c2ri(cmplx)
+$:"s" if cmplx=="c" else "d" if cmplx=="z" else "x" if cmplx=="y" else "q" if cmplx=="w" else "ERROR"
+#:enddef
+
 #! Complex types to be considered during templating
 #:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS]
 

src/CMakeLists.txt

@@ -27,7 +27,8 @@ set(fppFiles
     stdlib_linalg_outer_product.fypp
     stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp
-    stdlib_linalg_solve.fypp
+    stdlib_linalg_eigenvalues.fypp
+    stdlib_linalg_solve.fypp    
     stdlib_linalg_determinant.fypp
     stdlib_linalg_qr.fypp
     stdlib_linalg_state.fypp