| 435 |
""" |
""" |
| 436 |
pass |
pass |
| 437 |
|
|
| 438 |
def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)): |
class IterationHistory(object): |
| 439 |
|
""" |
| 440 |
|
The IterationHistory class is used to define a stopping criterium. It keeps track of the |
| 441 |
|
residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by |
| 442 |
|
a given tolerance. |
| 443 |
|
""" |
| 444 |
|
def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False): |
| 445 |
|
""" |
| 446 |
|
Initialization |
| 447 |
|
|
| 448 |
|
@param tolerance: tolerance |
| 449 |
|
@type tolerance: positive C{float} |
| 450 |
|
@param verbose: switches on the printing out some information |
| 451 |
|
@type verbose: C{bool} |
| 452 |
|
""" |
| 453 |
|
if not tolerance>0.: |
| 454 |
|
raise ValueError,"tolerance needs to be positive." |
| 455 |
|
self.tolerance=tolerance |
| 456 |
|
self.verbose=verbose |
| 457 |
|
self.history=[] |
| 458 |
|
def stoppingcriterium(self,norm_r,r,x): |
| 459 |
|
""" |
| 460 |
|
returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]} is the residual norm at the first call. |
| 461 |
|
|
| 462 |
|
|
| 463 |
|
@param norm_r: current residual norm |
| 464 |
|
@type norm_r: non-negative C{float} |
| 465 |
|
@param r: current residual (not used) |
| 466 |
|
@param x: current solution approximation (not used) |
| 467 |
|
@return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned. |
| 468 |
|
@rtype: C{bool} |
| 469 |
|
|
| 470 |
|
""" |
| 471 |
|
self.history.append(norm_r) |
| 472 |
|
if self.verbose: print "iter: %s: inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1]) |
| 473 |
|
return self.history[-1]<=self.tolerance * self.history[0] |
| 474 |
|
|
| 475 |
|
def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100): |
| 476 |
""" |
""" |
| 477 |
Solver for |
Solver for |
| 478 |
|
|
| 481 |
with a symmetric and positive definite operator A (more details required!). |
with a symmetric and positive definite operator A (more details required!). |
| 482 |
It uses the conjugate gradient method with preconditioner M providing an approximation of A. |
It uses the conjugate gradient method with preconditioner M providing an approximation of A. |
| 483 |
|
|
| 484 |
The iteration is terminated if |
The iteration is terminated if the C{stoppingcriterium} function return C{True}. |
|
|
|
|
M{norm(r) <= tolerance * norm(b)} |
|
|
|
|
|
where C{norm()} defines a norm and |
|
|
|
|
|
M{r = b-Ax} |
|
|
|
|
|
the residual. |
|
| 485 |
|
|
| 486 |
For details on the preconditioned conjugate gradient method see the book: |
For details on the preconditioned conjugate gradient method see the book: |
| 487 |
|
|
| 490 |
C. Romine, and H. van der Vorst. |
C. Romine, and H. van der Vorst. |
| 491 |
|
|
| 492 |
@param b: the right hand side of the liner system. C{b} is altered. |
@param b: the right hand side of the liner system. C{b} is altered. |
| 493 |
@type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y) |
@type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y) |
|
@param x: an initial guess for the solution |
|
|
@type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y) |
|
| 494 |
@param Aprod: returns the value Ax |
@param Aprod: returns the value Ax |
| 495 |
@type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R. |
@type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}. |
| 496 |
@param Msolve: solves Mx=r |
@param Msolve: solves Mx=r |
| 497 |
@type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S. |
@type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same |
| 498 |
|
type like argument C{x}. |
| 499 |
@param bilinearform: inner product C{<x,r>} |
@param bilinearform: inner product C{<x,r>} |
| 500 |
@type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}. |
@type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}. |
| 501 |
@param norm: norm calculation for the residual C{r=b-Ax}. |
@param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step. |
| 502 |
@type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}. |
@type stoppingcriterium: function that returns C{True} or C{False} |
| 503 |
@param verbose: switches on the printing out some information |
@param x: an initial guess for the solution. If no C{x} is given 0*b is used. |
| 504 |
@type verbose: C{bool} |
@type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y) |
| 505 |
@param iter_max: maximum number of iteration steps. |
@param iter_max: maximum number of iteration steps. |
| 506 |
@type iter_max: C{int} |
@type iter_max: C{int} |
| 507 |
@param tolerance: tolerance |
@return: the solution approximation and the corresponding residual |
| 508 |
@type tolerance: positive C{float} |
@rtype: C{tuple} |
| 509 |
@return: the solution apprximation and the corresponding residual |
@warning: C{b} and C{x} are altered. |
|
@rtype: C{tuple} of an S type and and an R type object.A |
|
|
@warning: C{b} ans C{x} are altered. |
|
| 510 |
""" |
""" |
|
if verbose: |
|
|
print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance =%e"%(iter_max,tolerance) |
|
| 511 |
iter=0 |
iter=0 |
| 512 |
normb = norm(b) |
if x==None: |
| 513 |
if normb<0: raise NegativeNorm |
x=0*b |
| 514 |
|
else: |
| 515 |
b += (-1)*Aprod(x) |
b += (-1)*Aprod(x) |
| 516 |
r=b |
r=b |
| 517 |
rhat=Msolve(r) |
rhat=Msolve(r) |
| 518 |
d = rhat; |
d = rhat |
| 519 |
rhat_dot_r = bilinearform(rhat, r) |
rhat_dot_r = bilinearform(rhat, r) |
| 520 |
|
if rhat_dot_r<0: raise NegativeNorm,"negative norm." |
| 521 |
|
|
| 522 |
while True: |
while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x): |
| 523 |
normr=norm(r) |
iter+=1 |
|
if normr<0: raise NegativeNorm |
|
|
if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb) |
|
|
if normr <= tolerance * normb: return x,r |
|
|
|
|
|
iter+=1 # k = iter = 1 first time through |
|
| 524 |
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
| 525 |
|
|
| 526 |
q=Aprod(d) |
q=Aprod(d) |
| 535 |
d=rhat |
d=rhat |
| 536 |
|
|
| 537 |
rhat_dot_r = rhat_dot_r_new |
rhat_dot_r = rhat_dot_r_new |
| 538 |
|
if rhat_dot_r<0: raise NegativeNorm,"negative norm." |
| 539 |
|
|
| 540 |
|
return x,r |
| 541 |
|
|
| 542 |
class SaddlePointProblem(object): |
class SaddlePointProblem(object): |
| 543 |
""" |
""" |
Parent Directory
| 
Revision Log
| 
Patch