Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - Printable Version +- Python Forum (https://python-forum.io) +-- Forum: Python Coding (https://python-forum.io/forum-7.html) +--- Forum: General Coding Help (https://python-forum.io/forum-8.html) +--- Thread: Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) (/thread-25618.html) |
Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - DimosG - Apr-05-2020 Hi guys I have a problem with defining A (its in line 33) in my Full matrix solution and i was wondering if you could help me . class Vector(object): def __init__(self, data): self.data = data self.rows = len(data) def __mul__(self, other): assert self.rows == other.rows return sum([vi * vj for vi, vj in zip(self.data, other.data)]) def __rmul__(self, a): return Vector([a * d for d in self.data]) def __add__(self, other): assert self.rows == other.rows return Vector([i + j for (i, j) in zip(self.data, other.data)]) def __sub__(self, other): assert self.rows == other.rows return Vector([i - j for (i, j) in zip(self.data, other.data)]) def norm(self): return math.sqrt(self * self) def __str__(self): return '{0}'.format(self.data) class Matrix(object): def __init__(self, rows, cols): self.rows = rows self.cols = cols def cg(self, b, epsilon=1e-4): ... # [1] xk = Vector([0.] * b.rows) rk = b - A * xk pk = 1. * rk for k in range(b.rows): alpha = (rk * rk) / (pk * (A * pk)) xkk = xk + alpha * pk rkk = rk - alpha * (A * pk) if rkk.norm() < epsilon: break beta = (rkk * rkk) / (rk * rk) pkk = rkk + beta * pk rk = rkk pk = pkk xk = xkk return xkk #exit() class FullMatrix(Matrix): def __init__(self, rows, cols, data): super().__init__(rows, cols) self.data = data def __mul__(self, v): assert self.cols == v.rows data = [0.] * self.rows ... # [2] for i in range(self.rows): for j in range(self.cols): data[i] += self.data[i+j*self.rows]*v.data[j] return Vector(data) #exit () class SparseMatrix(Matrix): def __init__(self, rows, cols, triplets): super().__init__(rows, cols) self.triplets = triplets def __mul__(self, v): assert self.cols == v.rows data = [0.] * self.rows ... # [3] for (i,j,d) in self.triplets: data[i] += d*v.data[j] return Vector(data) #exit() if __name__ == "__main__": sizes = [2 ** i for i in (2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)] t_full = [] t_sparse = [] for N in sizes: print('_' * 80) print("size : {0}".format(N)) # b - vector data = [0. for i in range(N)] data[-1] = 10 b = Vector(data) # Af - full matrix data = [0.] * (N * N) for i in range(N - 1): data[i + N * i ] = 10. data[i + 1 + N * i ] = -0.1 data[i + N * (i + 1)] = -0.1 data[N - 1 + N * (N - 1)] = 10. Af = FullMatrix(N, N, data) # solve t0 = time.time() #~ xf = cg(Af, b) xf = Af.cg(b) tf = time.time() - t0 t_full.append(tf) print("full : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(tf, xf.data[-1])) # As - sparse triplets = [] for i in range(N - 1): triplets.append((i, i, 10.)) triplets.append((i + 1, i, -0.1)) triplets.append((i, i + 1, -0.1)) triplets.append((N - 1, N - 1, 10.)) As = SparseMatrix(N, N, triplets) # solve t0 = time.time() #~ xs = cg(As, b) xs = As.cg(b) ts = time.time() - t0 t_sparse.append(ts) print("sparse : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(ts, xs.data[-1])) #exit() # An - numpy An = np.zeros((N, N)) for i in range(N - 1): An[i, i] = 10 An[i + 1, i] = -0.1 An[i, i + 1] = -0.1 # solve An[N - 1, N - 1] = 10. t0 = time.time() xn = np.linalg.solve(An, b.data) tn = time.time() - t0 print("numpy : time: {0:6.2f}s (x[-1] = {1:12.10g})".format(tn, xn[-1])) print() fig = plt.figure() ax = fig.add_subplot(111) ax.plot(sizes, t_full, 'b-o', label='full') ax.plot(sizes, t_sparse, 'r-o', label='sparse') ax.set_xlabel('size') ax.set_ylabel('times [s]') ax.grid() ax.legend() plt.show() RE: Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - micseydel - Apr-07-2020 Ideally you would be more specific about what you're looking for, and simplify your code so that it's only relevant to your question. RE: Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} ) - 1968Edwards - Sep-18-2021 # Who solve sparse matrix with A sparce: A.x = f # 2x + 3y + 0 = 6 # 0 + 0 + 2z = 4 # 4x + 0 + 1z = 2 #This is an example of sparse matrix above # So I'm going to name rows, columns and data import numpy as np from scipy.sparse import bsr_matrix row_bsr = np.array([0, 0, 1, 2, 2]) col_bsr = np.array([0, 1, 2, 0, 2]) data_bsr = np.array([2, 3, 2, 4, 1]) A = bsr_matrix((data_bsr, (row_bsr, col_bsr)), shape = (3,3) ) print(A) # how to declare "b" in this case # if I print A result OK print(A) # but I need x = A⁻1 f # how do i proceed to extract the x #=> This is my question because the scipy documentation doesn't show how to solve this. # have many way: # 1. bsr_matrix: Block Sparse Row matrix # 2. coo_matrix: COOrdinate format matrix # 4. csc_matrix: Compressed Sparce Column matrix # 5. dia_matrix: Sparce matrix with DAgoal storage # 6. dok_matrix: Dictionary Of Key based sparce matrix # 7. lil_matrixÇ Row-based linked list sparse matrix. # How to solve this matrix to extract the x. # Common array resolution I can solve in Python but sparse # arrays can't find the proper codes. The example above is #pretty simple but my problem is 29x29 with a lot of zeros #in A and only two values in f. |