In this .pdf (p. 172 and following) ways to use several linear solvers from MAGMA in Python in a straight forward manner are described. Unfortunately "Conjugate Gradient Solver" is not listed (but exists in MAGMA):
https://www.researchgate.net/profile/An ... Python.pdf
My question, is it possible to use the CG MAGMA implementation in Python the way it is described in the pdf? Or are more complex workarounds necessary?
Code: Select all
from time import time # import time
import numpy as np # import numpy
import skcuda.magma as magma # import scikit magma
4.2 LU decomposition and solving general linear systems
176
magma.magma_init () # init magma
m=4000 # matrix dimension
a=np.random.rand(m,m). astype(np.complex64) #random coeff.matr
nrhs=1 # number of rhs
b=np.ones((m,nrhs),np.complex64) # vector of ones
B=np.dot(a.T,b) # rhs: B=a*b
ipiv=np.ones(m,np.int32) # pivots
t0=time() # start timer
# solve the linear system A*x=B, B-mx1 matrix , a-mxm matrix;
# B is overwritten by the solution; syntax:
# magma_cgesv(n, nhrs , A, lda , ipiv , B, ldb)
status = magma.magma_cgesv(m, nrhs, a.ctypes.data, m, --------------------------------> Use CG here???
ipiv.ctypes.data, B.ctypes.data, m)
t=time()-t0 # stop timer
print "GPU time:", t, "sec." # print elapsed time
#print B[:5] # to see the real errors
print "max error:" # max difference expected/computed sol.
print np.max(np.abs(B-b))
magma.magma_finalize () # finalize magma