Cython lapack 不会覆盖参数

Question

我正在学习 Cython lapack 界面。我尝试将一个在 C 中使用 lapack 的小示例转换为在 Cython 中使用 cython_lapack 包的一个小示例。

这里是 C 代码：

#include <stdio.h>
#include <lapacke.h>

int main()
{
  int i, j;
  float norm;
  float A[] = {1.44, -7.84, -4.39, 4.53, 
               -9.96, -0.28, -3.24,  3.83,
               -7.55, 3.24, 6.27, -6.64,
               8.34, 8.09, 5.28, 2.06,
               7.08, 2.52, 0.74, -2.47,
               -5.45, -5.70, -1.19, 4.70};
  float B[] = {8.58, 9.35,
               8.26, -4.43,
               8.48, -0.70,
               -5.28, -0.26,
               5.72, -7.36,
               8.93, -2.52};
  lapack_int m = 6, n = 4, nrhs = 2;
  int info;

  info = LAPACKE_sgels(LAPACK_ROW_MAJOR, 'N', m, n, nrhs, A, n, B, nrhs);
  printf("Least squares solution:\n");
  for (i = 0; i < n; ++i) {
    for (j = 0; j < nrhs; ++j) printf(" %6.2f", B[nrhs * i + j]);
    printf("\n");
  }
  printf("Residual sum of squares for the solution:\n");
  for (i = 0; i < nrhs; ++i) {
    norm = 0.0;
    for (j = n; j < m; ++j) norm += B[nrhs * j + i] * B[nrhs * j + i];
  printf(" %6.2f", norm);
 }
 printf("\n");
 return 0;
}

我翻译成Cython，如下图：

import numpy as np
cimport numpy as cnp
from scipy.linalg.cython_lapack cimport dgels

a = np.array([1.44, -9.96, -7.55, 8.34, 7.08, -5.45,
              7.84, -0.28, 3.24, 8.09, 2.52, -5.70,
              -4.39, -3.24, 6.27, 5.28, 0.74, -1.19,
              4.53, 3.83, -6.64, 2.06, -2.47, 4.70], dtype=np.float64)
b = np.array([8.58, 8.26, 8.48, -5.28, 5.72, 8.93,
              9.35, -4.43, -.70, -0.26, -7.36, -2.52], dtype=np.float64)

cdef int m = 6, n = 4, nrhs = 2, lda = 6, ldb = 6, lwork = -1, info = 1
cdef double work = 0.0

cdef double* A = <double*>(cnp.PyArray_DATA(a))
cdef double* B = <double*>(cnp.PyArray_DATA(b))

dgels('N', &m, &n, &nrhs, A, &lda, B, &ldb, &work, &lwork, &info)
print "A:"
for j in range(4):
    for i in range(6):
        print "%6.2f " % A[6*j+i],
    print ""
print "B:"
for j in range(2):
    for i in range(6):
        print "%6.2f " % B[6*j+i],
    print ""

我期待 A 和 B 被修改，就像 C 代码中的行为一样。然而，奇怪的是它们根本没有改变。我想知道我应该怎样才能得到正确的结果。

Answer 1

来自the Netlib documentation

如果 LWORK = -1，则假定为工作区查询；例行公事 只计算 WORK 数组的最佳大小，返回 这个值作为 WORK 数组的第一个条目，并且没有错误 与 LWORK 相关的消息由 XERBLA 发出。

您没有要求进行计算 - 您只是询问了应该使工作区数组有多大。你可以

1. 使用适当大小的工作区数组进行第二次调用，或者
2. 只需根据文档分配一个工作区数组（最好是max( 1, MN + max( MN, NRHS )*NB )）并将其大小作为LWORK 传递。