I have been assigned to use the scipy sparse matrix solver, scipy.sparse.linalg.spsolve, to solve a sparse matrix-vector problem.
The problem is I do not understand how to set up the question properly. Once the matrix can be set up the data analysis and GPU acceleration of the code should be easy. I'm getting stuck at the first hurdle because of the notation.
Here is the Matrix notation problem I'm trying to solve
Here is what I have attempted so far:
def matrix_A(N):
nelements = 1 / N
row_ind = []
col_ind = []
data = []
f = np.zeros(N, dtype=np.float64)
count = 0
for i in range(N):
for j in range(N):
if i==j:
row_ind[count] = col_ind[count] = 1
if i==0 or i==N:
if i==N:
row_ind[count] = col_ind[count] = 1
if i==j:
row_ind[count] = col_ind[count] = 2-(h**2 * k**2)
if j==i+1 or j==i-1:
row_ind[count] = col_ind[count] = -1
else:
row_ind[count] = col_ind[count] = 0
return coo_matrix((data, (row_ind, col_ind)), shape=(N**2, N**2)).tocsr(), f
Well, the first part says that A[i,j] = 1 where i = j, which means your diagonal elements are (from what we know so far) all 1. In other words, you can start by creating the identity matrix.
import numpy as np
N = ...
h = ...
k = ...
A = np.eye(N, N)
The next part pertains to i = 0, N-1 (the image says N, but Python is 0-indexed and they probably mean the last value, which is index N-1). With that, it says A[i,j] = 2 - h^2*k^2 if i = j. Since we are referring to specific i values, this can be interpreted as when i = j = 0, N-1.
A[0, 0] = 2 - (h*k)**2
A[-1, -1] = 2 - (h*k)**2
It also says that A[i,j] = -1 if j = i + 1 or j = i - 1. Again, because we are dealing with i = 0, N-1, we only consider A[0,1] and A[N-1,N-2]. We ignore j = i - 1 when i = 0 and j = i + 1 when i = N - 1 because they would give values for j that are out of range.
A[0, 1] = A[-1, -2] = -1
After that, you can convert your matrix to a sparse matrix.