1 
2 #include <petsc/private/kspimpl.h>
3 
4 typedef struct {
5   PetscReal haptol;
6 } KSP_MINRES;
7 
KSPSetUp_MINRES(KSP ksp)8 static PetscErrorCode KSPSetUp_MINRES(KSP ksp)
9 {
10   PetscErrorCode ierr;
11 
12   PetscFunctionBegin;
13   if (ksp->pc_side == PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"No right preconditioning for KSPMINRES");
14   else if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"No symmetric preconditioning for KSPMINRES");
15   ierr = KSPSetWorkVecs(ksp,9);CHKERRQ(ierr);
16   PetscFunctionReturn(0);
17 }
18 
19 
KSPSolve_MINRES(KSP ksp)20 static PetscErrorCode  KSPSolve_MINRES(KSP ksp)
21 {
22   PetscErrorCode    ierr;
23   PetscInt          i;
24   PetscScalar       alpha,beta,ibeta,betaold,eta,c=1.0,ceta,cold=1.0,coold,s=0.0,sold=0.0,soold;
25   PetscScalar       rho0,rho1,irho1,rho2,rho3,dp = 0.0;
26   const PetscScalar none = -1.0;
27   PetscReal         np;
28   Vec               X,B,R,Z,U,V,W,UOLD,VOLD,WOLD,WOOLD;
29   Mat               Amat,Pmat;
30   KSP_MINRES        *minres = (KSP_MINRES*)ksp->data;
31   PetscBool         diagonalscale;
32 
33   PetscFunctionBegin;
34   ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
35   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
36 
37   X     = ksp->vec_sol;
38   B     = ksp->vec_rhs;
39   R     = ksp->work[0];
40   Z     = ksp->work[1];
41   U     = ksp->work[2];
42   V     = ksp->work[3];
43   W     = ksp->work[4];
44   UOLD  = ksp->work[5];
45   VOLD  = ksp->work[6];
46   WOLD  = ksp->work[7];
47   WOOLD = ksp->work[8];
48 
49   ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
50 
51   ksp->its = 0;
52 
53   ierr = VecSet(UOLD,0.0);CHKERRQ(ierr);          /*     u_old  <-   0   */
54   ierr = VecSet(VOLD,0.0);CHKERRQ(ierr);         /*     v_old  <-   0   */
55   ierr = VecSet(W,0.0);CHKERRQ(ierr);            /*     w      <-   0   */
56   ierr = VecSet(WOLD,0.0);CHKERRQ(ierr);         /*     w_old  <-   0   */
57 
58   if (!ksp->guess_zero) {
59     ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /*     r <- b - A*x    */
60     ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);
61   } else {
62     ierr = VecCopy(B,R);CHKERRQ(ierr);              /*     r <- b (x is 0) */
63   }
64   ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);       /*     z  <- B*r       */
65   ierr = VecNorm(Z,NORM_2,&np);CHKERRQ(ierr);      /*   np <- ||z||        */
66   KSPCheckNorm(ksp,np);
67   ierr = VecDot(R,Z,&dp);CHKERRQ(ierr);
68   KSPCheckDot(ksp,dp);
69 
70   if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
71     if (ksp->errorifnotconverged) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Detected indefinite operator %g tolerance %g",(double)PetscRealPart(dp),(double)minres->haptol);
72     ierr = PetscInfo2(ksp,"Detected indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);CHKERRQ(ierr);
73     ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
74     PetscFunctionReturn(0);
75   }
76 
77   ksp->rnorm = 0.0;
78   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
79   ierr = KSPLogResidualHistory(ksp,ksp->rnorm);CHKERRQ(ierr);
80   ierr = KSPMonitor(ksp,0,ksp->rnorm);CHKERRQ(ierr);
81   ierr = (*ksp->converged)(ksp,0,ksp->rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */
82   if (ksp->reason) PetscFunctionReturn(0);
83 
84   dp   = PetscAbsScalar(dp);
85   dp   = PetscSqrtScalar(dp);
86   beta = dp;                                        /*  beta <- sqrt(r'*z  */
87   eta  = beta;
88 
89   ierr  = VecCopy(R,V);CHKERRQ(ierr);
90   ierr  = VecCopy(Z,U);CHKERRQ(ierr);
91   ibeta = 1.0 / beta;
92   ierr  = VecScale(V,ibeta);CHKERRQ(ierr);        /*    v <- r / beta     */
93   ierr  = VecScale(U,ibeta);CHKERRQ(ierr);        /*    u <- z / beta     */
94 
95   i = 0;
96   do {
97     ksp->its = i+1;
98 
99     /*   Lanczos  */
100 
101     ierr = KSP_MatMult(ksp,Amat,U,R);CHKERRQ(ierr);   /*      r <- A*u   */
102     ierr = VecDot(U,R,&alpha);CHKERRQ(ierr);          /*  alpha <- r'*u  */
103     ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /*      z <- B*r   */
104 
105     ierr = VecAXPY(R,-alpha,V);CHKERRQ(ierr);     /*  r <- r - alpha v     */
106     ierr = VecAXPY(Z,-alpha,U);CHKERRQ(ierr);     /*  z <- z - alpha u     */
107     ierr = VecAXPY(R,-beta,VOLD);CHKERRQ(ierr);   /*  r <- r - beta v_old  */
108     ierr = VecAXPY(Z,-beta,UOLD);CHKERRQ(ierr);   /*  z <- z - beta u_old  */
109 
110     betaold = beta;
111 
112     ierr = VecDot(R,Z,&dp);CHKERRQ(ierr);
113     KSPCheckDot(ksp,dp);
114     dp   = PetscAbsScalar(dp);
115     beta = PetscSqrtScalar(dp);                               /*  beta <- sqrt(r'*z)   */
116 
117     /*    QR factorisation    */
118 
119     coold = cold; cold = c; soold = sold; sold = s;
120 
121     rho0 = cold * alpha - coold * sold * betaold;
122     rho1 = PetscSqrtScalar(rho0*rho0 + beta*beta);
123     rho2 = sold * alpha + coold * cold * betaold;
124     rho3 = soold * betaold;
125 
126     /*     Givens rotation    */
127 
128     c = rho0 / rho1;
129     s = beta / rho1;
130 
131     /*    Update    */
132 
133     ierr = VecCopy(WOLD,WOOLD);CHKERRQ(ierr);     /*  w_oold <- w_old      */
134     ierr = VecCopy(W,WOLD);CHKERRQ(ierr);         /*  w_old  <- w          */
135 
136     ierr  = VecCopy(U,W);CHKERRQ(ierr);           /*  w      <- u          */
137     ierr  = VecAXPY(W,-rho2,WOLD);CHKERRQ(ierr); /*  w <- w - rho2 w_old  */
138     ierr  = VecAXPY(W,-rho3,WOOLD);CHKERRQ(ierr); /*  w <- w - rho3 w_oold */
139     irho1 = 1.0 / rho1;
140     ierr  = VecScale(W,irho1);CHKERRQ(ierr);     /*  w <- w / rho1        */
141 
142     ceta = c * eta;
143     ierr = VecAXPY(X,ceta,W);CHKERRQ(ierr);      /*  x <- x + c eta w     */
144 
145     /*
146         when dp is really small we have either convergence or an indefinite operator so compute true
147         residual norm to check for convergence
148     */
149     if (PetscRealPart(dp) < minres->haptol) {
150       ierr = PetscInfo2(ksp,"Possible indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);CHKERRQ(ierr);
151       ierr = KSP_MatMult(ksp,Amat,X,VOLD);CHKERRQ(ierr);
152       ierr = VecAXPY(VOLD,none,B);CHKERRQ(ierr);
153       ierr = VecNorm(VOLD,NORM_2,&np);CHKERRQ(ierr);
154       KSPCheckNorm(ksp,np);
155     } else {
156       /* otherwise compute new residual norm via recurrence relation */
157       np *= PetscAbsScalar(s);
158     }
159 
160     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
161     ierr = KSPLogResidualHistory(ksp,ksp->rnorm);CHKERRQ(ierr);
162     ierr = KSPMonitor(ksp,i+1,ksp->rnorm);CHKERRQ(ierr);
163     ierr = (*ksp->converged)(ksp,i+1,ksp->rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */
164     if (ksp->reason) break;
165 
166     if (PetscRealPart(dp) < minres->haptol) {
167       if (ksp->errorifnotconverged) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_CONV_FAILED,"Detected indefinite operator %g tolerance %g",(double)PetscRealPart(dp),(double)minres->haptol);
168       ierr = PetscInfo2(ksp,"Detected indefinite operator %g tolerance %g\n",(double)PetscRealPart(dp),(double)minres->haptol);CHKERRQ(ierr);
169       ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
170       break;
171     }
172 
173     eta  = -s * eta;
174     ierr  = VecCopy(V,VOLD);CHKERRQ(ierr);
175     ierr  = VecCopy(U,UOLD);CHKERRQ(ierr);
176     ierr  = VecCopy(R,V);CHKERRQ(ierr);
177     ierr  = VecCopy(Z,U);CHKERRQ(ierr);
178     ibeta = 1.0 / beta;
179     ierr  = VecScale(V,ibeta);CHKERRQ(ierr);     /*  v <- r / beta       */
180     ierr  = VecScale(U,ibeta);CHKERRQ(ierr);     /*  u <- z / beta       */
181 
182     i++;
183   } while (i<ksp->max_it);
184   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
185   PetscFunctionReturn(0);
186 }
187 
188 /*MC
189      KSPMINRES - This code implements the MINRES (Minimum Residual) method.
190 
191    Options Database Keys:
192 .   see KSPSolve()
193 
194    Level: beginner
195 
196    Notes:
197     The operator and the preconditioner must be symmetric and the preconditioner must
198           be positive definite for this method.
199           Supports only left preconditioning.
200 
201    Reference: Paige & Saunders, 1975.
202 
203    Contributed by: Robert Scheichl: maprs@maths.bath.ac.uk
204 
205 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG, KSPCR
206 M*/
KSPCreate_MINRES(KSP ksp)207 PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
208 {
209   KSP_MINRES     *minres;
210   PetscErrorCode ierr;
211 
212   PetscFunctionBegin;
213   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr);
214   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
215   ierr = PetscNewLog(ksp,&minres);CHKERRQ(ierr);
216 
217   /* this parameter is arbitrary; but e-50 didn't work for __float128 in one example */
218 #if defined(PETSC_USE_REAL___FLOAT128)
219   minres->haptol = 1.e-100;
220 #elif defined(PETSC_USE_REAL_SINGLE)
221   minres->haptol = 1.e-25;
222 #else
223   minres->haptol = 1.e-50;
224 #endif
225   ksp->data      = (void*)minres;
226 
227   /*
228        Sets the functions that are associated with this data structure
229        (in C++ this is the same as defining virtual functions)
230   */
231   ksp->ops->setup          = KSPSetUp_MINRES;
232   ksp->ops->solve          = KSPSolve_MINRES;
233   ksp->ops->destroy        = KSPDestroyDefault;
234   ksp->ops->setfromoptions = NULL;
235   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
236   ksp->ops->buildresidual  = KSPBuildResidualDefault;
237   PetscFunctionReturn(0);
238 }
239