1 2 /* 3 4 Note that for the complex numbers version, the VecDot() arguments 5 within the code MUST remain in the order given for correct computation 6 of inner products. 7 */ 8 #include <private/kspimpl.h> 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "KSPSetUp_CGS" 12 static PetscErrorCode KSPSetUp_CGS(KSP ksp) 13 { 14 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 ierr = KSPDefaultGetWork(ksp,7);CHKERRQ(ierr); 18 PetscFunctionReturn(0); 19 } 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "KSPSolve_CGS" 23 static PetscErrorCode KSPSolve_CGS(KSP ksp) 24 { 25 PetscErrorCode ierr; 26 PetscInt i; 27 PetscScalar rho,rhoold,a,s,b; 28 Vec X,B,V,P,R,RP,T,Q,U,AUQ; 29 PetscReal dp = 0.0; 30 PetscBool diagonalscale; 31 32 PetscFunctionBegin; 33 /* not sure what residual norm it does use, should use for right preconditioning */ 34 35 ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr); 36 if (diagonalscale) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name); 37 38 X = ksp->vec_sol; 39 B = ksp->vec_rhs; 40 R = ksp->work[0]; 41 RP = ksp->work[1]; 42 V = ksp->work[2]; 43 T = ksp->work[3]; 44 Q = ksp->work[4]; 45 P = ksp->work[5]; 46 U = ksp->work[6]; 47 AUQ = V; 48 49 /* Compute initial preconditioned residual */ 50 ierr = KSPInitialResidual(ksp,X,V,T,R,B);CHKERRQ(ierr); 51 52 /* Test for nothing to do */ 53 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); 54 if (ksp->normtype == KSP_NORM_NATURAL) { 55 dp *= dp; 56 } 57 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 58 ksp->its = 0; 59 ksp->rnorm = dp; 60 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 61 KSPLogResidualHistory(ksp,dp); 62 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr); 63 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 64 if (ksp->reason) PetscFunctionReturn(0); 65 66 /* Make the initial Rp == R */ 67 ierr = VecCopy(R,RP);CHKERRQ(ierr); 68 /* added for Fidap */ 69 /* Penalize Startup - Isaac Hasbani Trick for CGS 70 Since most initial conditions result in a mostly 0 residual, 71 we change all the 0 values in the vector RP to the maximum. 72 */ 73 if (ksp->normtype == KSP_NORM_NATURAL) { 74 PetscReal vr0max; 75 PetscScalar *tmp_RP=0; 76 PetscInt numnp=0, *max_pos=0; 77 ierr = VecMax(RP, max_pos, &vr0max);CHKERRQ(ierr); 78 ierr = VecGetArray(RP, &tmp_RP);CHKERRQ(ierr); 79 ierr = VecGetLocalSize(RP, &numnp);CHKERRQ(ierr); 80 for (i=0; i<numnp; i++) { 81 if (tmp_RP[i] == 0.0) tmp_RP[i] = vr0max; 82 } 83 ierr = VecRestoreArray(RP, &tmp_RP);CHKERRQ(ierr); 84 } 85 /* end of addition for Fidap */ 86 87 /* Set the initial conditions */ 88 ierr = VecDot(R,RP,&rhoold);CHKERRQ(ierr); /* rhoold = (r,rp) */ 89 ierr = VecCopy(R,U);CHKERRQ(ierr); 90 ierr = VecCopy(R,P);CHKERRQ(ierr); 91 ierr = KSP_PCApplyBAorAB(ksp,P,V,T);CHKERRQ(ierr); 92 93 i = 0; 94 do { 95 96 ierr = VecDot(V,RP,&s);CHKERRQ(ierr); /* s <- (v,rp) */ 97 a = rhoold / s; /* a <- rho / s */ 98 ierr = VecWAXPY(Q,-a,V,U);CHKERRQ(ierr); /* q <- u - a v */ 99 ierr = VecWAXPY(T,1.0,U,Q);CHKERRQ(ierr); /* t <- u + q */ 100 ierr = VecAXPY(X,a,T);CHKERRQ(ierr); /* x <- x + a (u + q) */ 101 ierr = KSP_PCApplyBAorAB(ksp,T,AUQ,U);CHKERRQ(ierr); 102 ierr = VecAXPY(R,-a,AUQ);CHKERRQ(ierr); /* r <- r - a K (u + q) */ 103 ierr = VecDot(R,RP,&rho);CHKERRQ(ierr); /* rho <- (r,rp) */ 104 if (ksp->normtype == KSP_NORM_NATURAL) { 105 dp = PetscAbsScalar(rho); 106 } else { 107 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); 108 } 109 110 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 111 ksp->its++; 112 ksp->rnorm = dp; 113 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 114 KSPLogResidualHistory(ksp,dp); 115 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr); 116 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 117 if (ksp->reason) break; 118 119 b = rho / rhoold; /* b <- rho / rhoold */ 120 ierr = VecWAXPY(U,b,Q,R);CHKERRQ(ierr); /* u <- r + b q */ 121 ierr = VecAXPY(Q,b,P);CHKERRQ(ierr); 122 ierr = VecWAXPY(P,b,Q,U);CHKERRQ(ierr); /* p <- u + b(q + b p) */ 123 ierr = KSP_PCApplyBAorAB(ksp,P,V,Q);CHKERRQ(ierr); /* v <- K p */ 124 rhoold = rho; 125 i++; 126 } while (i<ksp->max_it); 127 if (i >= ksp->max_it) { 128 ksp->reason = KSP_DIVERGED_ITS; 129 } 130 131 ierr = KSPUnwindPreconditioner(ksp,X,T);CHKERRQ(ierr); 132 PetscFunctionReturn(0); 133 } 134 135 /*MC 136 KSPCGS - This code implements the CGS (Conjugate Gradient Squared) method. 137 138 Options Database Keys: 139 . see KSPSolve() 140 141 Level: beginner 142 143 References: Sonneveld, 1989. 144 145 Notes: Does not require a symmetric matrix. Does not apply transpose of the matrix. 146 Supports left and right preconditioning, but not symmetric. 147 148 Developer Notes: Has this weird support for doing the convergence test with the natural norm, I assume this works only with 149 no preconditioning and symmetric positive definite operator. 150 151 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBCGS, KSPSetPCSide() 152 M*/ 153 EXTERN_C_BEGIN 154 #undef __FUNCT__ 155 #define __FUNCT__ "KSPCreate_CGS" 156 PetscErrorCode KSPCreate_CGS(KSP ksp) 157 { 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ksp->data = (void*)0; 162 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 163 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);CHKERRQ(ierr); 164 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);CHKERRQ(ierr); 165 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_RIGHT,1);CHKERRQ(ierr); 166 167 ksp->ops->setup = KSPSetUp_CGS; 168 ksp->ops->solve = KSPSolve_CGS; 169 ksp->ops->destroy = KSPDefaultDestroy; 170 ksp->ops->buildsolution = KSPDefaultBuildSolution; 171 ksp->ops->buildresidual = KSPDefaultBuildResidual; 172 ksp->ops->setfromoptions = 0; 173 ksp->ops->view = 0; 174 PetscFunctionReturn(0); 175 } 176 EXTERN_C_END 177