def matvec(A,v):
M,N = A.shape
if N != M:
print(“Matrix is not square”)
return -1
N1 = v.shape[0]
if N1 != M:
print("The dimension of matrix and vector is not reasonable ")
result = np.zeros((N1,1),dtype = np.float32)
for i in range(N):
for j in range(N):
result[i] += A[i,j]*v[j]
return result
def product(v1,v2):
result = 0
M,N = v1.shape[0],v2.shape[0]
if N != M:
print(“The dimension of vectors is not reasonable”)
exit(1)
for i in range(M):
result += v1[i]*v2[i]
return result
def vecnorm2(v1):
N = v1.shape[0]
result = 0
for i in range(N):
result += v1[i]*v1[i]
result = np.sqrt(result)
return result
def cgsolver(A,b,x0,maxit:int,tol=1.0e-8):
if tol < 0:
print(“argument tol is not good (negitive)”)
if maxit <= 0:
maxit