Загрузка процессора при использовании Python Scipy разреженного решателя для системы Ax = b - PullRequest
Новая
       61

Загрузка процессора при использовании Python Scipy разреженного решателя для системы Ax = b

0 голосов
/ 26 сентября 2019

У меня есть решатель для PDE, который требует решения линейной системы вида Ax = b, где A - трехдиагональная матрица, а b - вектор.

Матрица A генерируется для 2 сценариев, в которых PDEпредставляет интерес в двухмерном пространстве и трехмерном пространстве.Код для построения двухмерного пространства: 

def linear_crank_nicolson(f0, wavelength, n0, dx = None, dz = None, layer = None, layer_width = 100):
        """
        Equation: d_t f = i d_xx f

        Method implementation M_left f(n+1) = M_right f(n)
        """
        if dx == None:
            dx = p.dx
        if dz == None:
            dz = p.dz
        if wavelength == None:
            wavelength = p.wavelength
        if np.size(f0) == p.N_x:
            Nx = p.N_x
        else:
            Nx = np.size(f0)

        matrix_minus = sp.lil_matrix((Nx, Nx), dtype = 'complex')
        matrix_plus = sp.lil_matrix((Nx, Nx), dtype = 'complex')


        if n0 == 0:
            print(r'n_0=0')
        else:
            r = 1j*dz/2./dx**2 * wavelength/n0/4./np.pi  

            matrix_minus[0,0] = 1.        
            matrix_minus[Nx-1,Nx-1] = 1.

            matrix_plus[0,0] = 1.        
            matrix_plus[Nx-1,Nx-1] = 1.

            for i in range(Nx-2):
                matrix_minus[i+1,i+1] = 1+2*r        
                matrix_plus[i+1,i+1] = 1-2*r

                matrix_minus[i+1,i]=-r
                matrix_minus[i+1,i+2]=-r

                matrix_plus[i+1,i]=r
                matrix_plus[i+1,i+2]=r

        matrix_minus = sp.csc_matrix(matrix_minus)
        matrix_plus = sp.csc_matrix(matrix_plus)

        return matrix_minus, matrix_plus

Код для трехмерного сценария: 

def linear_crank_nicolson(wavelength, n0, dx = None, dz = None):
        Nx = p.N_x
        Ny = p.N_y
        if dx == None:
            dx = p.dx
        if dz == None:
            dz = p.dz
        matrix_minus = sp.lil_matrix((Nx * Ny, Nx * Ny), dtype = 'complex')
        matrix_plus = sp.lil_matrix((Nx * Ny, Nx * Ny), dtype = 'complex')

        if n0 == 0:
            print(r'n_0=0')
        else:
            r = 1j*dz/2./dx**2 * wavelength/n0/4./np.pi  

            for i in range(Nx * Ny):
                matrix_minus[i,i] = 1+4*r        
                matrix_plus[i,i] = 1-4*r

            matrix_minus[0,0] = 1.        
            matrix_minus[Nx * Ny-1,Nx * Ny-1] = 1.

            matrix_plus[0,0] = 0.        
            matrix_plus[Nx * Ny-1,Nx * Ny-1] = 1.

            for i in range(Nx * Ny-2):
                #print('Matrices init step: ' + str(i) + '/' + str(Nx * Ny))
                matrix_minus[i+1,i]=-r
                matrix_minus[i+1,i+2]=-r

                if i>=Ny and i<(Nx-1) * Ny-1:
                    matrix_minus[i+1, i+1-Ny] = -r
                    matrix_minus[i+1, i+1+Ny] = -r

                matrix_plus[i+1,i]=r
                matrix_plus[i+1,i+2]=r

                if i>=Ny and i<(Nx-1) * Ny-1:
                    matrix_plus[i+1, i+1-Ny] = r
                    matrix_plus[i+1, i+1+Ny] = r

        matrix_minus = sp.csc_matrix(matrix_minus)
        matrix_plus = sp.csc_matrix(matrix_plus)

        return matrix_minus, matrix_plus

Теперь решатели для каждого сценария идентичны и выглядят так: 

while i<steps:
    solution[i,:] = f0
    f0 = np.transpose(np.matrix(f0))
    f0 = matrix_plus * f0
    f0 = spl.spsolve(matrix_minus, f0)
    f0 = np.array(np.transpose(f0))

Проблема, с которой я столкнулся, заключается в том, что всякий раз, когда я запускаю 2-мерный решатель, используется только 1 поток, тогда как для 3-мерного решателя используются все потоки.Кроме того, решатель работает значительно медленнее, когда я запускаю трехмерный решатель.

Я проверил на аналогичные размеры матрицы A, которая составляет около 10000 x 10000, и каждый раз появляются одинаковые результаты.

Я не против того, что решатель работает медленнее в трехмерном сценарии, поскольку матрица менее плотная, однако я не уверен, почему в двумерном случае для подобных или даже на 4 порядка больших размеров решатель быстрееИ это именно потому, что он реализован одинаково для обоих сценариев.

На данный момент я лучше всего предполагаю, что разреженность матрицы влияет на производительность, но это не объясняет разницу в количестве потоков, которые

РЕДАКТИРОВАТЬ: в отношении импорта:

  • spl = scipy.sparse.linalg
  • sp = scipy.sparse
  • p =какой-то файл, который я использую для импорта параметров моих симуляций
...