1 2 #include <petsc/private/kspimpl.h> 3 4 typedef struct { 5 PetscReal haptol; 6 } KSP_SYMMLQ; 7 8 PetscErrorCode KSPSetUp_SYMMLQ(KSP ksp) 9 { 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = KSPSetWorkVecs(ksp,9);CHKERRQ(ierr); 14 PetscFunctionReturn(0); 15 } 16 17 PetscErrorCode KSPSolve_SYMMLQ(KSP ksp) 18 { 19 PetscErrorCode ierr; 20 PetscInt i; 21 PetscScalar alpha,beta,ibeta,betaold,beta1,ceta = 0,ceta_oold = 0.0, ceta_old = 0.0,ceta_bar; 22 PetscScalar c = 1.0,cold=1.0,s=0.0,sold=0.0,coold,soold,rho0,rho1,rho2,rho3; 23 PetscScalar dp = 0.0; 24 PetscReal np,s_prod; 25 Vec X,B,R,Z,U,V,W,UOLD,VOLD,Wbar; 26 Mat Amat,Pmat; 27 KSP_SYMMLQ *symmlq = (KSP_SYMMLQ*)ksp->data; 28 PetscBool diagonalscale; 29 30 PetscFunctionBegin; 31 ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr); 32 if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name); 33 34 X = ksp->vec_sol; 35 B = ksp->vec_rhs; 36 R = ksp->work[0]; 37 Z = ksp->work[1]; 38 U = ksp->work[2]; 39 V = ksp->work[3]; 40 W = ksp->work[4]; 41 UOLD = ksp->work[5]; 42 VOLD = ksp->work[6]; 43 Wbar = ksp->work[7]; 44 45 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr); 46 47 ksp->its = 0; 48 49 ierr = VecSet(UOLD,0.0);CHKERRQ(ierr); /* u_old <- zeros; */ 50 ierr = VecCopy(UOLD,VOLD);CHKERRQ(ierr); /* v_old <- u_old; */ 51 ierr = VecCopy(UOLD,W);CHKERRQ(ierr); /* w <- u_old; */ 52 ierr = VecCopy(UOLD,Wbar);CHKERRQ(ierr); /* w_bar <- u_old; */ 53 if (!ksp->guess_zero) { 54 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /* r <- b - A*x */ 55 ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr); 56 } else { 57 ierr = VecCopy(B,R);CHKERRQ(ierr); /* r <- b (x is 0) */ 58 } 59 60 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- B*r */ 61 ierr = VecDot(R,Z,&dp);CHKERRQ(ierr); /* dp = r'*z; */ 62 if (PetscAbsScalar(dp) < symmlq->haptol) { 63 ierr = PetscInfo2(ksp,"Detected happy breakdown %g tolerance %g\n",(double)PetscAbsScalar(dp),(double)symmlq->haptol);CHKERRQ(ierr); 64 ksp->rnorm = 0.0; /* what should we really put here? */ 65 ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN; /* bugfix proposed by Lourens (lourens.vanzanen@shell.com) */ 66 PetscFunctionReturn(0); 67 } 68 69 #if !defined(PETSC_USE_COMPLEX) 70 if (dp < 0.0) { 71 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 72 PetscFunctionReturn(0); 73 } 74 #endif 75 dp = PetscSqrtScalar(dp); 76 beta = dp; /* beta <- sqrt(r'*z) */ 77 beta1 = beta; 78 s_prod = PetscAbsScalar(beta1); 79 80 ierr = VecCopy(R,V);CHKERRQ(ierr); /* v <- r; */ 81 ierr = VecCopy(Z,U);CHKERRQ(ierr); /* u <- z; */ 82 ibeta = 1.0 / beta; 83 ierr = VecScale(V,ibeta);CHKERRQ(ierr); /* v <- ibeta*v; */ 84 ierr = VecScale(U,ibeta);CHKERRQ(ierr); /* u <- ibeta*u; */ 85 ierr = VecCopy(U,Wbar);CHKERRQ(ierr); /* w_bar <- u; */ 86 ierr = VecNorm(Z,NORM_2,&np);CHKERRQ(ierr); /* np <- ||z|| */ 87 ierr = KSPLogResidualHistory(ksp,np);CHKERRQ(ierr); 88 ierr = KSPMonitor(ksp,0,np);CHKERRQ(ierr); 89 ksp->rnorm = np; 90 ierr = (*ksp->converged)(ksp,0,np,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */ 91 if (ksp->reason) PetscFunctionReturn(0); 92 93 i = 0; ceta = 0.; 94 do { 95 ksp->its = i+1; 96 97 /* Update */ 98 if (ksp->its > 1) { 99 ierr = VecCopy(V,VOLD);CHKERRQ(ierr); /* v_old <- v; */ 100 ierr = VecCopy(U,UOLD);CHKERRQ(ierr); /* u_old <- u; */ 101 102 ierr = VecCopy(R,V);CHKERRQ(ierr); 103 ierr = VecScale(V,1.0/beta);CHKERRQ(ierr); /* v <- ibeta*r; */ 104 ierr = VecCopy(Z,U);CHKERRQ(ierr); 105 ierr = VecScale(U,1.0/beta);CHKERRQ(ierr); /* u <- ibeta*z; */ 106 107 ierr = VecCopy(Wbar,W);CHKERRQ(ierr); 108 ierr = VecScale(W,c);CHKERRQ(ierr); 109 ierr = VecAXPY(W,s,U);CHKERRQ(ierr); /* w <- c*w_bar + s*u; (w_k) */ 110 ierr = VecScale(Wbar,-s);CHKERRQ(ierr); 111 ierr = VecAXPY(Wbar,c,U);CHKERRQ(ierr); /* w_bar <- -s*w_bar + c*u; (w_bar_(k+1)) */ 112 ierr = VecAXPY(X,ceta,W);CHKERRQ(ierr); /* x <- x + ceta * w; (xL_k) */ 113 114 ceta_oold = ceta_old; 115 ceta_old = ceta; 116 } 117 118 /* Lanczos */ 119 ierr = KSP_MatMult(ksp,Amat,U,R);CHKERRQ(ierr); /* r <- Amat*u; */ 120 ierr = VecDot(U,R,&alpha);CHKERRQ(ierr); /* alpha <- u'*r; */ 121 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- B*r; */ 122 123 ierr = VecAXPY(R,-alpha,V);CHKERRQ(ierr); /* r <- r - alpha* v; */ 124 ierr = VecAXPY(Z,-alpha,U);CHKERRQ(ierr); /* z <- z - alpha* u; */ 125 ierr = VecAXPY(R,-beta,VOLD);CHKERRQ(ierr); /* r <- r - beta * v_old; */ 126 ierr = VecAXPY(Z,-beta,UOLD);CHKERRQ(ierr); /* z <- z - beta * u_old; */ 127 betaold = beta; /* beta_k */ 128 ierr = VecDot(R,Z,&dp);CHKERRQ(ierr); /* dp <- r'*z; */ 129 if (PetscAbsScalar(dp) < symmlq->haptol) { 130 ierr = PetscInfo2(ksp,"Detected happy breakdown %g tolerance %g\n",(double)PetscAbsScalar(dp),(double)symmlq->haptol);CHKERRQ(ierr); 131 dp = 0.0; 132 } 133 134 #if !defined(PETSC_USE_COMPLEX) 135 if (dp < 0.0) { 136 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 137 break; 138 } 139 #endif 140 beta = PetscSqrtScalar(dp); /* beta = sqrt(dp); */ 141 142 /* QR factorization */ 143 coold = cold; cold = c; soold = sold; sold = s; 144 rho0 = cold * alpha - coold * sold * betaold; /* gamma_bar */ 145 rho1 = PetscSqrtScalar(rho0*rho0 + beta*beta); /* gamma */ 146 rho2 = sold * alpha + coold * cold * betaold; /* delta */ 147 rho3 = soold * betaold; /* epsilon */ 148 149 /* Givens rotation: [c -s; s c] (different from the Reference!) */ 150 c = rho0 / rho1; s = beta / rho1; 151 152 if (ksp->its==1) ceta = beta1/rho1; 153 else ceta = -(rho2*ceta_old + rho3*ceta_oold)/rho1; 154 155 s_prod = s_prod*PetscAbsScalar(s); 156 if (c == 0.0) np = s_prod*1.e16; 157 else np = s_prod/PetscAbsScalar(c); /* residual norm for xc_k (CGNORM) */ 158 159 ksp->rnorm = np; 160 ierr = KSPLogResidualHistory(ksp,np);CHKERRQ(ierr); 161 ierr = KSPMonitor(ksp,i+1,np);CHKERRQ(ierr); 162 ierr = (*ksp->converged)(ksp,i+1,np,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */ 163 if (ksp->reason) break; 164 i++; 165 } while (i<ksp->max_it); 166 167 /* move to the CG point: xc_(k+1) */ 168 if (c == 0.0) ceta_bar = ceta*1.e15; 169 else ceta_bar = ceta/c; 170 171 ierr = VecAXPY(X,ceta_bar,Wbar);CHKERRQ(ierr); /* x <- x + ceta_bar*w_bar */ 172 173 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 174 PetscFunctionReturn(0); 175 } 176 177 /*MC 178 KSPSYMMLQ - This code implements the SYMMLQ method. 179 180 Options Database Keys: 181 . see KSPSolve() 182 183 Level: beginner 184 185 Notes: The operator and the preconditioner must be symmetric for this method. The 186 preconditioner must be POSITIVE-DEFINITE. 187 188 Supports only left preconditioning. 189 190 Reference: Paige & Saunders, 1975. 191 192 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP 193 M*/ 194 PETSC_EXTERN PetscErrorCode KSPCreate_SYMMLQ(KSP ksp) 195 { 196 KSP_SYMMLQ *symmlq; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr); 201 202 ierr = PetscNewLog(ksp,&symmlq);CHKERRQ(ierr); 203 symmlq->haptol = 1.e-18; 204 ksp->data = (void*)symmlq; 205 206 /* 207 Sets the functions that are associated with this data structure 208 (in C++ this is the same as defining virtual functions) 209 */ 210 ksp->ops->setup = KSPSetUp_SYMMLQ; 211 ksp->ops->solve = KSPSolve_SYMMLQ; 212 ksp->ops->destroy = KSPDestroyDefault; 213 ksp->ops->setfromoptions = 0; 214 ksp->ops->buildsolution = KSPBuildSolutionDefault; 215 ksp->ops->buildresidual = KSPBuildResidualDefault; 216 PetscFunctionReturn(0); 217 } 218 219 220 221 222 223