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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
#!/usr/bin/env python
import numpy as np
from pyarpack import sparseBiCGILU as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex128(complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex128(complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex128(complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 200
arpackSlv.slvILUDropTol = 1.
arpackSlv.slvILUFillFactor = 2
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex64(complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex64(complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex64(complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 200
arpackSlv.slvILUDropTol = 1.
arpackSlv.slvILUFillFactor = 2
arpackSlv.sigmaReal = 1
arpackSlv.sigmaImag = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
|
