689 |
stopped=False |
stopped=False |
690 |
|
|
691 |
return x,stopped |
return x,stopped |
692 |
|
|
693 |
|
def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100): |
694 |
|
|
695 |
|
# |
696 |
|
# minres solves the system of linear equations Ax = b |
697 |
|
# where A is a symmetric matrix (possibly indefinite or singular) |
698 |
|
# and b is a given vector. |
699 |
|
# |
700 |
|
# "A" may be a dense or sparse matrix (preferably sparse!) |
701 |
|
# or the name of a function such that |
702 |
|
# y = A(x) |
703 |
|
# returns the product y = Ax for any given vector x. |
704 |
|
# |
705 |
|
# "M" defines a positive-definite preconditioner M = C C'. |
706 |
|
# "M" may be a dense or sparse matrix (preferably sparse!) |
707 |
|
# or the name of a function such that |
708 |
|
# solves the system My = x for any given vector x. |
709 |
|
# |
710 |
|
# |
711 |
|
|
712 |
|
# Initialize |
713 |
|
|
714 |
|
iter = 0 |
715 |
|
Anorm = 0 |
716 |
|
ynorm = 0 |
717 |
|
x=x*0 |
718 |
|
#------------------------------------------------------------------ |
719 |
|
# Set up y and v for the first Lanczos vector v1. |
720 |
|
# y = beta1 P' v1, where P = C**(-1). |
721 |
|
# v is really P' v1. |
722 |
|
#------------------------------------------------------------------ |
723 |
|
r1 = b |
724 |
|
y = Msolve(b) |
725 |
|
beta1 = bilinearform(b,y) |
726 |
|
|
727 |
|
if beta1< 0: raise NegativeNorm,"negative norm." |
728 |
|
|
729 |
|
# If b = 0 exactly, stop with x = 0. |
730 |
|
if beta1==0: return x*0. |
731 |
|
|
732 |
|
if beta1> 0: |
733 |
|
beta1 = math.sqrt(beta1) # Normalize y to get v1 later. |
734 |
|
|
735 |
|
#------------------------------------------------------------------ |
736 |
|
# Initialize other quantities. |
737 |
|
# ------------------------------------------------------------------ |
738 |
|
oldb = 0 |
739 |
|
beta = beta1 |
740 |
|
dbar = 0 |
741 |
|
epsln = 0 |
742 |
|
phibar = beta1 |
743 |
|
rhs1 = beta1 |
744 |
|
rhs2 = 0 |
745 |
|
rnorm = phibar |
746 |
|
tnorm2 = 0 |
747 |
|
ynorm2 = 0 |
748 |
|
cs = -1 |
749 |
|
sn = 0 |
750 |
|
w = b*0. |
751 |
|
w2 = b*0. |
752 |
|
r2 = r1 |
753 |
|
eps = 0.0001 |
754 |
|
|
755 |
|
#--------------------------------------------------------------------- |
756 |
|
# Main iteration loop. |
757 |
|
# -------------------------------------------------------------------- |
758 |
|
while not stoppingcriterium(rnorm,Anorm*ynorm): # ||r|| / (||A|| ||x||) |
759 |
|
|
760 |
|
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
761 |
|
iter = iter + 1 |
762 |
|
|
763 |
|
#----------------------------------------------------------------- |
764 |
|
# Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,... |
765 |
|
# The general iteration is similar to the case k = 1 with v0 = 0: |
766 |
|
# |
767 |
|
# p1 = Operator * v1 - beta1 * v0, |
768 |
|
# alpha1 = v1'p1, |
769 |
|
# q2 = p2 - alpha1 * v1, |
770 |
|
# beta2^2 = q2'q2, |
771 |
|
# v2 = (1/beta2) q2. |
772 |
|
# |
773 |
|
# Again, y = betak P vk, where P = C**(-1). |
774 |
|
#----------------------------------------------------------------- |
775 |
|
s = 1/beta # Normalize previous vector (in y). |
776 |
|
v = s*y # v = vk if P = I |
777 |
|
|
778 |
|
y = Aprod(v) |
779 |
|
|
780 |
|
if iter >= 2: |
781 |
|
y = y - (beta/oldb)*r1 |
782 |
|
|
783 |
|
alfa = bilinearform(v,y) # alphak |
784 |
|
y = (- alfa/beta)*r2 + y |
785 |
|
r1 = r2 |
786 |
|
r2 = y |
787 |
|
y = Msolve(r2) |
788 |
|
oldb = beta # oldb = betak |
789 |
|
beta = bilinearform(r2,y) # beta = betak+1^2 |
790 |
|
if beta < 0: raise NegativeNorm,"negative norm." |
791 |
|
|
792 |
|
beta = math.sqrt( beta ) |
793 |
|
tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta |
794 |
|
|
795 |
|
if iter==1: # Initialize a few things. |
796 |
|
gmax = abs( alfa ) # alpha1 |
797 |
|
gmin = gmax # alpha1 |
798 |
|
|
799 |
|
# Apply previous rotation Qk-1 to get |
800 |
|
# [deltak epslnk+1] = [cs sn][dbark 0 ] |
801 |
|
# [gbar k dbar k+1] [sn -cs][alfak betak+1]. |
802 |
|
|
803 |
|
oldeps = epsln |
804 |
|
delta = cs * dbar + sn * alfa # delta1 = 0 deltak |
805 |
|
gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k |
806 |
|
epsln = sn * beta # epsln2 = 0 epslnk+1 |
807 |
|
dbar = - cs * beta # dbar 2 = beta2 dbar k+1 |
808 |
|
|
809 |
|
# Compute the next plane rotation Qk |
810 |
|
|
811 |
|
gamma = math.sqrt(gbar*gbar+beta*beta) # gammak |
812 |
|
gamma = max(gamma,eps) |
813 |
|
cs = gbar / gamma # ck |
814 |
|
sn = beta / gamma # sk |
815 |
|
phi = cs * phibar # phik |
816 |
|
phibar = sn * phibar # phibark+1 |
817 |
|
|
818 |
|
# Update x. |
819 |
|
|
820 |
|
denom = 1/gamma |
821 |
|
w1 = w2 |
822 |
|
w2 = w |
823 |
|
w = (v - oldeps*w1 - delta*w2) * denom |
824 |
|
x = x + phi*w |
825 |
|
|
826 |
|
# Go round again. |
827 |
|
|
828 |
|
gmax = max(gmax,gamma) |
829 |
|
gmin = min(gmin,gamma) |
830 |
|
z = rhs1 / gamma |
831 |
|
ynorm2 = z*z + ynorm2 |
832 |
|
rhs1 = rhs2 - delta*z |
833 |
|
rhs2 = - epsln*z |
834 |
|
|
835 |
|
# Estimate various norms and test for convergence. |
836 |
|
|
837 |
|
Anorm = math.sqrt( tnorm2 ) |
838 |
|
ynorm = math.sqrt( ynorm2 ) |
839 |
|
|
840 |
|
rnorm = phibar |
841 |
|
|
842 |
|
# Return final answer. |
843 |
|
return x |
844 |
|
|
845 |
############################################# |
############################################# |
846 |
|
|
1037 |
def __stoppingcriterium_GMRES(self,norm_r,norm_b): |
def __stoppingcriterium_GMRES(self,norm_r,norm_b): |
1038 |
return self.stoppingcriterium_GMRES(norm_r,norm_b) |
return self.stoppingcriterium_GMRES(norm_r,norm_b) |
1039 |
|
|
1040 |
|
def __stoppingcriterium_MINRES(self,norm_r,norm_Ax): |
1041 |
|
return self.stoppingcriterium_MINRES(norm_r,norm_Ax) |
1042 |
|
|
1043 |
|
|
1044 |
def setTolerance(self,tolerance=1.e-8): |
def setTolerance(self,tolerance=1.e-8): |
1045 |
self.__tol=tolerance |
self.__tol=tolerance |
1046 |
def getTolerance(self): |
def getTolerance(self): |
1084 |
# A(u-v)=f-B^*p-Av |
# A(u-v)=f-B^*p-Av |
1085 |
# u=v+(u-v) |
# u=v+(u-v) |
1086 |
u=v+self.solve_A(v,p) |
u=v+self.solve_A(v,p) |
1087 |
|
|
1088 |
|
if solver=='MINRES': |
1089 |
|
if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter |
1090 |
|
p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.) |
1091 |
|
# solve Au=f-B^*p |
1092 |
|
# A(u-v)=f-B^*p-Av |
1093 |
|
# u=v+(u-v) |
1094 |
|
u=v+self.solve_A(v,p) |
1095 |
|
|
1096 |
else: |
if solver=='PCG': |
1097 |
if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter |
if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter |
1098 |
p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p) |
p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p) |
1099 |
u=r[0] |
u=r[0] |
1100 |
|
|
1101 |
print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u) |
print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u) |
1102 |
|
|
1103 |
return u,p |
return u,p |