Умножение разреженная матрица

@urung · Регистрация: 19.07.2014

Студворк — интернет-сервис помощи студентам

У меня проблема в модуле в котором производится функцией dot умножение матрицы эта функция не поддерживает разреженные матрицы(так мне кажется). Я глубокий чайник и при этом пенсионер так что случай тяжелый. Вот модуль которого надо заставить заработать

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
__author__ = 'kardashevskiy'
import numpy as np
from scipy import sparse
 
def sg1(B,b,tol=1.0e-8):
    """
    Solve the linear system Ax=b by conjugate gradient method
    """
    n=len(b)
    x=np.zeros((n),'float')
    r=np.copy(b)
    for i in range(n):
        r[i]=np.dot(B[i,0:n],x[0:n])-b[i]
    s=np.copy(r)
    As=np.zeros((n),'float')
    for k in range(n):
        for i in range(n):
            As[i]=np.dot(B[i,0:n],s[0:n])
        alpha=np.dot(r,r)/np.dot(s,As)
        x=x-alpha*s
        for i in range(n):
            r[i]=np.dot(B[i,0:n],x[0:n])-b[i]
        if np.dot(r,r)<tol**2:
            break
        else:
            beta=-np.dot(r,As)/np.dot(s,As)
            s=r+beta*s
    return x,k

а вот программа которая использует этот модуль

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
__author__ = 'kardashevskiy'
import numpy as np
from scipy import sparse
from scipy.sparse import linalg
from sg1 import sg1
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
def f(x, y, r):
    return 2*c*(y-y**2+r*(x**2-x))
def poisson(n, h, r):
    V = []
    I = []
    J = []
    b = []
    for j in range(n):
        for i in range(n):
            k = i+n*j
            b.append(h**2*f(i*h+h, j*h+h, r))
            V.append(2.-2*r)
            I.append(k)
            J.append(k)
            if i>0:
                V.append(-1.)
                I.append(k)
                J.append(k-1)
            if i<n-1:
                V.append(-1.)
                I.append(k)
                J.append(k+1)
            if j>0:
                V.append(r)
                I.append(k)
                J.append(k-n)
            if j<n-1:
                V.append(r)
                I.append(k)
                J.append(k+n)
    return sparse.coo_matrix((V, (I,J))), b
N =100
h = 1./N
c = 4.
r = 0.4
tst = time.clock()
A, b = poisson(N-1, h, r)
B = A.tocsr()
v, info = sg1(B, b,tol=1.e-08)
dt = time.clock() - tst
 
if iter == 0:
    print "N = %i, time solution = %1.3e" % (N, dt)
 
fig = plt.figure()
ax = Axes3D(fig)
x = np.linspace(h/2, 1.-h/2, N-1)
y = np.linspace(h/2, 1.-h/2, N-1)
X,Y = np.meshgrid(x, y)
Z = np.reshape(v, (N-1, N-1))
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='spectral')
plt.show()

Помогите

@Wolkodav · 20.07.2014, 00:49

urung, а что пишит? Какая ошибка?

@urung · 21.07.2014, 03:02 **[ТС]**

Python
1
2
3
4
5
6
7
8
9
10
11
/usr/bin/python /home/kardashevskiy/PycharmProjects/tests/qwe1.py
Traceback (most recent call last):
  File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module>
    v, info = sg1(B, b,tol=1.e-08)
  File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1
    r[i]=np.dot(B[i,0:n],x[0:n])-b[i]
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 193, in __sub__
    raise NotImplementedError('adding a scalar to a sparse '
NotImplementedError: adding a scalar to a sparse matrix is not supported
 
Process finished with exit code 1

Добавлено через 3 минуты
Спасибо что обратили внимание извините за задержку ответа у меня часовой пояс +6

@urung · 22.07.2014, 10:10 **[ТС]**

Я боюсь что у вас не хватит терпения дочитать до конца но все таки рискну: Дело в том что эта программа qwe1прекрасно работает с библиотечным модулем linalg.cg и этот модуль sg1 в другой тестовой программе также работает более того оба модуля sg1 и linalg.cg одинаково хорошо вычисляют в той другой тестовой программе. Я же хочу эту задачу решать применяя другие модули которые конечно же похожи на этот модуль. Я приведу ту правильную которая работает. Может быть это поможет. Я не думал что это так трудно.

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
__author__ = 'kardashevskiy'
#Hyperbolic equation
import numpy as np
from scipy import sparse
from scipy.sparse import linalg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
def f(x, y, r):
    return 2*c*(y-y**2+r*(x**2-x))
def poisson(n, h, r):
    V = []
    I = []
    J = []
    b = []
    for j in range(n):
        for i in range(n):
            k = i+n*j
            b.append(h**2*f(i*h+h, j*h+h, r))
            V.append(2.-2*r)
            I.append(k)
            J.append(k)
            if i>0:
                V.append(-1.)
                I.append(k)
                J.append(k-1)
            if i<n-1:
                V.append(-1.)
                I.append(k)
                J.append(k+1)
            if j>0:
                V.append(r)
                I.append(k)
                J.append(k-n)
            if j<n-1:
                V.append(r)
                I.append(k)
                J.append(k+n)
    return sparse.coo_matrix((V, (I,J))), b
 
N =100
h = 1./N
c = 4.
r = 0.4
tst = time.clock()
A, b = poisson(N-1, h, r)
B = A.tocsr()
v, info = linalg.cg(B, b, tol=1.e-08)
dt = time.clock() - tst
 
if iter == 0:
    print "N = %i, time solution = %1.3e" % (N, dt)
 
fig = plt.figure()
ax = Axes3D(fig)
x = np.linspace(h/2, 1.-h/2, N-1)
y = np.linspace(h/2, 1.-h/2, N-1)
X,Y = np.meshgrid(x, y)
Z = np.reshape(v, (N-1, N-1))
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='spectral')
plt.show()

@Wolkodav · 22.07.2014, 10:12

Не по теме:

urung, это возможно и не трудно, просто времени пока не хватает

@urung · 24.07.2014, 03:39 **[ТС]**

Может посоветуете как найти исходный код того самого библиотечного модуля linalg.cg который я не сумел найти

@urung · 01.08.2014, 04:50 **[ТС]**

Мне помогли на другом форуме на исходниках оказывается на 47 строке конвертируем матрицу обратно в простую то есть надо не В=А.tocsr() а В=А.toarray() и все. Спасибо всем.

Новые блоги и статьи Все статьи Все блоги /
Символические и жёсткие ссылки в Linux. algri14 15.03.2026 Существует два типа ссылок — символические и жёсткие. Ссылка в Linux — это дополнительная запись в каталоге, которая может указывать либо на inode «файла-ИСТОЧНИКА», тогда это будет «жёсткая. . .	[Owen Logic] Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ФедосеевПавел 14.03.2026 Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ВВЕДЕНИЕ Выполняя задание на управление насосной группой заполнения резервуара,. . .	делаю науч статью по влиянию грибов на сукцессию anaschu 13.03.2026 прикрепляю статью	SDL3 для Desktop (MinGW): Создаём пустое окно с нуля для 2D-графики на SDL3, Си и C++ 8Observer8 10.03.2026 Содержание блога Финальные проекты на Си и на C++: hello-sdl3-c. zip hello-sdl3-cpp. zip Результат:
Установка CMake и MinGW 13.1 для сборки С и C++ приложений из консоли и из Qt Creator в EXE 8Observer8 10.03.2026 Содержание блога MinGW - это коллекция инструментов для сборки приложений в EXE. CMake - это система сборки приложений. Здесь описаны базовые шаги для старта программирования с помощью CMake и. . .	Как дизайн сайта влияет на конверсию: 7 решений, которые реально повышают заявки Neotwalker 08.03.2026 Многие до сих пор воспринимают дизайн сайта как “красивую оболочку”. На практике всё иначе: дизайн напрямую влияет на то, оставит человек заявку или уйдёт через несколько секунд. Даже если у вас. . .	Модульная разработка через nuget packages DevAlt 07.03.2026 Сложившийся в . Net-среде способ разработки чаще всего предполагает монорепозиторий в котором находятся все исходники. При создании нового решения, мы просто добавляем нужные проекты и имеем. . .	Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .

@Wolkodav 842 / 480 / 58 Регистрация: 18.09.2012 Сообщений: 1,688
	20.07.2014, 00:49
	urung, а что пишит? Какая ошибка? 0

@Wolkodav
	22.07.2014, 10:12
	Не по теме: urung, это возможно и не трудно, просто времени пока не хватает 0

@urung 0 / 0 / 0 Регистрация: 19.07.2014 Сообщений: 11
	24.07.2014, 03:39 [ТС]
	Может посоветуете как найти исходный код того самого библиотечного модуля linalg.cg который я не сумел найти 0

@urung 0 / 0 / 0 Регистрация: 19.07.2014 Сообщений: 11
	01.08.2014, 04:50 [ТС]
	Мне помогли на другом форуме на исходниках оказывается на 47 строке конвертируем матрицу обратно в простую то есть надо не В=А.tocsr() а В=А.toarray() и все. Спасибо всем. 0