Conjugate Gradient Solver (CG) using Python
Posted: Mon Mar 11, 2019 5:39 am
I am using Python and currently solve my problem with the "scipy.sparse.linalg.cg" lib and "scipy.sparse.linalg.LinearOperator". To speed up my calculation (>>1 mio. DOF), I want to use MAGMA. I am quite new to Python.
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?
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