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