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
|
from __future__ import division, print_function, absolute_import
from scipy import *
from iterative import *
def test_fun(alpha, x, beta, y, A, n):
# compute z = alpha*A*x + beta*y
xx = x[:n]
yy = y[:n]
w = dot(A,xx)
z = alpha*w+beta*yy
y[:n] = z
return
def test_fun_t(alpha, x, beta, y, A, n):
# compute z = alpha*A*x + beta*y
xx = x[:n]
yy = y[:n]
AA = conj(transpose(A))
w = dot(AA,xx)
z = alpha*w+beta*yy
y[:n] = z
return
def test_psolve(x,b,n):
x[:n] = b[:n]
return
def test_psolve_t(x,b,n):
x[:n] = b[:n]
return
def test_psolveq(x,b,which,n):
x[:n] = b[:n]
return
def test_psolveq_t(x,b,which,n):
x[:n] = b[:n]
return
n = 5
dA = 1.0*array([[2, -1, 0, 0, 0],
[-1, 2, -1, 0, 0],
[0, -1, 2, -1, 0],
[0, 0, -1, 2, -1],
[0, 0, 0, 1, 2]])
db = 1.0*array([0,1,1,0,0])
##zA = (1.0+0j)*array([[ 2, -1+0.1j, 0, 0, 0],
## [-1+0.1j, 2, -1-0.1j, 0, 0],
## [ 0, -1-0.1j, 2, -1+0.1j, 0],
## [ 0, 0, -1+0.1j, 2, -1-0.1j],
## [ 0, 0, 0, -1, 2-0.1j]])
zA = (1.0+0j)*array([[2, -1 + 1j, 0, 0, 0],
[-1+0.1j, 2, -1-0.1j, 0, 0],
[0, -1 - 1j, 2, -1+0.1j, 0],
[0, 0, -1+0.1j, 2, -1-0.1j],
[0, 0, 0, -1, 2-0.1j]])
zb = (1.0+0j)*array([0,1,1,0,0])
dx = 0*db.copy()
zx = 0*zb.copy()
diter = 1000
dresid = 1e-6
ziter = 1000
zresid = 1e-6
drestrt = n
zrestrt = n
############### BiCG #######################
dx,diter,dresid,dinfor = dbicg(db,dx,diter,dresid,test_fun,test_fun_t,test_psolve,test_psolve_t,(dA,n),(dA,n),(n,),(n,))
zx,ziter,zresid,zinfor = zbicg(zb,zx,ziter,zresid,test_fun,test_fun_t,test_psolve,test_psolve_t,(zA,n),(zA,n),(n,),(n,))
############### BiCGSTAB ###################
#dx,diter,dresid,dinfor = dbicgstab(db,dx,diter,dresid,test_fun,test_psolve,(dA,n),(n,))
#zx,ziter,zresid,zinfor = zbicgstab(zb,zx,ziter,zresid,test_fun,test_psolve,(zA,n),(n,))
############### CG #########################
##dA = 1.0*array([[ 2, -1, 0, 0, 0],
## [-1, 2, -1, 0, 0],
## [ 0, -1, 2, -1, 0],
## [ 0, 0, -1, 2, -1],
## [ 0, 0, 0, -1, 2]])
##dx = db.copy()
##zA = (1.0+0j)*array([[ 2, -1+0.1j, 0, 0, 0],
## [-1+0.1j, 2, -1-0.1j, 0, 0],
## [ 0, -1-0.1j, 2, -1+0.1j, 0],
## [ 0, 0, -1+0.1j, 2, -1-0.1j],
## [ 0, 0, 0, -1, 2-0.1j]])
##zx = zb.copy()
##dx,diter,dresid,dinfor = dcg(db,dx,diter,dresid,test_fun,test_psolve,(dA,n),(n,))
##zx,ziter,zresid,zinfor = zcg(zb,zx,ziter,zresid,test_fun,test_psolve,(zA,n),(n,))
############### CGS ########################
#dx,diter,dresid,dinfor = dcgs(db,dx,diter,dresid,test_fun,test_psolve,(dA,n),(n,))
#zx,ziter,zresid,zinfor = zcgs(zb,zx,ziter,zresid,test_fun,test_psolve,(zA,n),(n,))
############### GMRES ######################
#dx,diter,dresid,dinfor = dgmres(db,dx,drestrt,diter,dresid,test_fun,test_psolve,(dA,n),(n,))
#zx,ziter,zresid,zinfor = zgmres(zb,zx,zrestrt,ziter,zresid,test_fun,test_psolve,(zA,n),(n,))
############### QMR ########################
#dx,diter,dresid,dinfor = dqmr(db,dx,diter,dresid,test_fun,test_fun_t,test_psolveq,test_psolveq_t,(dA,n),(dA,n),(n,),(n,))
#zx,ziter,zresid,zinfor = zqmr(zb,zx,ziter,zresid,test_fun,test_fun_t,test_psolveq,test_psolveq_t,(zA,n),(zA,n),(n,),(n,))
print()
print('**************** double *****************')
print('iter:',diter, 'resid:', dresid, 'info:',dinfor)
print('x=',dx)
print('*****************************************')
print()
print()
print('**************** complex ****************')
print('iter:',ziter, 'resid:',zresid, 'info:',zinfor)
print('x=',zx)
print('*****************************************')
print()
|
