我正在尝试使用最小二乘法来解决矩阵问题：b = A*x for x。系统是超定的，A 是一个稠密矩阵。
在 LAPACK 库中，我相信例程 DGELSY 最适合这个问题（或最接近 Matlab 的 LSQMINNORM 函数的任何东西）。但是，我是编码方面的相对业余爱好者，尤其是 Fortran，并且对此函数的输入存在问题。
你的语法如下：
call dgelsy(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info)
（这里是文档站点：https ://software.intel.com/en-us/mkl-developer-reference-fortran-gelsy ）。
例如，令 A 为 40,000 x 3,000 矩阵，b 为 40,000 元素向量，x 为 3,000 元素向量。
问题：在这种情况下，DGELSY 的输入应该是什么？
[旁白：如果有更适合解决此问题的 LAPACK 功能，请随时告诉我！]
你大部分都在那儿。
a是矩阵 A 并且b是右侧向量或矩阵 b。当子例程完成时，a已经完全改变（它现在包含 QR 分解）并且b现在包含结果 x。jpvt如您所想，它不包含输出矩阵 x，它包含您可能不需要的旋转信息。rcond了一些小的数字，比如1.e-8work所需的工作区数组。DGELSYLAPACK 应该非常快，而内存分配非常慢，但是这个子程序需要额外的空间来进行计算。它需要您分配空间，而不是分配空间本身。如果您要多次调用此子例程，那么您只分配一次内存是有意义的，而不是DGELSY每次都在内部进行。lwork（以及相应的大小work）可能会令人困惑。您可以DGELSY先在“工作区查询”模式下调用，以获取 的建议值lwork，或将其设置为足够大的值。对于您的情况lwork=30000应该很多。查看文档以获取更多信息。lda和ldb通常与a和相同b。如果您在更大矩阵的子块上操作，它们就在那里。由于数据必须采用列优先格式，这告诉 LAPACK 在读/写数据列之间必须“跳过”多少值。