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
62
63
64
65
66
67
68
69
|
from scipy import *
import multigridtools
import scipy
import numpy
from utils import diag_sparse,inf_norm
def rs_strong_connections(A,theta):
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
Sp,Sj,Sx = multigridtools.rs_strong_connections(A.shape[0],theta,A.indptr,A.indices,A.data)
return scipy.sparse.csr_matrix((Sx,Sj,Sp),A.shape)
def rs_interpolation(A,theta=0.25):
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
S = rs_strong_connections(A,theta)
T = S.T.tocsr()
Ip,Ij,Ix = multigridtools.rs_interpolation(A.shape[0],\
A.indptr,A.indices,A.data,\
S.indptr,S.indices,S.data,\
T.indptr,T.indices,T.data)
return scipy.sparse.csr_matrix((Ix,Ij,Ip))
def sa_strong_connections(A,epsilon):
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
return scipy.sparse.csr_matrix((Sx,Sj,Sp),A.shape)
def sa_constant_interpolation(A,epsilon=None):
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
if epsilon is not None:
S = sa_strong_connections(A,epsilon)
else:
S = A
#tentative (non-smooth) interpolation operator I
Ij = multigridtools.sa_get_aggregates(A.shape[0],S.indptr,S.indices)
Ip = numpy.arange(len(Ij)+1)
Ix = numpy.ones(len(Ij))
return scipy.sparse.csr_matrix((Ix,Ij,Ip))
def sa_interpolation(A,epsilon,omega=4.0/3.0):
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
I = sa_constant_interpolation(A,epsilon)
D_inv = diag_sparse(1.0/diag_sparse(A))
D_inv_A = D_inv * A
D_inv_A *= -omega/inf_norm(D_inv_A)
P = I + (D_inv_A*I) #same as P=S*I, (faster?)
return P
|
