关于 Fortran 中 SGESV 语法的问题

Question

我对这个子程序有点困惑。我已经阅读了文档，但我有点困惑 IPIV 向量到底是做什么的，以及我是如何设置我的前导维度的。我读到前导维度有助于在数组的每个连续列中找到矩阵元素的起点。例如，假设我们要解决

Ax = B

在哪里

integer, parameter :: sp = selected_real_kind(6,37)
real(kind=sp),dimension(:,:),intent(inout) :: A
real(kind=sp),dimension(:),intent(inout) :: B
integer, dimension(10) :: IPIV

其中 sp 是我在主程序中设置的单精度 尺寸是

A(10,10)
B(10)

在我的主程序中设置并传递给这个子程序 我应该将我的子程序设置为

integer :: n,INFO
n = size(A,1)


IPIV = 0
call SGESV(n,n,A,2*n,IPIV,B,2*n,INFO)

或

call SGESV(n,n,A,n,IPIV,B,n,INFO)

对于 IPIV，我应该只创建一个大小为 10 的向量并用零​​初始化它？

---

编辑：我用过

call sgesv(n, n, A, n, ipiv, B, n, INFO)

按照建议，但我收到分段错误程序收到信号 SIGSEGV: Segmentation fault - invalid memory reference。

我已经打印了矩阵大小，它们是正确的，矩阵 A 的大小是 100，向量的大小是 10

Edit2：所以在我的 main 中，我有一个循环，在我的循环中，它在每次迭代时计算 A (10,10) 的矩阵和向量 B(10)。然后我调用我的子程序来解决系统

call solver(A,B)

但是，由于尺寸正确，我得到了我不理解的分段错误。 （为了检查它，我打印了矩阵和向量的大小并注释掉了对我的子程序的调用，它们是 100 和 10）

也许我应该让我的矩阵可分配？但我看不出有什么问题，因为在每次迭代中，我都会通过一系列计算来计算矩阵和向量并覆盖它们。

基本上我如下声明矩阵和向量

   real(sp) , dimension (10) :: B
   real(sp) , dimension (10,10) :: A

然后在我的循环中执行一系列计算以用值填充它们 然后我调用我的子程序 然后用新值重复

Answer 1

您正在使用旧界面进行 lapack。请注意我对现代/通用例程的较低答案。

---

旧界面 你会这样称呼它

call sgesv(n, n, A, n, ipiv, B, n, info)

推理：

• 主要尺寸是n 而不是2n
• ipiv 是一个输出变量 s.t.你不需要用 0 初始化它

---

现代界面：LAPACK95 只使用提供通用调用的现代接口会容易得多

call gesv(A, B, ipiv=ipiv, info=info)

您不需要指定数据类型（例如，没有更多的sgesv）或矩阵维度。

确保您需要使用适当的模块

use lapack95

Answer 2

下面是调用gesv 的示例，它是sgesv 和dgesv 的通用Lapack95 等价物（而且更简单）。

subroutine test_lapack95(n)
use BLAS95
use LAPACK95
use f95_precision
implicit none

integer, intent(in) :: n
real(float), allocatable :: A(:,:), LU(:,:)
real(float), allocatable :: b(:), x(:)
integer, allocatable :: ipiv(:)

allocate(A(n,n))
allocate(b(n))
allocate(ipiv(n))

! Fill values in A and b
call prepare_values(n, A, b)
        
LU = A
x = b
call gesv(LU,x,ipiv)
! Solve to A*x=b, for x

end subroutine

不用担心辅助函数prepare_values，它只需填写A和b。