线性代数

线性代数

stdlib 线性代数库提供了用于处理常见线性代数运算的高级 API。

BLAS 和 LAPACK

状态

实验性的

描述

BLAS 和 LAPACK 后端为许多线性代数算法提供了高效的低级实现，并被用于非平凡的运算符。一个现代 Fortran 版本的 Reference-LAPACK 3.10.1 实现被提供为后端。经过自动转换遗留代码后，提供了具有完整显式类型特征的现代 Fortran 模块： - [stdlib_linalg_blas(module)]、[stdlib_linalg_lapack(module)] 为所有函数提供与 KIND 无关的接口。 - 这两个库都可用于 32 位 (sp)、64 位 (dp) 和 128 位 (qp) real 和 complex 数（如果在当前构建中可用）。 - 自由格式，小写样式 - implicit none(type, external) 应用于所有过程和模块 - intent 已添加，并且所有 pure 过程尽可能使用 - stdlib 以两种不同的风格提供所有过程：（a）原始 BLAS/LAPACK 名称，前缀为 stdlib_?<name>（例如：stdlib_dgemv、stdlib_sgemv）；（b）泛型，与 KIND 无关的 <name>，即 gemv。 - F77 样式的 parameter 已移除，并且所有数字常量都已使用依赖于 KIND 的 Fortran 内在函数进行泛化。 - 保留了基于预处理器的 OpenMP 指令。单源模块结构有望允许跨过程内联，否则没有链接时优化将无法实现。

如果可用，应优先使用利用专用处理器指令的高度优化库，而不是 stdlib 实现。此类库的示例包括：OpenBLAS、MKL (TM)、Accelerate 和 ATLAS。为了启用它们的用法，只需确保定义了以下预处理器宏

STDLIB_EXTERNAL_BLAS 将所有 BLAS 过程（除了 128 位过程）封装到外部库
STDLIB_EXTERNAL_LAPACK 将所有 LAPACK 过程（除了 128 位过程）封装到外部库

这些可以在构建过程中启用。例如，使用 CMake，可以使用 add_compile_definitions(STDLIB_EXTERNAL_BLAS STDLIB_EXTERNAL_LAPACK) 启用这些预处理器指令。从 fpm 分支也是可能的，其中默认情况下启用了 cpp 预处理器。例如，宏可以添加到项目的清单中

# Link against appropriate external BLAS and LAPACK libraries, if necessary
[build]
link = ["blas", "lapack"]  

[dependencies]
stdlib="*"

# Macros are only needed if using an external library
[preprocess]
[preprocess.cpp]
macros = ["STDLIB_EXTERNAL_BLAS", "STDLIB_EXTERNAL_LAPACK"]

或直接通过编译器标志

fpm build --flag "-DSTDLIB_EXTERNAL_BLAS -DSTDLIB_EXTERNAL_LAPACK -lblas -llapack".

语法

BLAS 和 LAPACK 后端中的所有过程都遵循来自 Reference LAPACK 的标准接口。因此，应参考在线用户指南以了解完整 API 和过程参数及其用法的描述。

stdlib 实现通过模块 [stdlib_linalg_blas(module)] 和 [stdlib_linalg_lapack(module)] 提供与 KIND 无关和特定过程接口。因为所有过程都以一个字母开头表示基本数据类型，所以 stdlib 泛型接口省略了开头字母，并包含所有依赖于 KIND 的实现。例如，axpy 函数的泛型接口如下所示

!> AXPY: constant times a vector plus a vector.
interface axpy
    module procedure stdlib_saxpy
    module procedure stdlib_daxpy
    module procedure stdlib_qaxpy
    module procedure stdlib_caxpy
    module procedure stdlib_zaxpy
    module procedure stdlib_waxpy
end interface axpy

泛型接口是使用外部库的终点。只要使用外部库，对内部 module procedure 的引用就会被替换为对外部库的接口，例如

!> AXPY: constant times a vector plus a vector.
interface axpy
    pure subroutine caxpy(n,ca,cx,incx,cy,incy)
        import sp,dp,qp,ilp,lk 
        implicit none(type,external) 
        complex(sp), intent(in) :: ca,cx(*)
        integer(ilp), intent(in) :: incx,incy,n
        complex(sp), intent(inout) :: cy(*)
    end subroutine caxpy
    ! [....]
    module procedure stdlib_qaxpy
end interface axpy

请注意，128 位函数仅由 stdlib 提供，并且始终指向内部实现。由于 128 位精度被识别为 [stdlib_kinds(module):qp]，因此 128 位过程的开头字母被标记为 q（四精度实数）和 w（“宽”或四精度复数）。扩展精度 ([stdlib_kinds(module):xdp]) 计算被标记为 x（扩展精度实数）。和 y（扩展精度复数）。

示例

program example_gemv
  use stdlib_linalg, only: eye
  use stdlib_linalg_blas, only: sp,gemv
  implicit none(type,external)
  real(sp) :: A(2, 2), B(2), C(2)
  B = [1.0,2.0]
  A = eye(2)

  ! Use legacy BLAS interface 
  call gemv('No transpose',m=size(A,1),n=size(A,2),alpha=1.0,a=A,lda=size(A,1),x=B,incx=1,beta=0.0,y=C,incy=1)

  print *, C ! returns 1.0 2.0

end program example_gemv

program example_getrf
  use stdlib_linalg, only: eye
  use stdlib_linalg_lapack, only: dp,ilp,getrf
  implicit none(type,external)
  real(dp) :: A(3, 3)
  integer(ilp) :: ipiv(3),info

  A = eye(3)

  ! LAPACK matrix factorization interface (overwrite result)
  call getrf(size(A,1),size(A,2),A,size(A,1),ipiv,info)
  print *, info ! info==0: Success!

end program example_getrf

许可

Fortran 标准库是在 MIT 许可证下分发的。LAPACK 及其包含的 BLAS 是一个免费提供的软件包。它们可通过匿名 ftp 和万维网从 netlib 获取。因此，它们可以包含在商业软件包中（并且已经包含）。用于 BLAS 和 LAPACK 后端的许可证是修改后的 BSD 许可证。

LICENSE.txt 文件的标题将其许可要求写为

Copyright (c) 1992-2013 The University of Tennessee and The University
                        of Tennessee Research Foundation.  All rights
                        reserved.
Copyright (c) 2000-2013 The University of California Berkeley. All
                        rights reserved.
Copyright (c) 2006-2013 The University of Colorado Denver.  All rights
                        reserved.

$COPYRIGHT$

Additional copyrights may follow

$HEADER$

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

- Redistributions of source code must retain the above copyright
  notice, this list of conditions and the following disclaimer.

- Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions and the following disclaimer listed
  in this license in the documentation and/or other materials
  provided with the distribution.

- Neither the name of the copyright holders nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

The copyright holders provide no reassurances that the source code
provided does not infringe any patent, copyright, or any other
intellectual property rights of third parties.  The copyright holders
disclaim any liability to any recipient for claims brought against
recipient by any third party for infringement of that parties
intellectual property rights.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

因此，LICENSE.txt 代码的许可证与在 MIT 许可证下在 Fortran 标准库中使用修改后的代码版本兼容。
应将 BLAS、LAPACK 库的功劳归于 LAPACK 作者。根据原始许可证，我们还更改了例程的名称，并对原始代码所做的更改进行了注释。

`diag` - 创建对角线数组或提取数组的对角线元素

状态

实验性的

描述

创建对角线数组或提取数组的对角线元素

语法

d = diag (a [, k])

参数

a: 应为秩 1 或秩 2 数组。如果 a 是秩 1 数组（即向量），则 diag 返回一个秩 2 数组，其元素位于对角线上。如果 a 是秩 2 数组（即矩阵），则 diag 返回对角线元素的秩 1 数组。

k (可选): 应为类型 integer 的标量，并指定对角线。默认的 k = 0 表示主对角线，k > 0 是主对角线以上对角线，k < 0 是主对角线以下对角线。

返回值

返回一个对角线数组或一个包含提取的对角线元素的向量。

示例

program example_diag1
  use stdlib_linalg, only: diag
  implicit none
  real, allocatable :: A(:, :)
  integer :: i
  A = diag([(1, i=1, 10)]) ! creates a 10 by 10 identity matrix
end program example_diag1

program example_diag2
  use stdlib_linalg, only: diag
  implicit none
  real, allocatable :: v(:)
  real, allocatable :: A(:, :)
  v = [1, 2, 3, 4, 5]
  A = diag(v) ! creates a 5 by 5 matrix with elements of v on the diagonal
end program example_diag2

program example_diag3
  use stdlib_linalg, only: diag
  implicit none
  integer, parameter :: n = 10
  real :: c(n), ul(n - 1)
  real :: A(n, n)
  c = 2
  ul = -1
  A = diag(ul, -1) + diag(c) + diag(ul, 1) ! Gil Strang's favorite matrix
end program example_diag3

program example_diag4
  use stdlib_linalg, only: diag
  implicit none
  integer, parameter :: n = 12
  real :: A(n, n)
  real :: v(n)
  call random_number(A)
  v = diag(A) ! v contains diagonal elements of A
end program example_diag4

program example_diag5
  use stdlib_linalg, only: diag
  implicit none
  integer, parameter :: n = 3
  real :: A(n, n)
  real, allocatable :: v(:)
  A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [n, n])
  v = diag(A, -1) ! v is [2,6]
  v = diag(A, 1)  ! v is [4,8]
end program example_diag5

`eye` - 构造单位矩阵

状态

实验性的

类

纯函数。

描述

构造单位矩阵。

语法

I = eye (dim1 [, dim2])

参数

dim1: 应为默认类型 integer 的标量。这是一个 intent(in) 参数。

dim2: 应为默认类型 integer 的标量。这是一个 intent(in) 和 optional 参数。

返回值

返回单位矩阵，即主对角线上为 1，其他位置为 0 的矩阵。返回值的类型为 integer(int8)。建议使用 int8 来节省存储空间。

警告

由于 eye 的结果类型为 integer(int8)，因此在算术表达式中使用它时应小心。例如

!> Be careful
A = eye(2,2)/2     !! A == 0.0
!> Recommend
A = eye(2,2)/2.0   !! A == diag([0.5, 0.5])

示例

program example_eye1
  use stdlib_linalg, only: eye
  implicit none
  integer :: i(2, 2)
  real :: a(3, 3)
  real :: b(2, 3)  !! Matrix is non-square.
  complex :: c(2, 2)
  I = eye(2)              !! [1,0; 0,1]
  A = eye(3)              !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
  A = eye(3, 3)            !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
  B = eye(2, 3)            !! [1.0,0.0,0.0; 0.0,1.0,0.0]
  C = eye(2, 2)            !! [(1.0,0.0),(0.0,0.0); (0.0,0.0),(1.0,0.0)]
  C = (1.0, 1.0)*eye(2, 2)  !! [(1.0,1.0),(0.0,0.0); (0.0,0.0),(1.0,1.0)]
end program example_eye1

program example_eye2
  use stdlib_linalg, only: eye, diag
  implicit none
  print *, all(eye(4) == diag([1, 1, 1, 1])) ! prints .true.
end program example_eye2

`trace` - 矩阵的迹

状态

实验性的

描述

矩阵（秩 2 数组）的迹

语法

result = trace (A)

参数

A: 应为秩 2 数组。如果 A 不是方阵，则 trace(A) 将返回 A 的方阵子部分的对角线值的总和。

返回值

返回矩阵的迹，即对角线元素的总和。

示例

program example_trace
  use stdlib_linalg, only: trace
  implicit none
  real :: A(3, 3)
  A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
  print *, trace(A) ! 1 + 5 + 9
end program example_trace

`outer_product` - 计算两个向量的外积

状态

实验性的

描述

计算两个向量的外积

语法

d = outer_product (u, v)

参数

u: 应为秩 1 数组

v: 应为秩 1 数组

返回值

返回一个秩 2 数组，等于 u v^T（其中 u, v 被视为列向量）。返回数组的形状为 [size(u), size(v)]。

示例

program example_outer_product
  use stdlib_linalg, only: outer_product
  implicit none
  real, allocatable :: A(:, :), u(:), v(:)
  u = [1., 2., 3.]
  v = [3., 4.]
  A = outer_product(u, v)
!A = reshape([3., 6., 9., 4., 8., 12.], [3,2])
end program example_outer_product

`kronecker_product` - 计算两个秩 2 数组的克罗内克积

状态

实验性的

描述

计算两个秩 2 数组的克罗内克积

语法

C = kronecker_product (A, B)

参数

A: 应为维度为 M1、N1 的秩 2 数组

B: 应为维度为 M2、N2 的秩 2 数组

返回值

返回一个秩 2 数组，等于 A \otimes B。返回数组的形状为 [M1*M2, N1*N2]。

示例

program example_kronecker_product
  use stdlib_linalg, only: kronecker_product
  implicit none
  integer, parameter :: m1 = 1, n1 = 2, m2 = 2, n2 = 3
  integer :: i, j
  real :: A(m1, n1), B(m2,n2)
  real, allocatable :: C(:,:)

  do j = 1, n1
     do i = 1, m1
        A(i,j) = i*j ! A = [1, 2]
     end do
  end do

  do j = 1, n2
     do i = 1, m2      ! B = [ 1, 2, 3 ]
        B(i,j) = i*j !     [ 2, 4, 6 ]
     end do
  end do

  C = kronecker_product(A, B)
  ! C =     [ a(1,1) * B(:,:) | a(1,2) * B(:,:) ]
  ! or in other words, 
  ! C =     [  1.00      2.00      3.00      2.00      4.00      6.00  ]
  !         [  2.00      4.00      6.00      4.00      8.00     12.00  ]
end program example_kronecker_product

`cross_product` - 计算两个向量的叉积

状态

实验性的

描述

计算两个向量的叉积

语法

c = cross_product (a, b)

参数

a: 应为秩 1 且大小为 3 的数组

b: 应为秩 1 且大小为 3 的数组

返回值

返回一个秩 1 且大小为 3 的数组，该数组与 a 和 b 都垂直。

示例

program demo_cross_product
    use stdlib_linalg, only: cross_product
    implicit none
    real :: a(3), b(3), c(3)
    a = [1., 0., 0.]
    b = [0., 1., 0.]
    c = cross_product(a, b)
    !c = [0., 0., 1.]
end program demo_cross_product

`is_square` - 检查矩阵是否为方阵

状态

实验性的

描述

检查矩阵是否为方阵

语法

d = is_square (A)

参数

A: 应为秩为 2 的数组

返回值

如果输入矩阵为方阵，则返回一个logical标量，其值为.true.；否则返回.false.。

示例

program example_is_square
  use stdlib_linalg, only: is_square
  implicit none
  real :: A(2, 2), B(3, 2)
  logical :: res
  A = reshape([1., 2., 3., 4.], shape(A))
  B = reshape([1., 2., 3., 4., 5., 6.], shape(B))
  res = is_square(A) ! returns .true.
  res = is_square(B) ! returns .false.
end program example_is_square

`is_diagonal` - 检查矩阵是否为对角矩阵

状态

实验性的

描述

检查矩阵是否为对角矩阵

语法

d = is_diagonal (A)

参数

A: 应为秩为 2 的数组

返回值

如果输入矩阵为对角矩阵，则返回一个logical标量，其值为.true.；否则返回.false.。请注意，只要a_ij = 0当i /= j时，非方阵也可能为对角矩阵。

示例

program example_is_diagonal
  use stdlib_linalg, only: is_diagonal
  implicit none
  real :: A(2, 2), B(2, 2)
  logical :: res
  A = reshape([1., 0., 0., 4.], shape(A))
  B = reshape([1., 0., 3., 4.], shape(B))
  res = is_diagonal(A) ! returns .true.
  res = is_diagonal(B) ! returns .false.
end program example_is_diagonal

`is_symmetric` - 检查矩阵是否为对称矩阵

状态

实验性的

描述

检查矩阵是否为对称矩阵

语法

d = is_symmetric (A)

参数

A: 应为秩为 2 的数组

返回值

如果输入矩阵为对称矩阵，则返回一个logical标量，其值为.true.；否则返回.false.。

示例

program example_is_symmetric
  use stdlib_linalg, only: is_symmetric
  implicit none
  real :: A(2, 2), B(2, 2)
  logical :: res
  A = reshape([1., 3., 3., 4.], shape(A))
  B = reshape([1., 0., 3., 4.], shape(B))
  res = is_symmetric(A) ! returns .true.
  res = is_symmetric(B) ! returns .false.
end program example_is_symmetric

`is_skew_symmetric` - 检查矩阵是否为反对称矩阵

状态

实验性的

描述

检查矩阵是否为反对称矩阵

语法

d = is_skew_symmetric (A)

参数

A: 应为秩为 2 的数组

返回值

如果输入矩阵为反对称矩阵，则返回一个logical标量，其值为.true.；否则返回.false.。

示例

program example_is_skew_symmetric
  use stdlib_linalg, only: is_skew_symmetric
  implicit none
  real :: A(2, 2), B(2, 2)
  logical :: res
  A = reshape([0., -3., 3., 0.], shape(A))
  B = reshape([0., 3., 3., 0.], shape(B))
  res = is_skew_symmetric(A) ! returns .true.
  res = is_skew_symmetric(B) ! returns .false.
end program example_is_skew_symmetric

`is_hermitian` - 检查矩阵是否为厄米特矩阵

状态

实验性的

描述

检查矩阵是否为厄米特矩阵

语法

d = is_hermitian (A)

参数

A: 应为秩为 2 的数组

返回值

如果输入矩阵为厄米特矩阵，则返回一个logical标量，其值为.true.；否则返回.false.。

示例

program example_is_hermitian
  use stdlib_linalg, only: is_hermitian
  implicit none
  complex :: A(2, 2), B(2, 2)
  logical :: res
  A = reshape([cmplx(1., 0.), cmplx(3., -1.), cmplx(3., 1.), cmplx(4., 0.)], shape(A))
  B = reshape([cmplx(1., 0.), cmplx(3., 1.), cmplx(3., 1.), cmplx(4., 0.)], shape(B))
  res = is_hermitian(A) ! returns .true.
  res = is_hermitian(B) ! returns .false.
end program example_is_hermitian

`is_triangular` - 检查矩阵是否为三角矩阵

状态

实验性的

描述

检查矩阵是否为三角矩阵

语法

d = is_triangular (A,uplo)

参数

A: 应为秩为 2 的数组

uplo: 应为来自{'u','U','l','L'}的单个字符

返回值

如果输入矩阵为uplo（上三角或下三角）指定的三角类型，则返回一个logical标量，其值为.true.；否则返回.false.。请注意，此实现中使用的三角形定义允许非方阵为三角形。具体来说，上三角矩阵满足a_ij = 0当j < i时，而下三角矩阵满足a_ij = 0当j > i时。

示例

program example_is_triangular
  use stdlib_linalg, only: is_triangular
  implicit none
  real :: A(3, 3), B(3, 3)
  logical :: res
  A = reshape([1., 0., 0., 4., 5., 0., 7., 8., 9.], shape(A))
  B = reshape([1., 0., 3., 4., 5., 0., 7., 8., 9.], shape(B))
  res = is_triangular(A, 'u') ! returns .true.
  res = is_triangular(B, 'u') ! returns .false.
end program example_is_triangular

`is_hessenberg` - 检查矩阵是否为 Hessenberg 矩阵

状态

实验性的

描述

检查矩阵是否为 Hessenberg 矩阵

语法

d = is_hessenberg (A,uplo)

参数

A: 应为秩为 2 的数组

uplo: 应为来自{'u','U','l','L'}的单个字符

返回值

如果输入矩阵为uplo（上三角或下三角）指定的 Hessenberg 类型，则返回一个logical标量，其值为.true.；否则返回.false.。请注意，此实现中使用的 Hessenberg 定义允许非方阵为 Hessenberg。具体来说，上 Hessenberg 矩阵满足a_ij = 0当j < i-1时，而下 Hessenberg 矩阵满足a_ij = 0当j > i+1时。

示例

program example_is_hessenberg
  use stdlib_linalg, only: is_hessenberg
  implicit none
  real :: A(3, 3), B(3, 3)
  logical :: res
  A = reshape([1., 2., 0., 4., 5., 6., 7., 8., 9.], shape(A))
  B = reshape([1., 2., 3., 4., 5., 6., 7., 8., 9.], shape(B))
  res = is_hessenberg(A, 'u') ! returns .true.
  res = is_hessenberg(B, 'u') ! returns .false.
end program example_is_hessenberg

`solve` - 求解线性矩阵方程或线性方程组。

状态

实验性的

描述

此函数计算线性矩阵方程的解 $A \cdot x = b$ ，其中 $A$ 是一个方阵、满秩、real或complex矩阵。

结果向量或数组x返回在数值精度范围内对线性方程组的精确解，前提是矩阵不病态。如果矩阵在工作精度下是秩亏或奇异的，则会返回错误。求解器基于 LAPACK 的*GESV后端。

语法

Pure接口

x = solve (a, b)

专家接口

x = solve (a, b [, overwrite_a], err)

参数

a: 应为一个秩为 2 的real或complex方阵，包含系数矩阵。它通常是intent(in)参数。如果overwrite_a=.true.，它是一个intent(inout)参数，并将被调用破坏。

b: 应为与a同种类型的秩为 1 或 2 的数组，包含右侧向量。它是一个intent(in)参数。

overwrite_a（可选）：应为一个输入逻辑标记。如果为.true.，则输入矩阵a将被用作临时存储并被覆盖。这避免了内部数据分配。它是一个intent(in)参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。如果提供此参数，则函数不为pure。

返回值

对于一个满秩矩阵，返回一个数组值，表示线性方程组的解。

如果矩阵在工作精度下是奇异的，则引发LINALG_ERROR。如果矩阵和右侧向量的大小无效或不兼容，则引发LINALG_VALUE_ERROR。如果未提供err，则异常会触发error stop。

示例

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


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

`solve_lu` - 求解线性矩阵方程或线性方程组（子例程接口）。

状态

实验性的

描述

此子例程计算线性矩阵方程的解 $A \cdot x = b$ ，其中 $A$ 是一个方阵、满秩、real或complex矩阵。

结果向量或数组x返回在数值精度范围内对线性方程组的精确解，前提是矩阵不病态。如果矩阵在工作精度下是秩亏或奇异的，则会返回错误。如果用户提供了所有可选数组，则不会发生内部分配。求解器基于 LAPACK 的*GESV后端。

语法

简单（Pure）接口

call solve_lu (a, b, x)

专家（Pure）接口

call solve_lu (a, b, x [, pivot, overwrite_a, err])

参数

a: 应为一个秩为 2 的real或complex方阵，包含系数矩阵。它通常是intent(in)参数。如果overwrite_a=.true.，它是一个intent(inout)参数，并将被调用破坏。

b: 应为与a同种类型的秩为 1 或 2 的数组，包含右侧向量。它是一个intent(in)参数。

x: 应为与b同种类型和大小的秩为 1 或 2 的数组，返回系统解。它是一个intent(inout)参数，并且必须具有contiguous属性。

pivot（可选）：应为与a同种类型且矩阵维数相同的秩为 1 的数组，为对角枢轴索引提供存储空间。它是一个intent(inout)参数，并返回对角枢轴索引。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

对于一个满秩矩阵，返回一个数组值，表示线性方程组的解。

示例

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

`lstsq` - 计算线性矩阵方程的最小二乘解。

状态

实验性的

描述

此函数计算线性矩阵方程的最小二乘解 $A \cdot x = b$ .

结果向量x返回最小化 2 范数的近似解 $|| A \cdot x - b ||_2$ ，即，它包含问题的最小二乘解。矩阵A可以是满秩、超定或欠定。求解器基于 LAPACK 的*GELSD后端。

语法

x = lstsq (a, b, [, cond, overwrite_a, rank, err])

参数

a: 应为一个秩为 2 的real或complex数组，包含系数矩阵。它是一个intent(inout)参数。

b: 应为与a同种类型的秩为 1 或 2 的数组，包含一个或多个右侧向量，每个向量都在其前导维数中。它是一个intent(in)参数。

cond（可选）：应为一个标量real值，用于秩评估的截止阈值：s_i >= cond*maxval(s), i=1:rank。应为标量，intent(in)参数。

overwrite_a（可选）：应为一个输入logical标记。如果为.true.，则输入矩阵A将被用作临时存储并被覆盖。这避免了内部数据分配。它是一个intent(in)参数。

rank（可选）：应为一个integer标量值，包含输入矩阵A的秩。它是一个intent(out)参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

返回与b同种类型和秩的数组值，包含最小二乘系统的解。

如果底层的奇异值分解过程未收敛，则引发LINALG_ERROR。如果矩阵和右侧向量的大小无效或不兼容，则引发LINALG_VALUE_ERROR。异常会触发error stop。

示例

! Least-squares solver: functional interface
program example_lstsq1
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: lstsq
  implicit none

  integer, allocatable :: x(:),y(:)
  real(dp), allocatable :: A(:,:),b(:),coef(:)

  ! Data set
  x = [1, 2, 2]
  y = [5, 13, 25]

  ! Fit three points using a parabola, least squares method
  ! A = [1 x x**2]
  A = reshape([[1,1,1],x,x**2],[3,3])
  b = y

  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
  coef = lstsq(A,b)

  print *, 'parabola: ',coef
  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 


end program example_lstsq1

`solve_lstsq` - 计算线性矩阵方程的最小二乘解（子例程接口）。

状态

实验性的

描述

此子例程计算线性矩阵方程的最小二乘解 $A \cdot x = b$ .

语法

call solve_lstsq (a, b, x, [, real_storage, int_storage, [cmpl_storage, ] cond, singvals, overwrite_a, rank, err])

参数

a: 应为一个秩为 2 的real或complex数组，包含系数矩阵。它是一个intent(inout)参数。

b: 应为与a同种类型的秩为 1 或 2 的数组，包含一个或多个右侧向量，每个向量都在其前导维数中。它是一个intent(in)参数。

x: 应为与b同种类型和秩的数组，前导维数至少为n，包含最小二乘系统的解。它是一个intent(inout)参数。

real_storage（可选）：应为与a同种类型的real秩为 1 的数组，为求解器提供工作存储空间。其最小大小可以通过调用lstsq_space来确定。它是一个intent(inout)参数。

int_storage（可选）：应为一个integer秩为 1 的数组，为求解器提供工作存储空间。其最小大小可以通过调用lstsq_space来确定。它是一个intent(inout)参数。

cmpl_storage（可选）：对于complex系统，它应为一个complex秩为 1 的数组，为求解器提供工作存储空间。其最小大小可以通过调用lstsq_space来确定。它是一个intent(inout)参数。

cond（可选）：应为一个标量real值，用于秩评估的截止阈值：s_i >= cond*maxval(s), i=1:rank。应为标量，intent(in)参数。

singvals（可选）：应为与a同种类型的real秩为 1 的数组，大小至少为min(m,n)，返回内部 SVD 中的奇异值列表s(i)>=cond*maxval(s)，按幅度降序排列。它是一个intent(out)参数。

rank（可选）：应为一个integer标量值，包含输入矩阵A的秩。它是一个intent(out)参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

返回一个数组值，表示最小二乘系统的解。

如果底层的奇异值分解过程未收敛，则引发LINALG_ERROR。如果矩阵和右侧向量的大小无效或不兼容，则引发LINALG_VALUE_ERROR。异常会触发error stop。

示例

! Demonstrate expert subroutine interface with pre-allocated arrays
program example_lstsq2
  use stdlib_linalg_constants, only: dp,ilp
  use stdlib_linalg, only: solve_lstsq, lstsq_space, linalg_state_type
  implicit none

  integer, allocatable :: x(:),y(:)
  real(dp), allocatable :: A(:,:),b(:),coef(:),real_space(:),singvals(:)
  integer(ilp), allocatable :: int_space(:)
  integer(ilp) :: lrwork,liwork,arank

  ! Data set
  x = [1, 2, 2]
  y = [5, 13, 25]

  ! Fit three points using a parabola, least squares method
  ! A = [1 x x**2]
  A = reshape([[1,1,1],x,x**2],[3,3])
  b = y

  ! Get storage sizes for the arrays and pre-allocate data
  call lstsq_space(A,b,lrwork,liwork)
  allocate(coef(size(x)),real_space(lrwork),int_space(liwork),singvals(minval(shape(A))))

  ! Solve coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
  ! with no internal allocations
  call solve_lstsq(A,b,x=coef,              &
                   real_storage=real_space, &
                   int_storage=int_space,   &
                   singvals=singvals,       &
                   overwrite_a=.true.,      &
                   rank=arank)

  print *, 'parabola: ',coef
  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 
  print *, 'rank: ',arank
  ! rank: 2


end program example_lstsq2

`lstsq_space` - 计算最小二乘求解器的内部工作空间要求。

状态

实验性的

描述

此子例程计算最小二乘求解器solve_lstsq的内部工作空间要求。

语法

call lstsq_space (a, b, lrwork, liwork [, lcwork])

参数

a: 应为一个秩为 2 的real或complex数组，包含线性系统系数矩阵。它是一个intent(in)参数。

b: 应为与a同种类型的秩为 1 或 2 的数组，包含系统的右侧向量。它是一个intent(in)参数。

lrwork: 应为一个integer标量，返回此系统所需的real工作存储空间的最小数组大小。

liwork: 应为一个integer标量，返回此系统所需的integer工作存储空间的最小数组大小。

lcwork（complex a，b）：对于complex系统，应为一个integer标量，返回此系统所需的complex工作存储空间的最小数组大小。

`det` - 计算方阵的行列式

状态

实验性的

描述

此函数计算real或complex方阵的行列式。

此接口带有一个pure版本det(a)，以及一个非pure版本det(a,overwrite_a,err)，允许更专业的控制。

语法

c = det (a [, overwrite_a, err])

参数

a: 应为一个秩为 2 的方阵

overwrite_a（可选）：应为一个输入logical标记。如果为.true.，则输入矩阵a将被用作临时存储并被覆盖。这避免了内部数据分配。它是一个intent(in)参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

返回与a同种类型的real标量值，表示矩阵的行列式。

如果矩阵是奇异的，则引发LINALG_ERROR。如果矩阵是非方阵，则引发LINALG_VALUE_ERROR。如果提供了异常，则将返回到err参数；否则，将触发error stop。

示例

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

`.det.` - 方阵的行列式运算符

状态

实验性的

描述

此运算符返回实方阵的行列式。

此接口等效于行列式 det 的 pure 版本。

语法

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

参数

a：应为任何 real 或 complex 类型的秩 2 方阵。它是一个 intent(in) 参数。

返回值

返回表示矩阵行列式的实标量值。

如果矩阵是奇异的，则会引发 LINALG_ERROR。如果矩阵是非方阵，则会引发 LINALG_VALUE_ERROR。异常会触发 error stop。

示例

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

`eig` - 方阵的特征值和特征向量

状态

实验性的

描述

此子例程计算特征问题解 $A \cdot \bar{v} - \lambda \cdot \bar{v}$ ，其中 $A$ 是一个方阵、满秩、real或complex矩阵。

结果数组 lambda 返回 $A$ 的特征值。用户可以请求返回特征向量：如果提供，则在输出时 left 将包含左特征向量，right 将包含 $A$ 的右特征向量。这两个 left 和 right 都是秩 2 数组，其中特征向量存储为列。求解器基于 LAPACK 的 *GEEV 后端。

语法

call eig (a, lambda [, right] [,left] [,overwrite_a] [,err])

参数

a : real 或 complex 方阵，包含系数矩阵。如果 overwrite_a=.false.，则它是一个 intent(in) 参数。否则，它是一个 intent(inout) 参数，并会被调用破坏。

lambda：应为与 a 相同类型的 complex 或 real 秩 1 数组，包含特征值或仅包含其 real 部分。它是一个 intent(out) 参数。

right（可选）：应为与 a 相同大小和类型的 complex 秩 2 数组，包含 a 的右特征向量。它是一个 intent(out) 参数。

left（可选）：应为与 a 相同大小和类型的 complex 秩 2 数组，包含 a 的左特征向量。它是一个 intent(out) 参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

如果计算未收敛，则会引发 LINALG_ERROR。如果任何矩阵或数组具有无效/不兼容的大小，则会引发 LINALG_VALUE_ERROR。如果仅请求 real 部分，但特征值具有非平凡的虚部，则会引发 LINALG_VALUE_ERROR。如果不存在 err，则异常会触发 error stop。

示例

! 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

`eigh` - 实对称或复厄米特方阵的特征值和特征向量

状态

实验性的

描述

此子例程计算特征分解解 $A \cdot \bar{v} - \lambda \cdot \bar{v}$ ，其中 $A$ 是一个方阵、满秩、real 对称 $A = A^T$ 或 complex 厄米特 $A = A^H$ 矩阵。

结果数组 lambda 返回 $A$ 的 real 特征值。用户可以请求返回正交特征向量：在输出时，vectors 可以包含作为列返回的特征向量矩阵。

通常，仅访问 $A$ 的下三角部分。在输入时，logical 标志 upper_a 允许用户请求应使用矩阵的哪个三角部分。

求解器基于 LAPACK 的 *SYEV 和 *HEEV 后端。

语法

call eigh (a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])

参数

a : real 或 complex 方阵，包含系数矩阵。它是一个 intent(in) 参数。如果 overwrite_a=.true.，则它是一个 intent(inout) 参数，并会被调用破坏。

lambda：应为与 a 相同精度的 complex 秩 1 数组，包含特征值。它是一个 intent(out) 参数。

vectors（可选）：应为与 a 相同类型、大小和类型的秩 2 数组，包含 a 的特征向量。它是一个 intent(out) 参数。

upper_a（可选）：应为输入 logical 标志。如果 .true.，将访问 a 的上三角部分。否则，将访问下三角部分。它是一个 intent(in) 参数。

overwrite_a（可选）：应为输入 logical 标志。如果 .true.，输入矩阵 a 将用作临时存储并被覆盖。这将避免内部数据分配。这是一个 intent(in) 参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

如果计算未收敛，则会引发 LINALG_ERROR。如果任何矩阵或数组具有无效/不兼容的大小，则会引发 LINALG_VALUE_ERROR。如果不存在 err，则异常会触发 error stop。

示例

! 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

`eigvals` - 方阵的特征值

状态

实验性的

描述

此函数返回矩阵 $A$ 的特征值：一个方阵、满秩、real 或 complex 矩阵。特征值是特征问题 $A \cdot \bar{v} - \lambda \cdot \bar{v}$ .

的解。结果数组 lambda 是 complex，并返回 $A$ 的特征值。求解器基于 LAPACK 的 *GEEV 后端。

语法

lambda = eigvals (a, [,err])

参数

a : real 或 complex 方阵，包含系数矩阵。它是一个 intent(in) 参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

返回一个包含 a 的特征值的 complex 数组。

示例

! 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

`eigvalsh` - 实对称或复厄米特方阵的特征值

状态

实验性的

描述

此函数返回矩阵 $A$ 其中， $A$ 是一个方阵、满秩、real 对称 $A = A^T$ 或 complex 厄米特 $A = A^H$ 矩阵。特征值是特征问题 $A \cdot \bar{v} - \lambda \cdot \bar{v}$ .

的解。结果数组 lambda 是 real，并返回 $A$ 的特征值。求解器基于 LAPACK 的 *SYEV 和 *HEEV 后端。

语法

lambda = eigvalsh (a, [, upper_a] [,err])

参数

a : real 或 complex 方阵，包含系数矩阵。它是一个 intent(in) 参数。

upper_a（可选）：应为输入逻辑标志。如果 .true.，将使用 a 的上三角部分。否则，将访问下三角部分（默认）。它是一个 intent(in) 参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

返回一个包含 a 的特征值的 real 数组。

示例

! 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

`svd` - 计算秩 2 数组（矩阵）的奇异值分解。

状态

实验性的

描述

此子例程计算 real 或 complex 秩 2 数组（矩阵） $A = U \cdot S \cdot \V^T$ 的奇异值分解。求解器基于 LAPACK 的 *GESDD 后端。

结果向量 s 在 $S$ 的对角线上返回奇异值数组。如果请求，u 包含左奇异向量，作为 $U$ 的列。如果请求，vt 包含右奇异向量，作为 $V^T$ .

语法

的行。call svd (a, s, [, u, vt, overwrite_a, full_matrices, err])

类

子例程

参数

a：应为包含大小为 [m,n] 的系数矩阵的秩 2 real 或 complex 数组。它是一个 intent(inout) 参数，但除非 overwrite_a=.true.，否则它将保持不变。

s：应为秩 1 real 数组，返回 k = min(m,n) 个奇异值的列表。它是一个 intent(out) 参数。

u（可选）：应为与 a 相同类型的秩 2 数组，返回 a 的左奇异向量作为列。除非 full_matrices=.false.，否则其大小应为 [m,m]，在这种情况下，其大小可以为 [m,min(m,n)]。它是一个 intent(out) 参数。

vt（可选）：应为与 a 相同类型的秩 2 数组，返回 a 的右奇异向量作为行。除非 full_matrices=.false.，否则其大小应为 [n,n]，在这种情况下，其大小可以为 [min(m,n),n]。它是一个 intent(out) 参数。

overwrite_a（可选）：应为输入 logical 标志。如果 .true.，输入矩阵 A 将用作临时存储并被覆盖。这将避免内部数据分配。默认情况下，overwrite_a=.false.。它是一个 intent(in) 参数。

full_matrices（可选）：应为输入 logical 标志。如果 .true.（默认），矩阵 u 和 vt 应为全尺寸。否则，它们的次要维度可以调整为 min(m,n)。有关详细信息，请参阅 u、v。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

返回一个包含矩阵 a 的奇异值列表的数组 s。如果请求，返回一个包含 a 的左奇异向量（沿其列）的秩 2 数组 u。如果请求，返回一个包含 a 的右奇异向量（沿其行）的秩 2 数组 vt。

如果底层奇异值分解过程未收敛，则会引发 LINALG_ERROR。如果矩阵或任何输出数组具有无效/不兼容的大小，则会引发 LINALG_VALUE_ERROR。异常会触发 error stop，除非存在参数 err。

示例

`svdvals` - 计算秩 2 数组（矩阵）的奇异值。

状态

实验性的

描述

此子例程根据其奇异值分解计算 real 或 complex 秩 2 数组（矩阵）的奇异值 $A = U \cdot S \cdot \V^T$ 的奇异值分解。求解器基于 LAPACK 的 *GESDD 后端。

结果向量 s 在 $S$ .

语法

s = svdvals (a [, err])

参数

a：应为包含大小为 [m,n] 的系数矩阵的秩 2 real 或 complex 数组。它是一个 intent(in) 参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

返回一个包含矩阵 a 的奇异值列表的数组 s。

示例

! Singular Values
program example_svdvals
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: svdvals
  implicit none

  real(dp), allocatable :: A(:,:),s(:)
  character(*), parameter :: fmt="(a,*(1x,f12.8))"

  ! We want to find the singular values of matrix: 
  ! 
  !    A = [ 3  2  2]
  !        [ 2  3 -2]  
  !   
  A = transpose(reshape([ 3, 2, 2, &
                          2, 3,-2], [3,2]))

  ! Get singular values
  s = svdvals(A)

  ! Singular values: [5, 3]
  print fmt, '    '
  print fmt, 'S = ',s
  print fmt, '    '

end program example_svdvals

`cholesky` - 计算秩 2 方阵（矩阵）的 Cholesky 分解

状态

实验性的

描述

此子例程计算 real 或 complex 秩 2 方阵（矩阵）的 Cholesky 分解， $A = L \cdot L^T$ 或 $A = U^T \cdot U$ . $A$ 是对称或复厄米特矩阵，并且 $L$ , $U$ 分别是下三角或上三角。求解器基于 LAPACK 的 *POTRF 后端。

语法

call cholesky (a, c, lower [, other_zeroed] [, err])

类

子例程

参数

a：应为包含大小为 [n,n] 的系数矩阵的秩 2 方阵 real 或 complex 数组。它是一个 intent(inout) 参数，但如果存在参数 c，则它将保持不变。

c（可选）：应为与 a 相同大小和类型的秩 2 方阵 real 或 complex。它是一个 intent(out) 参数，它返回三角 Cholesky 矩阵 L 或 U。

lower：应为输入 logical 标志。如果 .true.，将执行下三角分解 $A = L \cdot L^T$ 。如果 .false.，将执行上三角分解 $A = U^T \cdot U$ 。

other_zeroed（可选）：应为输入 logical 标志。如果 .true.，输出矩阵的未用部分将包含零。否则，不会访问它。这将节省 CPU 时间。默认情况下，other_zeroed=.true.。它是一个 intent(in) 参数。

err（可选）：应为 type(linalg_state_type) 值。它是一个 intent(out) 参数。

返回值

如果未提供其他参数，则分解后的矩阵将返回到 a 中，覆盖 a。否则，它可以作为第二个参数 c 提供。在这种情况下，不会覆盖 a。logical 标志 lower 决定是否应执行下三角分解或上三角分解。

结果将返回到输出矩阵的适用三角区域，而默认情况下，未用三角区域将填充零。可选参数 other_zeroed（如果为 .false.）允许专家用户避免将未用部分置零；但是，在这种情况下，不会访问矩阵的未用区域，并且通常会包含无效值。

如果底层过程未收敛，则会引发 LINALG_ERROR。如果矩阵或任何输出数组具有无效/不兼容的大小，则会引发 LINALG_VALUE_ERROR。异常会触发 error stop，除非存在参数 err。

示例

! Cholesky factorization: subroutine interface 
program example_cholesky
  use stdlib_linalg, only: cholesky
  implicit none

  real, dimension(3,3) :: A,L,U 

  ! Set real matrix
  A = reshape( [ [6, 15, 55], &
                 [15, 55, 225], &
                 [55, 225, 979] ], [3,3] )

  ! Decompose (lower)
  call cholesky(A, L, lower=.true.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(L,transpose(L))))

  ! Decompose (upper)
  call cholesky(A, U, lower=.false.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(transpose(U),U)))

end program example_cholesky

`chol` - 计算秩 2 方阵（矩阵）的 Cholesky 分解

状态

实验性的

描述

这是一个 pure 函数接口，用于计算 real 或 complex 秩 2 方阵（矩阵）的 Cholesky 分解，计算方法如下： $A = L \cdot L^T$ 或 $A = U^T \cdot U$ . $A$ 是对称或复厄米特矩阵，并且 $L$ , $U$ 分别是下三角或上三角。求解器基于 LAPACK 的 *POTRF 后端。

结果矩阵c与a的大小和类型相同，并返回下三角或上三角因子。

语法

c = chol (a, lower [, other_zeroed])

参数

a: 应为一个大小为[n,n]的秩2方阵，包含系数矩阵，其类型为real或complex。它是intent(inout)参数，但如果参数c存在则返回不变。

lower：应为输入 logical 标志。如果 .true.，将执行下三角分解 $A = L \cdot L^T$ 。如果 .false.，将执行上三角分解 $A = U^T \cdot U$ 。

返回值

返回一个与a大小和类型相同的秩2数组c，其中包含三角形 Cholesky 矩阵L或U。

示例

! Cholesky factorization: function interface 
program example_chol
  use stdlib_linalg, only: chol
  implicit none

  real, allocatable, dimension(:,:) :: A,L,U 

  ! Set real matrix
  A = reshape( [ [6, 15, 55], &
                 [15, 55, 225], &
                 [55, 225, 979] ], [3,3] )

  ! Decompose (lower)
  L = chol(A, lower=.true.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(L,transpose(L))))

  ! Decompose (upper)
  U = chol(A, lower=.false.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(transpose(U),U)))

end program example_chol

`.inv.` - 方阵的逆运算符

状态

实验性的

描述

此运算符返回real或complex方阵的逆。 $A$ 逆 $A^{-1}$ 定义如下： $A \cdot A^{-1} = A^{-1} \cdot A = I_n$ .

此接口等效于函数 inv。

语法

b = [[stdlib_linalg(module):operator(.inv.)(interface)]] a

参数

a：应为任何 real 或 complex 类型的秩 2 方阵。它是一个 intent(in) 参数。

返回值

返回一个与a类型、大小和秩相同的秩2方阵，其中包含a的逆。

如果输入错误或奇异矩阵发生异常，将返回NaN。对于奇异矩阵的细粒度错误控制，请优先使用subroutine和function接口。

示例

! Matrix inversion example: operator interface
program example_inverse_operator
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: operator(.inv.),eye
  implicit none

  real(dp) :: A(2,2), Am1(2,2)

  ! Input matrix (NB Fortran is column major! input columns then transpose)
  A = transpose(reshape( [4, 3, &
                          3, 2], [2,2] ))

  ! Invert matrix
  Am1 = .inv.A

  print *, '         |',Am1(1,:),'|' ! | -2  3 |
  print *, ' inv(A)= |',Am1(2,:),'|' ! |  3 -4 |

  ! Final check 
  print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)

end program example_inverse_operator

`invert` - 方阵求逆

状态

实验性的

描述

此子例程就地求逆一个方阵，其类型为real或complex。逆 $A^{-1}$ 定义如下： $A \cdot A^{-1} = A^{-1} \cdot A = I_n$ .

返回后，输入矩阵a被替换为其逆。求解器基于 LAPACK 的*GETRF和*GETRI后端。

语法

call invert (a, [,inva] [, pivot] [, err])

参数

a: 应为一个秩2方阵，包含系数矩阵，其类型为real或complex。如果提供inva，它是intent(in)参数。如果未提供inva，它是intent(inout)参数：在输出时，它被替换为a的逆。

inva（可选）：应为一个秩2方阵，其大小和类型与a相同。在输出时，它包含a的逆。

pivot（可选）：应为一个秩1数组，其类型和矩阵维度与a相同，在返回时包含对角线主元索引。它是intent(inout)参数。

err（可选）：应为一个type(linalg_state_type)值。它是一个intent(out)参数。

返回值

计算矩阵的逆 $A$ , $A^{-1}\, and returns it either in \( A$ 或在另一个矩阵中。

如果矩阵奇异或大小无效，则引发LINALG_ERROR。如果inva和a大小不同，则引发LINALG_VALUE_ERROR。如果未提供err，异常将触发error stop。

示例

! Matrix inversion example: in-place inversion
program example_inverse_inplace
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: invert,eye
  implicit none

  real(dp) :: A(2,2), Am1(2,2)

  ! Input matrix (NB Fortran is column major! input columns then transpose)
  A = transpose(reshape( [4, 3, &
                          3, 2], [2,2] ))
  Am1 = A

  ! Invert matrix
  call invert(Am1)

  print *, '         |',Am1(1,:),'|' ! | -2  3 |
  print *, ' inv(A)= |',Am1(2,:),'|' ! |  3 -4 |

  ! Final check 
  print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)

end program example_inverse_inplace

! Matrix inversion example: subroutine interface
program example_inverse_subroutine
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: invert,eye
  implicit none

  real(dp) :: A(2,2), Am1(2,2)

  ! Input matrix (NB Fortran is column major! input columns then transpose)
  A = transpose(reshape( [4, 3, &
                          3, 2], [2,2] ))

  ! Invert matrix
  call invert(A,Am1)

  print *, '         |',Am1(1,:),'|' ! | -2  3 |
  print *, ' inv(A)= |',Am1(2,:),'|' ! |  3 -4 |

  ! Final check 
  print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)

end program example_inverse_subroutine

`inv` - 方阵求逆。

状态

实验性的

描述

此函数返回一个方阵的逆，其类型为real或complex，就地求逆。逆 $A^{-1}$ ，定义如下： $A \cdot A^{-1} = A^{-1} \cdot A = I_n$ .

求解器基于 LAPACK 的*GETRF和*GETRI后端。

语法

b inv (a, [, err])

参数

a: 应为一个秩2方阵，包含系数矩阵，其类型为real或complex。它是intent(inout)参数。

err（可选）：应为 type(linalg_state_type) 值。它是一个 intent(out) 参数。

返回值

返回一个与a类型、大小和秩相同的数组值，其中包含逆矩阵 $A^{-1}$ .

如果矩阵奇异或大小无效，则引发LINALG_ERROR。如果未提供err，异常将触发error stop。

示例

! Matrix inversion example: function interface
program example_inverse_function
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: inv,eye
  implicit none

  real(dp) :: A(2,2), Am1(2,2)

  ! Input matrix (NB Fortran is column major! input columns then transpose)
  A = transpose(reshape( [4, 3, &
                          3, 2], [2,2] ))

  ! Invert matrix
  Am1 = inv(A)

  print *, '         |',Am1(1,:),'|' ! | -2  3 |
  print *, ' inv(A)= |',Am1(2,:),'|' ! |  3 -4 |

  ! Final check 
  print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)

end program example_inverse_function

linalg

线性代数

BLAS 和 LAPACK

状态

描述

语法

示例

许可

diag - 创建对角线数组或提取数组的对角线元素

状态

描述

语法

参数

返回值

示例

eye - 构造单位矩阵

状态

类

描述

语法

参数

返回值

警告

示例

trace - 矩阵的迹

状态

描述

语法

参数

返回值

示例

outer_product - 计算两个向量的外积

状态

描述

语法

参数

返回值

示例

kronecker_product - 计算两个秩 2 数组的克罗内克积

状态

描述

语法

参数

返回值

示例

cross_product - 计算两个向量的叉积

状态

描述

语法

参数

返回值

示例

is_square - 检查矩阵是否为方阵

状态

描述

语法

参数

返回值

示例

is_diagonal - 检查矩阵是否为对角矩阵

状态

描述

语法

参数

返回值

示例

is_symmetric - 检查矩阵是否为对称矩阵

状态

描述

语法

参数

返回值

示例

is_skew_symmetric - 检查矩阵是否为反对称矩阵

状态

描述

语法

参数

返回值

示例

`diag` - 创建对角线数组或提取数组的对角线元素

`eye` - 构造单位矩阵

`trace` - 矩阵的迹

`outer_product` - 计算两个向量的外积

`kronecker_product` - 计算两个秩 2 数组的克罗内克积

`cross_product` - 计算两个向量的叉积

`is_square` - 检查矩阵是否为方阵

`is_diagonal` - 检查矩阵是否为对角矩阵

`is_symmetric` - 检查矩阵是否为对称矩阵

`is_skew_symmetric` - 检查矩阵是否为反对称矩阵

`is_hermitian` - 检查矩阵是否为厄米特矩阵

`is_triangular` - 检查矩阵是否为三角矩阵

`is_hessenberg` - 检查矩阵是否为 Hessenberg 矩阵

`solve` - 求解线性矩阵方程或线性方程组。

`solve_lu` - 求解线性矩阵方程或线性方程组（子例程接口）。

`lstsq` - 计算线性矩阵方程的最小二乘解。

`solve_lstsq` - 计算线性矩阵方程的最小二乘解（子例程接口）。

`lstsq_space` - 计算最小二乘求解器的内部工作空间要求。

`det` - 计算方阵的行列式

`.det.` - 方阵的行列式运算符

`eig` - 方阵的特征值和特征向量