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