Python Forum
Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} )
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Conjugate Gradient having issues with defining A (function to solve [A]{x} = {b} )
#1
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()
Reply
#2
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.
Reply
#3
# 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.
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  how can I solve fsolve function error? troddydeeneeeeeey 3 2,386 Oct-14-2021, 07:07 PM
Last Post: deanhystad
  Defining a function with input abcd 5 3,041 Feb-21-2021, 02:34 AM
Last Post: NullAdmin
  subprogram issues: cannot unpack non-iterable function object error djwilson0495 13 5,881 Aug-20-2020, 05:53 PM
Last Post: deanhystad
  Issues with storing variables outside of a function cerulean747 7 3,640 Apr-30-2020, 08:46 AM
Last Post: DeaD_EyE
  When Defining a Function with an Equation as a Default Argument, which Value Is Used? OJGeorge4 4 2,609 Apr-09-2020, 08:48 AM
Last Post: DeaD_EyE
  Gradient background swisha 0 2,340 Mar-04-2020, 03:07 AM
Last Post: swisha
  How to solve a function and get x, y and z? TheZenMan 2 2,242 Mar-22-2019, 01:10 PM
Last Post: scidam
  Problem with defining a function snow_y 4 3,156 Nov-26-2018, 02:13 AM
Last Post: snippsat
  main function scope issues wak_stephanie 1 2,431 Aug-29-2018, 02:53 AM
Last Post: Larz60+
  IDLE indenting 29 columns after defining function Cosmo_Kane 1 2,399 Jun-03-2018, 08:53 AM
Last Post: Larz60+

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020