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