我试图使用python中的ctypes在BLAS中使用函数sgemm
。试图解决C = A x B下面的代码工作得很好:
no_trans = c_char("n")
m = c_int(number_of_rows_of_A)
n = c_int(number_of_columns_of_B)
k = c_int(number_of_columns_of_A)
one = c_float(1.0)
zero = c_float(0.0)
blaslib.sgemm_(byref(no_trans), byref(no_trans), byref(m), byref(n), byref(k),
byref(one), A, byref(m), B, byref(k), byref(zero), C, byref(m))
现在我想解决这个方程:C = A' x A其中A'是A的转置,下面的代码运行没有异常,但返回的结果是错误的:
trans = c_char("t")
no_trans = c_char("n")
m = c_int(number_of_rows_of_A)
n = c_int(number_of_columns_of_A)
one = c_float(1.0)
zero = c_float(0.0)
blaslib.sgemm_(byref(trans), byref(no_trans), byref(n), byref(n), byref(m),
byref(one), A, byref(m), A, byref(m), byref(zero), C, byref(n))
对于一个检验,我插入一个矩阵 a = [1 2;3 4] 。正确的结果是C = [10 14;[14 20]但sgemm
例程吐出C = [5 11];11 25) . 据我所知,矩阵A不必由我转置,因为算法会处理它。在第二种情况下,我的参数传递有什么问题?
任何帮助,链接,文章,建议是感激的!
Blas通常使用列主矩阵(如Fortran),因此A = [1 2; 3 4]
表示
|1 3|
A = | |
|2 4|
,结果是正确的(假设您的Python库也这样做)。read-me
您得到的结果表明sgemm
计算的是A*A'而不是您想要的A'*A。简单的解决方案是将两个输入切换到函数。