| 1 |
|
| 2 |
######################################################## |
| 3 |
# |
| 4 |
# Copyright (c) 2003-2008 by University of Queensland |
| 5 |
# Earth Systems Science Computational Center (ESSCC) |
| 6 |
# http://www.uq.edu.au/esscc |
| 7 |
# |
| 8 |
# Primary Business: Queensland, Australia |
| 9 |
# Licensed under the Open Software License version 3.0 |
| 10 |
# http://www.opensource.org/licenses/osl-3.0.php |
| 11 |
# |
| 12 |
######################################################## |
| 13 |
|
| 14 |
__copyright__="""Copyright (c) 2003-2008 by University of Queensland |
| 15 |
Earth Systems Science Computational Center (ESSCC) |
| 16 |
http://www.uq.edu.au/esscc |
| 17 |
Primary Business: Queensland, Australia""" |
| 18 |
__license__="""Licensed under the Open Software License version 3.0 |
| 19 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 20 |
__url__="http://www.uq.edu.au/esscc/escript-finley" |
| 21 |
|
| 22 |
""" |
| 23 |
Test suite for the pdetools module |
| 24 |
|
| 25 |
The tests must be linked with a Domain class object in the setUp method: |
| 26 |
|
| 27 |
from esys.finley import Rectangle |
| 28 |
class Test_LinearPDEOnFinley(Test_LinearPDE): |
| 29 |
RES_TOL=1.e-8 |
| 30 |
def setUp(self): |
| 31 |
self.domain = Rectangle(10,10,2) |
| 32 |
def tearDown(self): |
| 33 |
del self.domain |
| 34 |
suite = unittest.TestSuite() |
| 35 |
suite.addTest(unittest.makeSuite(Test_LinearPDEOnFinley)) |
| 36 |
unittest.TextTestRunner(verbosity=2).run(suite) |
| 37 |
|
| 38 |
@var __author__: name of author |
| 39 |
@var __copyright__: copyrights |
| 40 |
@var __license__: licence agreement |
| 41 |
@var __url__: url entry point on documentation |
| 42 |
@var __version__: version |
| 43 |
@var __date__: date of the version |
| 44 |
""" |
| 45 |
|
| 46 |
__author__="Lutz Gross, l.gross@uq.edu.au" |
| 47 |
|
| 48 |
import unittest |
| 49 |
from esys.escript import * |
| 50 |
from esys.escript.pdetools import Locator,Projector,TimeIntegrationManager,NoPDE,PCG, ArithmeticTuple, GMRES, MINRES, TFQMR, HomogeneousSaddlePointProblem |
| 51 |
from esys.escript.pdetools import Defect, NewtonGMRES |
| 52 |
from numarray.linear_algebra import solve_linear_equations |
| 53 |
|
| 54 |
class Test_pdetools_noLumping(unittest.TestCase): |
| 55 |
DEBUG=False |
| 56 |
VERBOSE=False |
| 57 |
def test_TimeIntegrationManager_scalar(self): |
| 58 |
t=0. |
| 59 |
dt=0.1 |
| 60 |
tm=TimeIntegrationManager(0.,p=1) |
| 61 |
while t<1.: |
| 62 |
t+=dt |
| 63 |
tm.checkin(dt,t) |
| 64 |
v_guess=tm.extrapolate(dt) |
| 65 |
self.failUnless(abs(v_guess-(tm.getTime()+dt))<self.RES_TOL,"extrapolation is wrong") |
| 66 |
|
| 67 |
def test_TimeIntegrationManager_vector(self): |
| 68 |
t=0. |
| 69 |
dt=0.3 |
| 70 |
tm=TimeIntegrationManager(0.,0.,p=1) |
| 71 |
while t<1.: |
| 72 |
t+=dt |
| 73 |
tm.checkin(dt,t,3*t) |
| 74 |
v_guess=tm.extrapolate(dt) |
| 75 |
e=max(abs(v_guess[0]-(tm.getTime()+dt)),abs(v_guess[1]-(tm.getTime()+dt)*3.)) |
| 76 |
self.failUnless(e<self.RES_TOL,"extrapolation is wrong") |
| 77 |
|
| 78 |
def test_Locator(self): |
| 79 |
x=self.domain.getX() |
| 80 |
l=Locator(self.domain,numarray.ones((self.domain.getDim(),))) |
| 81 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain") |
| 82 |
|
| 83 |
l=Locator(ContinuousFunction(self.domain),numarray.ones((self.domain.getDim(),))) |
| 84 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space") |
| 85 |
|
| 86 |
xx=l.getX() |
| 87 |
self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type") |
| 88 |
self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"location wrong") |
| 89 |
xx=l(x) |
| 90 |
self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type") |
| 91 |
self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"value wrong vector") |
| 92 |
xx=l(x[0]+x[1]) |
| 93 |
self.failUnless(isinstance(xx,float),"wrong scalar type") |
| 94 |
self.failUnless(abs(xx-2.)<self.RES_TOL,"value wrong scalar") |
| 95 |
|
| 96 |
def test_Locator_withList(self): |
| 97 |
x=self.domain.getX() |
| 98 |
arg=[numarray.ones((self.domain.getDim(),)), numarray.zeros((self.domain.getDim(),))] |
| 99 |
l=Locator(self.domain,arg) |
| 100 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain") |
| 101 |
|
| 102 |
l=Locator(ContinuousFunction(self.domain),arg) |
| 103 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space") |
| 104 |
|
| 105 |
xx=l.getX() |
| 106 |
self.failUnless(isinstance(xx,list),"list expected") |
| 107 |
for i in range(len(xx)): |
| 108 |
self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item"%i) |
| 109 |
self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong"%i) |
| 110 |
xx=l(x) |
| 111 |
self.failUnless(isinstance(xx,list),"list expected (2)") |
| 112 |
for i in range(len(xx)): |
| 113 |
self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item (2)"%i) |
| 114 |
self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong (2)"%i) |
| 115 |
xx=l(x[0]+x[1]) |
| 116 |
self.failUnless(isinstance(xx,list),"list expected (3)") |
| 117 |
for i in range(len(xx)): |
| 118 |
self.failUnless(isinstance(xx[i],float),"wrong scalar type") |
| 119 |
self.failUnless(abs(xx[i]-(arg[i][0]+arg[i][1]))<self.RES_TOL,"value wrong scalar") |
| 120 |
|
| 121 |
def testProjector_rank0(self): |
| 122 |
x=ContinuousFunction(self.domain).getX() |
| 123 |
p=Projector(self.domain,reduce=False,fast=False) |
| 124 |
td_ref=x[0] |
| 125 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 126 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 127 |
|
| 128 |
def testProjector_rank1(self): |
| 129 |
x=ContinuousFunction(self.domain).getX() |
| 130 |
p=Projector(self.domain,reduce=False,fast=False) |
| 131 |
td_ref=x |
| 132 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 133 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 134 |
|
| 135 |
def testProjector_rank2(self): |
| 136 |
x=ContinuousFunction(self.domain).getX() |
| 137 |
p=Projector(self.domain,reduce=False,fast=False) |
| 138 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 139 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 140 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 141 |
|
| 142 |
def testProjector_rank3(self): |
| 143 |
x=ContinuousFunction(self.domain).getX() |
| 144 |
p=Projector(self.domain,reduce=False,fast=False) |
| 145 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 146 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 147 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 148 |
|
| 149 |
def testProjector_rank4(self): |
| 150 |
x=ContinuousFunction(self.domain).getX() |
| 151 |
p=Projector(self.domain,reduce=False,fast=False) |
| 152 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 153 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 154 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 155 |
|
| 156 |
|
| 157 |
def testProjector_rank0_reduced(self): |
| 158 |
x=ContinuousFunction(self.domain).getX() |
| 159 |
p=Projector(self.domain,reduce=True,fast=False) |
| 160 |
td_ref=x[0] |
| 161 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 162 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 163 |
|
| 164 |
def testProjector_rank1_reduced(self): |
| 165 |
x=ContinuousFunction(self.domain).getX() |
| 166 |
p=Projector(self.domain,reduce=True,fast=False) |
| 167 |
td_ref=x |
| 168 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 169 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 170 |
|
| 171 |
def testProjector_rank2_reduced(self): |
| 172 |
x=ContinuousFunction(self.domain).getX() |
| 173 |
p=Projector(self.domain,reduce=True,fast=False) |
| 174 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 175 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 176 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 177 |
|
| 178 |
def testProjector_rank3_reduced(self): |
| 179 |
x=ContinuousFunction(self.domain).getX() |
| 180 |
p=Projector(self.domain,reduce=True,fast=False) |
| 181 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 182 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 183 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 184 |
|
| 185 |
def testProjector_rank4_reduced(self): |
| 186 |
x=ContinuousFunction(self.domain).getX() |
| 187 |
p=Projector(self.domain,reduce=True,fast=False) |
| 188 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 189 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 190 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 191 |
|
| 192 |
def testProjector_rank0_with_reduced_input(self): |
| 193 |
x=ContinuousFunction(self.domain).getX() |
| 194 |
p=Projector(self.domain,reduce=False,fast=False) |
| 195 |
td_ref=x[0] |
| 196 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 197 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 198 |
|
| 199 |
def testProjector_rank1_with_reduced_input(self): |
| 200 |
x=ContinuousFunction(self.domain).getX() |
| 201 |
p=Projector(self.domain,reduce=False,fast=False) |
| 202 |
td_ref=x |
| 203 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 204 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 205 |
|
| 206 |
def testProjector_rank2_with_reduced_input(self): |
| 207 |
x=ContinuousFunction(self.domain).getX() |
| 208 |
p=Projector(self.domain,reduce=False,fast=False) |
| 209 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 210 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 211 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 212 |
|
| 213 |
def testProjector_rank3_with_reduced_input(self): |
| 214 |
x=ContinuousFunction(self.domain).getX() |
| 215 |
p=Projector(self.domain,reduce=False,fast=False) |
| 216 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 217 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 218 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 219 |
|
| 220 |
def testProjector_rank4_with_reduced_input(self): |
| 221 |
x=ContinuousFunction(self.domain).getX() |
| 222 |
p=Projector(self.domain,reduce=False,fast=False) |
| 223 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 224 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 225 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 226 |
|
| 227 |
|
| 228 |
def testProjector_rank0_reduced_with_reduced_input(self): |
| 229 |
x=ContinuousFunction(self.domain).getX() |
| 230 |
p=Projector(self.domain,reduce=True,fast=False) |
| 231 |
td_ref=1. |
| 232 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 233 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 234 |
|
| 235 |
def testProjector_rank1_reduced_with_reduced_input(self): |
| 236 |
x=ContinuousFunction(self.domain).getX() |
| 237 |
p=Projector(self.domain,reduce=True,fast=False) |
| 238 |
td_ref=numarray.array([1.,2.,3.]) |
| 239 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 240 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 241 |
|
| 242 |
def testProjector_rank2_reduced_with_reduced_input(self): |
| 243 |
x=ContinuousFunction(self.domain).getX() |
| 244 |
p=Projector(self.domain,reduce=True,fast=False) |
| 245 |
td_ref=numarray.array([[11.,12.],[21,22.]]) |
| 246 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 247 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 248 |
|
| 249 |
def testProjector_rank3_reduced_with_reduced_input(self): |
| 250 |
x=ContinuousFunction(self.domain).getX() |
| 251 |
p=Projector(self.domain,reduce=True,fast=False) |
| 252 |
td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]) |
| 253 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 254 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 255 |
|
| 256 |
def testProjector_rank4_reduced_with_reduced_input(self): |
| 257 |
x=ContinuousFunction(self.domain).getX() |
| 258 |
p=Projector(self.domain,reduce=True,fast=False) |
| 259 |
td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]) |
| 260 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 261 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 262 |
|
| 263 |
|
| 264 |
def test_NoPDE_scalar_missing_r(self): |
| 265 |
p=NoPDE(self.domain) |
| 266 |
x=self.domain.getX() |
| 267 |
msk=whereZero(x[0]) |
| 268 |
p.setValue(D=1.,Y=1.,q=msk) |
| 269 |
u=p.getSolution() |
| 270 |
u_ex=(1.-msk) |
| 271 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 272 |
|
| 273 |
def test_NoPDE_scalar_missing_Y(self): |
| 274 |
p=NoPDE(self.domain) |
| 275 |
x=self.domain.getX() |
| 276 |
msk=whereZero(x[0]) |
| 277 |
p.setValue(D=1.,q=msk,r=2.) |
| 278 |
u=p.getSolution() |
| 279 |
u_ex=msk*2. |
| 280 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 281 |
|
| 282 |
def test_NoPDE_scalar_constant(self): |
| 283 |
p=NoPDE(self.domain) |
| 284 |
x=self.domain.getX() |
| 285 |
msk=whereZero(x[0]) |
| 286 |
p.setValue(D=1.,Y=1.,q=msk,r=2.) |
| 287 |
u=p.getSolution() |
| 288 |
u_ex=(1.-msk)+msk*2. |
| 289 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 290 |
|
| 291 |
def test_NoPDE_scalar_variable(self): |
| 292 |
p=NoPDE(self.domain) |
| 293 |
x=self.domain.getX() |
| 294 |
msk=whereZero(x[0]) |
| 295 |
p.setValue(D=x[0],Y=2*x[0],q=msk,r=2.) |
| 296 |
u=p.getSolution() |
| 297 |
u_ex=2. |
| 298 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 299 |
|
| 300 |
def test_NoPDE_vector_missing_Y(self): |
| 301 |
p=NoPDE(self.domain) |
| 302 |
x=self.domain.getX() |
| 303 |
msk=whereZero(x[0])*[1.,0.] |
| 304 |
p.setValue(D=numarray.ones([2]),q=msk,r=2.) |
| 305 |
u=p.getSolution() |
| 306 |
u_ex=msk*2. |
| 307 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 308 |
|
| 309 |
def test_NoPDE_vector_missing_r(self): |
| 310 |
p=NoPDE(self.domain) |
| 311 |
x=self.domain.getX() |
| 312 |
msk=whereZero(x[0])*[1.,0.] |
| 313 |
p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk) |
| 314 |
u=p.getSolution() |
| 315 |
u_ex=(1.-msk) |
| 316 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 317 |
|
| 318 |
def test_NoPDE_vector_constant(self): |
| 319 |
p=NoPDE(self.domain) |
| 320 |
x=self.domain.getX() |
| 321 |
msk=whereZero(x[0])*[1.,0.] |
| 322 |
p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk,r=2.) |
| 323 |
u=p.getSolution() |
| 324 |
u_ex=(1.-msk)+msk*2. |
| 325 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 326 |
|
| 327 |
def test_NoPDE_vector_variable(self): |
| 328 |
p=NoPDE(self.domain) |
| 329 |
x=self.domain.getX() |
| 330 |
msk=whereZero(x[0])*[1.,0.] |
| 331 |
p.setValue(D=x[:2]+1,Y=2*(x[:2]+1),q=msk,r=2.) |
| 332 |
u=p.getSolution() |
| 333 |
u_ex=2. |
| 334 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 335 |
#===== |
| 336 |
def testPCG(self): |
| 337 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 338 |
from math import sqrt |
| 339 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 340 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 341 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 342 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 343 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 344 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 345 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 346 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 347 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 348 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 349 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 350 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 351 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 352 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 353 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 354 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 355 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 356 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 357 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 358 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 359 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 360 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 361 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 362 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 363 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 364 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 365 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 366 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 367 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 368 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 369 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 370 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 371 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 372 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 373 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 374 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 375 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 376 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 377 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 378 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 379 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 380 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 381 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 382 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 383 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 384 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 385 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 386 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 387 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 388 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 389 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 390 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 391 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 392 |
0.314596738052894]) |
| 393 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 394 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 395 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 396 |
-8.7934289814322 ]) |
| 397 |
|
| 398 |
def Ap(x): |
| 399 |
return matrixmultiply(A,x) |
| 400 |
def Ms(b): |
| 401 |
out=zeros((size(b),),Float64) |
| 402 |
for i in xrange(size(b)): |
| 403 |
out[i]=b[i]/A[i,i] |
| 404 |
return out |
| 405 |
|
| 406 |
tol=1.e-4 |
| 407 |
x,r,a_norm=PCG(b*1.,Ap,x_ref*0.,Ms,dot, atol=0, rtol=tol, iter_max=12) |
| 408 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 409 |
self.failUnless(Lsup(r-(b-matrixmultiply(A,x)))<=Lsup(b)*EPSILON*100.,"wrong solution") |
| 410 |
|
| 411 |
def testMINRES(self): |
| 412 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 413 |
from math import sqrt |
| 414 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 415 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 416 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 417 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 418 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 419 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 420 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 421 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 422 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 423 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 424 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 425 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 426 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 427 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 428 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 429 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 430 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 431 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 432 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 433 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 434 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 435 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 436 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 437 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 438 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 439 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 440 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 441 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 442 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 443 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 444 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 445 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 446 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 447 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 448 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 449 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 450 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 451 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 452 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 453 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 454 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 455 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 456 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 457 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 458 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 459 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 460 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 461 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 462 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 463 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 464 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 465 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 466 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 467 |
0.314596738052894]) |
| 468 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 469 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 470 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 471 |
-8.7934289814322 ]) |
| 472 |
|
| 473 |
def Ap(x): |
| 474 |
return matrixmultiply(A,x) |
| 475 |
def Ms(b): |
| 476 |
out=zeros((size(b),),Float64) |
| 477 |
for i in xrange(size(b)): |
| 478 |
out[i]=b[i]/A[i,i] |
| 479 |
return out |
| 480 |
|
| 481 |
tol=1.e-4 |
| 482 |
x=MINRES(b*1.,Ap,x_ref*0,Ms,dot, atol=0, rtol=tol, iter_max=12) |
| 483 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 484 |
|
| 485 |
def testTFQMR(self): |
| 486 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 487 |
from math import sqrt |
| 488 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 489 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 490 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 491 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 492 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 493 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 494 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 495 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 496 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 497 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 498 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 499 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 500 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 501 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 502 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 503 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 504 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 505 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 506 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 507 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 508 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 509 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 510 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 511 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 512 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 513 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 514 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 515 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 516 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 517 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 518 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 519 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 520 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 521 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 522 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 523 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 524 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 525 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 526 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 527 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 528 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 529 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 530 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 531 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 532 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 533 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 534 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 535 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 536 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 537 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 538 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 539 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 540 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 541 |
0.314596738052894]) |
| 542 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 543 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 544 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 545 |
-8.7934289814322 ]) |
| 546 |
|
| 547 |
def Ap(x): |
| 548 |
out=matrixmultiply(A,x) |
| 549 |
for i in xrange(size(x)): |
| 550 |
out[i]/=A[i,i] |
| 551 |
return out |
| 552 |
|
| 553 |
tol=1.e-5 |
| 554 |
for i in xrange(size(b)): b[i]/=A[i,i] |
| 555 |
x=TFQMR(b,Ap,x_ref*0,dot, atol=0, rtol=tol, iter_max=12) |
| 556 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 557 |
|
| 558 |
def testGMRES(self): |
| 559 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 560 |
from math import sqrt |
| 561 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 562 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 563 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 564 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 565 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 566 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 567 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 568 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 569 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 570 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 571 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 572 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 573 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 574 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 575 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 576 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 577 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 578 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 579 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 580 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 581 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 582 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 583 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 584 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 585 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 586 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 587 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 588 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 589 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 590 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 591 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 592 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 593 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 594 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 595 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 596 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 597 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 598 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 599 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 600 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 601 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 602 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 603 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 604 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 605 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 606 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 607 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 608 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 609 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 610 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 611 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 612 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 613 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 614 |
0.314596738052894]) |
| 615 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 616 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 617 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 618 |
-8.7934289814322 ]) |
| 619 |
|
| 620 |
def Ap(x): |
| 621 |
b=matrixmultiply(A,x) |
| 622 |
for i in xrange(size(b)): |
| 623 |
b[i]/=A[i,i] |
| 624 |
return b |
| 625 |
|
| 626 |
tol=1.e-4 |
| 627 |
for i in xrange(size(b)): b[i]/=A[i,i] |
| 628 |
x=GMRES(b,Ap,x_ref*0,dot,atol=0, rtol=tol, iter_max=12) |
| 629 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 630 |
|
| 631 |
def testGMRES_P_R(self): |
| 632 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 633 |
from math import sqrt |
| 634 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 635 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 636 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 637 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 638 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 639 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 640 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 641 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 642 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 643 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 644 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 645 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 646 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 647 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 648 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 649 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 650 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 651 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 652 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 653 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 654 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 655 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 656 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 657 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 658 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 659 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 660 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 661 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 662 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 663 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 664 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 665 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 666 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 667 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 668 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 669 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 670 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 671 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 672 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 673 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 674 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 675 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 676 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 677 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 678 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 679 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 680 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 681 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 682 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 683 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 684 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 685 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 686 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 687 |
0.314596738052894]) |
| 688 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 689 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 690 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 691 |
-8.7934289814322 ]) |
| 692 |
|
| 693 |
def Ap(x): |
| 694 |
return matrixmultiply(A,x) |
| 695 |
def P_Rp(x): |
| 696 |
out=zeros(size(x)) |
| 697 |
for i in xrange(size(b)): |
| 698 |
out[i]=x[i]/A[i,i] |
| 699 |
return out |
| 700 |
|
| 701 |
tol=1.e-4 |
| 702 |
x=GMRES(b,Ap,x_ref*0,dot,atol=0, rtol=tol, iter_max=12,P_R=P_Rp) |
| 703 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 704 |
|
| 705 |
def testNewtonGMRES(self): |
| 706 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 707 |
from math import sqrt |
| 708 |
class LL(Defect): |
| 709 |
def __init__(self,*kwargs): |
| 710 |
super(LL, self).__init__(*kwargs) |
| 711 |
self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 712 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 713 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 714 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 715 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 716 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 717 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 718 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 719 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 720 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 721 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 722 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 723 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 724 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 725 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 726 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 727 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 728 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 729 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 730 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 731 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 732 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 733 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 734 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 735 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 736 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 737 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 738 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 739 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 740 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 741 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 742 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 743 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 744 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 745 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 746 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 747 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 748 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 749 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 750 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 751 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 752 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 753 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 754 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 755 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 756 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 757 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 758 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 759 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 760 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 761 |
self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 762 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 763 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 764 |
0.314596738052894]) |
| 765 |
self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 766 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 767 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 768 |
-8.7934289814322 ]) |
| 769 |
def eval(self,x): |
| 770 |
out=matrixmultiply(self.A,x)-self.b |
| 771 |
for i in xrange(size(self.b)): |
| 772 |
out[i]/=self.A[i,i] |
| 773 |
return out |
| 774 |
def bilinearform(self,x0,x1): |
| 775 |
return dot(x0,x1) |
| 776 |
|
| 777 |
tol=1.e-8 |
| 778 |
ll=LL() |
| 779 |
x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE) |
| 780 |
self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution") |
| 781 |
|
| 782 |
def testNewtonGMRES(self): |
| 783 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 784 |
from math import sqrt |
| 785 |
class LL(Defect): |
| 786 |
def __init__(self,*kwargs): |
| 787 |
super(LL, self).__init__(*kwargs) |
| 788 |
self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 789 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 790 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 791 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 792 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 793 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 794 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 795 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 796 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 797 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 798 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 799 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 800 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 801 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 802 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 803 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 804 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 805 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 806 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 807 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 808 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 809 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 810 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 811 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 812 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 813 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 814 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 815 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 816 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 817 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 818 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 819 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 820 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 821 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 822 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 823 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 824 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 825 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 826 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 827 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 828 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 829 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 830 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 831 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 832 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 833 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 834 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 835 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 836 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 837 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 838 |
self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 839 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 840 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 841 |
0.314596738052894]) |
| 842 |
self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 843 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 844 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 845 |
-8.7934289814322 ]) |
| 846 |
def eval(self,x): |
| 847 |
out=matrixmultiply(self.A,x)-self.b |
| 848 |
for i in xrange(size(self.b)): |
| 849 |
out[i]/=self.A[i,i] |
| 850 |
return out |
| 851 |
def bilinearform(self,x0,x1): |
| 852 |
return dot(x0,x1) |
| 853 |
|
| 854 |
tol=1.e-8 |
| 855 |
ll=LL() |
| 856 |
x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE) |
| 857 |
self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution") |
| 858 |
|
| 859 |
def testHomogeneousSaddlePointProblem_PCG(self): |
| 860 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 861 |
from math import sqrt |
| 862 |
class LL(HomogeneousSaddlePointProblem): |
| 863 |
def initialize(self): |
| 864 |
self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 865 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 866 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 867 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 868 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 869 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 870 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 871 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 872 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 873 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 874 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 875 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 876 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 877 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 878 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 879 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 880 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 881 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 882 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 883 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 884 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 885 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 886 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 887 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 888 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 889 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 890 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 891 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 892 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 893 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 894 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 895 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 896 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 897 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 898 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 899 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 900 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 901 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 902 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 903 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 904 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 905 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 906 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 907 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 908 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 909 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 910 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 911 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 912 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 913 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 914 |
self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721, |
| 915 |
0.253012436539738, -0.346223308561905, 0.2425508275422, |
| 916 |
-0.194695862196008, 0.09451439391473, 0.302961126826511, |
| 917 |
-0.236043777597633] ) |
| 918 |
|
| 919 |
self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614], |
| 920 |
[ 0.062879587145773 ,0.038798770300146, 0.022155850155616], |
| 921 |
[ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ], |
| 922 |
[ 0.059000597728388 ,0.090986953740106, 0.035316011834982], |
| 923 |
[ 0.091209362659698 ,0.13205572801294 , 0.069462874306956], |
| 924 |
[ 0.077790176986096 ,0.133626423045765, 0.011149969846981], |
| 925 |
[ 0.01407283482513 ,0.094910926488907, 0.133498532648644], |
| 926 |
[ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ], |
| 927 |
[ 0.071254288170748 ,0.071738715618163, 0.078005951991733], |
| 928 |
[ 0.049463014576779 ,0.103559223780991, 0.003356415647637]]) |
| 929 |
self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546]) |
| 930 |
|
| 931 |
self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913, |
| 932 |
6.512480714694167 , -5.727371407390975 , 4.802494840775022, |
| 933 |
-4.171606044721161 , -1.862366353566293 ,74.850226163257105, |
| 934 |
-118.602464657076439]) |
| 935 |
|
| 936 |
self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696], |
| 937 |
[-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ], |
| 938 |
[ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]]) |
| 939 |
def inner_pBv(self,p,v): |
| 940 |
return dot(p,matrixmultiply(transpose(self.Bt),v)) |
| 941 |
def inner_p(self,p0,p1): |
| 942 |
return dot(p0,p1) |
| 943 |
def norm_v(self,v): |
| 944 |
return sqrt(dot(v,v)) |
| 945 |
def getV(self,p,v): |
| 946 |
out=self.b-matrixmultiply(self.Bt,p) |
| 947 |
return solve_linear_equations(self.A,out)+v*self.getSubProblemTolerance() |
| 948 |
def norm_Bv(self,v): |
| 949 |
t=matrixmultiply(transpose(self.Bt),v) |
| 950 |
return sqrt(dot(t,t)) |
| 951 |
def solve_AinvBt(self,p): |
| 952 |
out=matrixmultiply(self.Bt,p) |
| 953 |
return solve_linear_equations(self.A,out) |
| 954 |
def solve_precB(self,v): |
| 955 |
out=matrixmultiply(transpose(self.Bt),v) |
| 956 |
for i in xrange(size(out)): out[i]*=self.Sinv[i,i] |
| 957 |
return out |
| 958 |
|
| 959 |
tol=1.e-8 |
| 960 |
ll=LL() |
| 961 |
ll.initialize() |
| 962 |
ll.setTolerance(tol) |
| 963 |
# ll.setSubToleranceReductionFactor(0.1) |
| 964 |
x,p=ll.solve(ll.x_ref*1.20,ll.p_ref*(-2),max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20,max_correction_steps=3) |
| 965 |
self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution") |
| 966 |
self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution") |
| 967 |
|
| 968 |
def testHomogeneousSaddlePointProblem_GMRES(self): |
| 969 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 970 |
from math import sqrt |
| 971 |
class LL(HomogeneousSaddlePointProblem): |
| 972 |
def initialize(self): |
| 973 |
self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 974 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 975 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 976 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 977 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 978 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 979 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 980 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 981 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 982 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 983 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 984 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 985 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 986 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 987 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 988 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 989 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 990 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 991 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 992 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 993 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 994 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 995 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 996 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 997 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 998 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 999 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 1000 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 1001 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 1002 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 1003 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 1004 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 1005 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 1006 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 1007 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 1008 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 1009 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 1010 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 1011 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 1012 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 1013 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 1014 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 1015 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 1016 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 1017 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 1018 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 1019 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 1020 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 1021 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 1022 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 1023 |
self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721, |
| 1024 |
0.253012436539738, -0.346223308561905, 0.2425508275422, |
| 1025 |
-0.194695862196008, 0.09451439391473, 0.302961126826511, |
| 1026 |
-0.236043777597633] ) |
| 1027 |
|
| 1028 |
self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614], |
| 1029 |
[ 0.062879587145773 ,0.038798770300146, 0.022155850155616], |
| 1030 |
[ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ], |
| 1031 |
[ 0.059000597728388 ,0.090986953740106, 0.035316011834982], |
| 1032 |
[ 0.091209362659698 ,0.13205572801294 , 0.069462874306956], |
| 1033 |
[ 0.077790176986096 ,0.133626423045765, 0.011149969846981], |
| 1034 |
[ 0.01407283482513 ,0.094910926488907, 0.133498532648644], |
| 1035 |
[ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ], |
| 1036 |
[ 0.071254288170748 ,0.071738715618163, 0.078005951991733], |
| 1037 |
[ 0.049463014576779 ,0.103559223780991, 0.003356415647637]]) |
| 1038 |
self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546]) |
| 1039 |
|
| 1040 |
self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913, |
| 1041 |
6.512480714694167 , -5.727371407390975 , 4.802494840775022, |
| 1042 |
-4.171606044721161 , -1.862366353566293 ,74.850226163257105, |
| 1043 |
-118.602464657076439]) |
| 1044 |
|
| 1045 |
self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696], |
| 1046 |
[-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ], |
| 1047 |
[ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]]) |
| 1048 |
def inner_pBv(self,p,v): |
| 1049 |
return dot(p,matrixmultiply(transpose(self.Bt),v)) |
| 1050 |
def inner_p(self,p0,p1): |
| 1051 |
return dot(p0,p1) |
| 1052 |
def norm_v(self,v): |
| 1053 |
return sqrt(dot(v,v)) |
| 1054 |
def getV(self,p,v): |
| 1055 |
out=self.b-matrixmultiply(self.Bt,p) |
| 1056 |
return solve_linear_equations(self.A,out)+v*self.getSubProblemTolerance() |
| 1057 |
def norm_Bv(self,v): |
| 1058 |
t=matrixmultiply(transpose(self.Bt),v) |
| 1059 |
return sqrt(dot(t,t)) |
| 1060 |
def solve_AinvBt(self,p): |
| 1061 |
out=matrixmultiply(self.Bt,p) |
| 1062 |
return solve_linear_equations(self.A,out) |
| 1063 |
def solve_precB(self,v): |
| 1064 |
out=matrixmultiply(transpose(self.Bt),v) |
| 1065 |
for i in xrange(size(out)): out[i]*=self.Sinv[i,i] |
| 1066 |
return out |
| 1067 |
|
| 1068 |
tol=1.e-8 |
| 1069 |
ll=LL() |
| 1070 |
ll.initialize() |
| 1071 |
ll.setTolerance(tol) |
| 1072 |
# ll.setSubToleranceReductionFactor(0.1) |
| 1073 |
x,p=ll.solve(ll.x_ref*1.20,ll.p_ref*(-2),max_iter=20, verbose=False, show_details=False, usePCG=False, iter_restart=20,max_correction_steps=3) |
| 1074 |
self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution") |
| 1075 |
self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution") |
| 1076 |
|
| 1077 |
def testArithmeticTuple(self): |
| 1078 |
a=ArithmeticTuple(1.,2.) |
| 1079 |
self.failUnless(len(a)==2,"wrong length") |
| 1080 |
self.failUnless(a[0]==1.,"wrong first item") |
| 1081 |
self.failUnless(a[1]==2.,"wrong second item") |
| 1082 |
c=a*6. |
| 1083 |
self.failUnless(isinstance(c,ArithmeticTuple),"c is not an instance of ArithmeticTuple") |
| 1084 |
self.failUnless(len(c)==2,"c has wrong length") |
| 1085 |
self.failUnless(c[0]==6.,"c has wrong first item") |
| 1086 |
self.failUnless(c[1]==12.,"c has wrong second item") |
| 1087 |
b=5.*a |
| 1088 |
self.failUnless(isinstance(b,ArithmeticTuple),"b is not an instance of ArithmeticTuple") |
| 1089 |
self.failUnless(len(b)==2,"b has wrong length") |
| 1090 |
self.failUnless(b[0]==5.,"b has wrong first item") |
| 1091 |
self.failUnless(b[1]==10.,"b has wrong second item") |
| 1092 |
a+=ArithmeticTuple(3.,4.) |
| 1093 |
self.failUnless(a[0]==4.,"wrong first item of inplace update") |
| 1094 |
self.failUnless(a[1]==6.,"wrong second item of inplace update") |
| 1095 |
|
| 1096 |
|
| 1097 |
|
| 1098 |
class Test_pdetools(Test_pdetools_noLumping): |
| 1099 |
def testProjector_rank0_fast_reduced(self): |
| 1100 |
x=ContinuousFunction(self.domain).getX() |
| 1101 |
h=Lsup(self.domain.getSize()) |
| 1102 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1103 |
td_ref=x[0] |
| 1104 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 1105 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1106 |
|
| 1107 |
def testProjector_rank1_fast_reduced(self): |
| 1108 |
x=ContinuousFunction(self.domain).getX() |
| 1109 |
h=Lsup(self.domain.getSize()) |
| 1110 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1111 |
td_ref=x |
| 1112 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 1113 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1114 |
|
| 1115 |
def testProjector_rank2_fast_reduced(self): |
| 1116 |
x=ContinuousFunction(self.domain).getX() |
| 1117 |
h=Lsup(self.domain.getSize()) |
| 1118 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1119 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 1120 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 1121 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1122 |
|
| 1123 |
def testProjector_rank3_fast_reduced(self): |
| 1124 |
x=ContinuousFunction(self.domain).getX() |
| 1125 |
h=Lsup(self.domain.getSize()) |
| 1126 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1127 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 1128 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 1129 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1130 |
|
| 1131 |
def testProjector_rank4_fast_reduced(self): |
| 1132 |
x=ContinuousFunction(self.domain).getX() |
| 1133 |
h=Lsup(self.domain.getSize()) |
| 1134 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1135 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 1136 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 1137 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1138 |
|
| 1139 |
def testProjector_rank0_fast_reduced_with_reduced_input(self): |
| 1140 |
x=ContinuousFunction(self.domain).getX() |
| 1141 |
h=Lsup(self.domain.getSize()) |
| 1142 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1143 |
td_ref=1. |
| 1144 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 1145 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1146 |
|
| 1147 |
def testProjector_rank1_fast_reduced_with_reduced_input(self): |
| 1148 |
x=ContinuousFunction(self.domain).getX() |
| 1149 |
h=Lsup(self.domain.getSize()) |
| 1150 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1151 |
td_ref=numarray.array([1.,2.,3.]) |
| 1152 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 1153 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1154 |
|
| 1155 |
def testProjector_rank2_fast_reduced_with_reduced_input(self): |
| 1156 |
x=ContinuousFunction(self.domain).getX() |
| 1157 |
h=Lsup(self.domain.getSize()) |
| 1158 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1159 |
td_ref=numarray.array([[11.,12.],[21,22.]]) |
| 1160 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 1161 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1162 |
|
| 1163 |
def testProjector_rank3_fast_reduced_with_reduced_input(self): |
| 1164 |
x=ContinuousFunction(self.domain).getX() |
| 1165 |
h=Lsup(self.domain.getSize()) |
| 1166 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1167 |
td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]) |
| 1168 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 1169 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1170 |
|
| 1171 |
def testProjector_rank4_fast_reduced_with_reduced_input(self): |
| 1172 |
x=ContinuousFunction(self.domain).getX() |
| 1173 |
h=Lsup(self.domain.getSize()) |
| 1174 |
p=Projector(self.domain,reduce=True,fast=True) |
| 1175 |
td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]) |
| 1176 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 1177 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 1178 |
|
Parent Directory
| 
Revision Log