我正在尝试用 Fortran 重新编码旧的 C++ 程序以使用 LAPACK（我知道 C++ 确实有 LAPACK++，但我在安装它时遇到了很多麻烦，所以我放弃了）。
我最初没有任何编译问题，但那是当我将所有变量声明为REAL. 当我开始编写需要 LAPACK 的程序部分时，我发现传递给的所有参数都DSYEV需要是DOUBLE PRECISION. 因此，我尝试将所有内容更改为双精度（包括将所有硬编码数字更改为双精度对应数字，即 0.0 -> 0.0D0） 现在，当我尝试编译时，我得到的所有函数和子例程均出现以下错误已经写道：
Error: Return type mismatch of function <function> at (1) (REAL(4)/REAL(8))
我不确定这是从哪里来的，因为程序中的所有内容都已更改为双精度。
例如，我声明了以下内容：
double precision :: alpha(3),d(3),zeta1,zeta2
double precision :: A1(3),A2(3),D1(3),D2(3)
double precision :: PI
PI = 3.14159265359D0
alpha = (/0.109818D0, 0.405771D0, 2.22766D0/)
d = (/0.444635D0, 0.535328D0, 0.154329D0 /)
do 10 i=1,3
A1(i) = alpha(i)*zeta1**2.0D0
A2(i) = alpha(i)*zeta2**2.0D0
D1(i) = d(i)*(2.0D0*A1(i)/PI)**(3.0D0/4.0D0)
D2(i) = d(i)*(2.0D0*A2(i)/PI)**(3.0D0/4.0D0)
10 continue …在 MatLab 中，命令 lu(A) 给出两个矩阵 L 和 U 作为输出，即 A 的 LU 分解。我想知道 Fortran 中是否有某个命令执行完全相同的操作。我一直没能在任何地方找到它。我发现很多 LAPACK 子例程通过首先执行 LU 分解来求解线性系统，但出于我的目的，我需要专门执行 LU 分解并存储 L 和 U 矩阵。
是否有一个命令或子程序以方阵 A 作为输入并以 LU 分解的矩阵 L 和 U 作为输出？
我想使用 dgemm 或任何其他 lapack/blas 函数将矩阵 A 的子集乘以另一个矩阵。我认为由于子矩阵的元素可能不连续，所以我不能直接使用 dgemm 而不将子矩阵复制到另一个空间。所以，当这个子矩阵本身很大时，它可能效率很低，我认为用 C 编写这个特定问题的乘法代码可能会更好。 由于复制然后使用 lapack/blas 本身，可能根本没有效率。我在 matlab 中使用 lapack/blas 作为 mex 文件。
我的问题是
1- lapack/blas 是否有任何功能可以在乘法中处理子矩阵？2-如果不是，直接编写乘法代码更好还是将子矩阵复制到另一个矩阵并使用dgemm更好？
解决Ax = b。真双。 A是超定的 Mx2，其中 M >> 2。 b是 Mx1。我已经针对 运行了大量数据mldivide，结果非常好。我用 MKL 写了一个 mex 例程，LAPACKE_dgels它远没有那么好。结果有大量的噪音，而潜在的信号几乎不存在。我首先根据 MKL 示例结果检查了例程。我已经搜索了mldivide文档（流程图）和 SO 问题。我发现的只是 Matlab 对超定矩形使用 QR 分解。
我接下来应该尝试什么？我是否使用了错误的 LAPACK 例程？请帮助指导我正确的方向。
更新： 在求解向量上的 E-15 浮点差异内，英特尔 MKL LAPACKE_dgels 与 Matlab mldivide 具有相同的结果，用于处理真正的双超定（矩形）问题。据我所知，这是使用的 QR 方法。
当心从这个 dgels 返回的残差。它们不等于 b - Ax。他们中的许多人接近这个值，而有些则远离它。
我试图通过 LAPACK 的 DGETRF 和 DGETRI 例程对矩阵求逆，但以下代码：
program Tester
!use LapackMatrixOps
use MathematicalResources
implicit none
real :: B(2, 2), A(2, 2), WORK(2)
integer :: i, j, SIZ, IPIV(2), INFO
SIZ=2
do i=1, SIZ
do j=1, SIZ
!if(i==j) then
! A(i, j)=1
!else
! A(i, j)=0
!end if
A(i, j)=rand(0)/2+0.5
B(i, j)=0
end do
end do
B=A
call PrintMatrix(A)
call dgetrf(size(B, 1), size(B, 2), B, size(B, 1), IPIV, INFO)
print *, "========="
call PrintMatrix(B)
print *, IPIV
call dgetri(size(B, 1), B, …下面的 C 程序 (dgesv_ex.c)
#include <stdlib.h>
#include <stdio.h>
/* DGESV prototype */
extern void dgesv( int* n, int* nrhs, double* a, int* lda, int* ipiv,
double* b, int* ldb, int* info );
/* Main program */
int main() {
/* Locals */
int n = 10000, info;
/* Local arrays */
/* Initialization */
double *a = malloc(n*n*sizeof(double));
double *b = malloc(n*n*sizeof(double));
int *ipiv = malloc(n*sizeof(int));
for (int i = 0; i < n*n; i++ )
{
a[i] = ((double) …我知道dgemv是矩阵向量,但哪个更有效？使用dgemm直接矩阵乘法,或使用dgemv通过与矩阵B的每一列乘以使用矩阵A做矩阵乘法dgemv？
我有一个问题需要我做许多(~4k)小(~3x3)平方Hermitian矩阵的特征分解和矩阵乘法.特别是,我需要每个工作项来执行一个这样的矩阵的特征分解,然后执行两个矩阵乘法.因此,每个线程必须做的工作相当少,并且完整的工作应该是高度可并行化的.
不幸的是,似乎所有可用的OpenCL LAPACK都用于将大矩阵上的操作委托给GPU,而不是用于在OpenCL内核中进行较小的线性代数操作.由于我不想在OpenCL中为任意大小的矩阵实现矩阵乘法和特征分解,我希望这里有人可能知道这个工作适合的库？
我知道OpenCL可能会在某些时候获得内置矩阵运算,因为矩阵类型是保留的,但现在这并没有多大用处.这里有一个类似的问题,从2011年开始,但它只是说要自己动手,所以我希望从那时起情况有所改善.
我有一个有效的LAPACK实现,据我所知,它包含BLAS.
我想使用稀疏BLAS而据我了解这个网站,稀疏BLAS是BLAS的一部分.
但是当我尝试使用稀疏blas手册运行下面的代码时
g ++ -o sparse.x sparse_blas_example.c -L/usr/local/lib -lblas && ./sparse_ex.x
编译器(或链接器？)要求blas_sparse.h.当我把那个文件放在工作目录中时,我得到了:
ludi@ludi-M17xR4:~/Desktop/tests$ g++ -o sparse.x sparse_blas_example.c -L/usr/local/lib -lblas && ./sparse_ex.x
In file included from sparse_blas_example.c:3:0:
blas_sparse.h:4:23: fatal error: blas_enum.h: No such file or directory
#include "blas_enum.h"
使用SPARSE BLAS和LAPACK我该怎么办？我可以开始将很多头文件移动到工作目录中,但我收集了我已经将它们与lapack一起使用了!
/* C example: sparse matrix/vector multiplication */
#include "blas_sparse.h"
int main()
{
const int n = 4;
const int nz = 6;
double val[] = { 1.1, 2.2, 2.4, 3.3, 4.1, 4.4 };
int indx[] = { …对于方阵，可以得到SVD
X= USV'
通过简单地使用numpy.linalg.svd进行分解
u,s,vh = numpy.linalg.svd(X)
例程或numpy.linalg.eigh，以计算埃尔米特矩阵X'X 和XX'的eig分解
他们使用相同的算法吗？调用相同的Lapack例程？
在速度方面有什么不同吗？和稳定性？