我用两种不同的方式进行QR分解:使用标准numpy方法和使用CULA库中实现的GEQRF-LAPACK函数。下面是python中的一个简单示例(PyCULA用于访问CULA):
from PyCULA.cula import culaInitialize,culaShutdown
from PyCULA.cula import gpu_geqrf, gpu_orgqr
import numpy as np
import sys
def test_numpy(A):
Q, R = np.linalg.qr(A)
print "Q"
print Q
print "R"
print R
print "transpose(Q)*Q"
print np.dot(np.transpose(Q), Q)
print "Q*R"
print np.dot(Q,R)
def test_cula(A):
culaInitialize()
QR, TAU = gpu_geqrf(A)
R = np.triu(QR)
Q = gpu_orgqr(QR, A.shape[0], TAU)
culaShutdown()
print "Q"
print Q
print "R"
print R
print "transpose(Q)*Q"
print np.dot(np.transpose(Q), Q)
print "Q*R"
print np.dot(Q,R)
def main():
rows = int(sys.argv[1])
cols = int(sys.argv[2])
A = np.array(np.ones((rows,cols)).astype(np.float64))
print "A"
print A
print "NUMPY"
test_numpy(A.copy())
print "CULA"
test_cula(A.copy())
if __name__ == '__main__':
main()
它产生以下输出:
^{pr2}$我的代码怎么了?在
这是一个棘手的问题,这里的问题是Python使用Row主顺序,但是CULA使用列major order作为R。请查看CULA文档以了解更多详细信息。在
以下是scikit cuda的示例:
如果改为使用(这是Python中的默认设置)
^{pr2}$你会得到和你在上面发布的相同的错误结果。在
这个问题可能导致几个问题,所以你必须非常小心,并注意矩阵的顺序。在
干杯, 本
我在R.CULA中测试了您的示例。CULA似乎提供了与R相同的结果。下面是我的代码:
你的矩阵:
^{pr2}$以下是来自CULA的结果:
以下是R的结果:
我希望这有帮助!在
相关问题 更多 >
编程相关推荐