1 #include <../src/tao/bound/impls/bqnk/bqnk.h> /*I "petsctao.h" I*/ /*I "petscmat.h" I*/ 2 #include <petscksp.h> 3 4 static const char *BQNK_INIT[64] = {"constant", "direction"}; 5 static const char *BNK_UPDATE[64] = {"step", "reduction", "interpolation"}; 6 static const char *BNK_AS[64] = {"none", "bertsekas"}; 7 8 static PetscErrorCode TaoBQNKComputeHessian(Tao tao) 9 { 10 TAO_BNK *bnk = (TAO_BNK *)tao->data; 11 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 12 PetscErrorCode ierr; 13 PetscReal gnorm2, delta; 14 15 PetscFunctionBegin; 16 /* Alias the LMVM matrix into the TAO hessian */ 17 if (tao->hessian) { 18 ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr); 19 } 20 if (tao->hessian_pre) { 21 ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr); 22 } 23 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 24 tao->hessian = bqnk->B; 25 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 26 tao->hessian_pre = bqnk->B; 27 /* Update the Hessian with the latest solution */ 28 if (bqnk->is_spd) { 29 gnorm2 = bnk->gnorm*bnk->gnorm; 30 if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON; 31 if (bnk->f == 0.0) { 32 delta = 2.0 / gnorm2; 33 } else { 34 delta = 2.0 * PetscAbsScalar(bnk->f) / gnorm2; 35 } 36 ierr = MatLMVMSymBroydenSetDelta(bqnk->B, delta);CHKERRQ(ierr); 37 } 38 ierr = MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 39 ierr = MatLMVMResetShift(tao->hessian);CHKERRQ(ierr); 40 /* Prepare the reduced sub-matrices for the inactive set */ 41 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 42 if (bnk->active_idx) { 43 ierr = MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);CHKERRQ(ierr); 44 ierr = PCLMVMSetIS(bqnk->pc, bnk->inactive_idx);CHKERRQ(ierr); 45 } else { 46 ierr = PetscObjectReference((PetscObject)tao->hessian);CHKERRQ(ierr); 47 bnk->H_inactive = tao->hessian; 48 ierr = PCLMVMClearIS(bqnk->pc);CHKERRQ(ierr); 49 } 50 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 51 ierr = PetscObjectReference((PetscObject)bnk->H_inactive);CHKERRQ(ierr); 52 bnk->Hpre_inactive = bnk->H_inactive; 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type) 57 { 58 TAO_BNK *bnk = (TAO_BNK *)tao->data; 59 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 ierr = TaoBNKComputeStep(tao, shift, ksp_reason, step_type);CHKERRQ(ierr); 64 if (*ksp_reason < 0) { 65 /* Krylov solver failed to converge so reset the LMVM matrix */ 66 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 67 ierr = MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 68 } 69 PetscFunctionReturn(0); 70 } 71 72 PetscErrorCode TaoSolve_BQNK(Tao tao) 73 { 74 TAO_BNK *bnk = (TAO_BNK *)tao->data; 75 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 76 Mat_LMVM *lmvm = (Mat_LMVM*)bqnk->B->data; 77 Mat_LMVM *J0; 78 Mat_SymBrdn *diag_ctx; 79 PetscBool flg = PETSC_FALSE; 80 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 if (!tao->recycle) { 84 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 85 lmvm->nresets = 0; 86 if (lmvm->J0) { 87 ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &flg);CHKERRQ(ierr); 88 if (flg) { 89 J0 = (Mat_LMVM*)lmvm->J0->data; 90 J0->nresets = 0; 91 } 92 } 93 flg = PETSC_FALSE; 94 ierr = PetscObjectTypeCompareAny((PetscObject)bqnk->B, &flg, MATLMVMSYMBROYDEN, MATLMVMSYMBADBROYDEN, MATLMVMBFGS, MATLMVMDFP, "");CHKERRQ(ierr); 95 if (flg) { 96 diag_ctx = (Mat_SymBrdn*)lmvm->ctx; 97 J0 = (Mat_LMVM*)diag_ctx->D->data; 98 J0->nresets = 0; 99 } 100 } 101 ierr = (*bqnk->solve)(tao);CHKERRQ(ierr); 102 PetscFunctionReturn(0); 103 } 104 105 PetscErrorCode TaoSetUp_BQNK(Tao tao) 106 { 107 TAO_BNK *bnk = (TAO_BNK *)tao->data; 108 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 109 PetscErrorCode ierr; 110 PetscInt n, N; 111 PetscBool is_lmvm, is_sym, is_spd; 112 113 PetscFunctionBegin; 114 ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr); 115 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 116 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 117 ierr = MatSetSizes(bqnk->B, n, n, N, N);CHKERRQ(ierr); 118 ierr = MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);CHKERRQ(ierr); 119 ierr = PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);CHKERRQ(ierr); 120 if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type"); 121 ierr = MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);CHKERRQ(ierr); 122 if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric"); 123 ierr = MatGetOption(bqnk->B, MAT_SPD, &is_spd);CHKERRQ(ierr); 124 ierr = KSPGetPC(tao->ksp, &bqnk->pc);CHKERRQ(ierr); 125 ierr = PCSetType(bqnk->pc, PCLMVM);CHKERRQ(ierr); 126 ierr = PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 130 static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject,Tao tao) 131 { 132 TAO_BNK *bnk = (TAO_BNK *)tao->data; 133 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 134 PetscErrorCode ierr; 135 136 PetscFunctionBegin; 137 ierr = PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");CHKERRQ(ierr); 138 ierr = PetscOptionsEList("-tao_bqnk_init_type", "radius initialization type", "", BQNK_INIT, BQNK_INIT_TYPES, BQNK_INIT[bnk->init_type], &bnk->init_type, NULL);CHKERRQ(ierr); 139 ierr = PetscOptionsEList("-tao_bqnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, NULL);CHKERRQ(ierr); 140 ierr = PetscOptionsEList("-tao_bqnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, NULL);CHKERRQ(ierr); 141 ierr = PetscOptionsReal("-tao_bqnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval,NULL);CHKERRQ(ierr); 142 ierr = PetscOptionsReal("-tao_bqnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin,NULL);CHKERRQ(ierr); 143 ierr = PetscOptionsReal("-tao_bqnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax,NULL);CHKERRQ(ierr); 144 ierr = PetscOptionsReal("-tao_bqnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac,NULL);CHKERRQ(ierr); 145 ierr = PetscOptionsReal("-tao_bqnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin,NULL);CHKERRQ(ierr); 146 ierr = PetscOptionsReal("-tao_bqnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax,NULL);CHKERRQ(ierr); 147 ierr = PetscOptionsReal("-tao_bqnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac,NULL);CHKERRQ(ierr); 148 ierr = PetscOptionsReal("-tao_bqnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac,NULL);CHKERRQ(ierr); 149 ierr = PetscOptionsReal("-tao_bqnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac,NULL);CHKERRQ(ierr); 150 ierr = PetscOptionsReal("-tao_bqnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac,NULL);CHKERRQ(ierr); 151 ierr = PetscOptionsReal("-tao_bqnk_eta1", "(developer) threshold for rejecting step (-tao_bqnk_update_type reduction)", "", bnk->eta1, &bnk->eta1,NULL);CHKERRQ(ierr); 152 ierr = PetscOptionsReal("-tao_bqnk_eta2", "(developer) threshold for accepting marginal step (-tao_bqnk_update_type reduction)", "", bnk->eta2, &bnk->eta2,NULL);CHKERRQ(ierr); 153 ierr = PetscOptionsReal("-tao_bqnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bqnk_update_type reduction)", "", bnk->eta3, &bnk->eta3,NULL);CHKERRQ(ierr); 154 ierr = PetscOptionsReal("-tao_bqnk_eta4", "(developer) threshold for accepting good step (-tao_bqnk_update_type reduction)", "", bnk->eta4, &bnk->eta4,NULL);CHKERRQ(ierr); 155 ierr = PetscOptionsReal("-tao_bqnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bqnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1,NULL);CHKERRQ(ierr); 156 ierr = PetscOptionsReal("-tao_bqnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bqnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2,NULL);CHKERRQ(ierr); 157 ierr = PetscOptionsReal("-tao_bqnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3,NULL);CHKERRQ(ierr); 158 ierr = PetscOptionsReal("-tao_bqnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4,NULL);CHKERRQ(ierr); 159 ierr = PetscOptionsReal("-tao_bqnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5,NULL);CHKERRQ(ierr); 160 ierr = PetscOptionsReal("-tao_bqnk_nu1", "(developer) threshold for small line-search step length (-tao_bqnk_update_type step)", "", bnk->nu1, &bnk->nu1,NULL);CHKERRQ(ierr); 161 ierr = PetscOptionsReal("-tao_bqnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bqnk_update_type step)", "", bnk->nu2, &bnk->nu2,NULL);CHKERRQ(ierr); 162 ierr = PetscOptionsReal("-tao_bqnk_nu3", "(developer) threshold for large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu3, &bnk->nu3,NULL);CHKERRQ(ierr); 163 ierr = PetscOptionsReal("-tao_bqnk_nu4", "(developer) threshold for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu4, &bnk->nu4,NULL);CHKERRQ(ierr); 164 ierr = PetscOptionsReal("-tao_bqnk_omega1", "(developer) radius reduction factor for very small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega1, &bnk->omega1,NULL);CHKERRQ(ierr); 165 ierr = PetscOptionsReal("-tao_bqnk_omega2", "(developer) radius reduction factor for small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega2, &bnk->omega2,NULL);CHKERRQ(ierr); 166 ierr = PetscOptionsReal("-tao_bqnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bqnk_update_type step)", "", bnk->omega3, &bnk->omega3,NULL);CHKERRQ(ierr); 167 ierr = PetscOptionsReal("-tao_bqnk_omega4", "(developer) radius increase factor for large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega4, &bnk->omega4,NULL);CHKERRQ(ierr); 168 ierr = PetscOptionsReal("-tao_bqnk_omega5", "(developer) radius increase factor for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega5, &bnk->omega5,NULL);CHKERRQ(ierr); 169 ierr = PetscOptionsReal("-tao_bqnk_mu1", "(developer) threshold for accepting very good step (-tao_bqnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1,NULL);CHKERRQ(ierr); 170 ierr = PetscOptionsReal("-tao_bqnk_mu2", "(developer) threshold for accepting good step (-tao_bqnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2,NULL);CHKERRQ(ierr); 171 ierr = PetscOptionsReal("-tao_bqnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1,NULL);CHKERRQ(ierr); 172 ierr = PetscOptionsReal("-tao_bqnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2,NULL);CHKERRQ(ierr); 173 ierr = PetscOptionsReal("-tao_bqnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3,NULL);CHKERRQ(ierr); 174 ierr = PetscOptionsReal("-tao_bqnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4,NULL);CHKERRQ(ierr); 175 ierr = PetscOptionsReal("-tao_bqnk_theta", "(developer) trust region interpolation factor (-tao_bqnk_update_type interpolation)", "", bnk->theta, &bnk->theta,NULL);CHKERRQ(ierr); 176 ierr = PetscOptionsReal("-tao_bqnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL);CHKERRQ(ierr); 177 ierr = PetscOptionsReal("-tao_bqnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL);CHKERRQ(ierr); 178 ierr = PetscOptionsReal("-tao_bqnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);CHKERRQ(ierr); 179 ierr = PetscOptionsReal("-tao_bqnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);CHKERRQ(ierr); 180 ierr = PetscOptionsReal("-tao_bqnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);CHKERRQ(ierr); 181 ierr = PetscOptionsInt("-tao_bqnk_max_cg_its", "number of BNCG iterations to take for each Newton step", "", bnk->max_cg_its, &bnk->max_cg_its,NULL);CHKERRQ(ierr); 182 ierr = PetscOptionsTail();CHKERRQ(ierr); 183 ierr = TaoSetFromOptions(bnk->bncg);CHKERRQ(ierr); 184 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 185 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 186 if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION; 187 ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr); 188 ierr = MatGetOption(bqnk->B, MAT_SPD, &bqnk->is_spd);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer) 193 { 194 TAO_BNK *bnk = (TAO_BNK*)tao->data; 195 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 196 PetscErrorCode ierr; 197 PetscBool isascii; 198 199 PetscFunctionBegin; 200 ierr = TaoView_BNK(tao, viewer);CHKERRQ(ierr); 201 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 202 if (isascii) { 203 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 204 ierr = MatView(bqnk->B, viewer);CHKERRQ(ierr); 205 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 206 } 207 PetscFunctionReturn(0); 208 } 209 210 static PetscErrorCode TaoDestroy_BQNK(Tao tao) 211 { 212 TAO_BNK *bnk = (TAO_BNK*)tao->data; 213 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 214 PetscErrorCode ierr; 215 216 PetscFunctionBegin; 217 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 218 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 219 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 220 ierr = PetscFree(bnk->ctx);CHKERRQ(ierr); 221 ierr = TaoDestroy_BNK(tao);CHKERRQ(ierr); 222 PetscFunctionReturn(0); 223 } 224 225 PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao) 226 { 227 TAO_BNK *bnk; 228 TAO_BQNK *bqnk; 229 PetscErrorCode ierr; 230 231 PetscFunctionBegin; 232 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 233 ierr = KSPSetOptionsPrefix(tao->ksp, "tao_bqnk_");CHKERRQ(ierr); 234 tao->ops->solve = TaoSolve_BQNK; 235 tao->ops->setfromoptions = TaoSetFromOptions_BQNK; 236 tao->ops->destroy = TaoDestroy_BQNK; 237 tao->ops->view = TaoView_BQNK; 238 tao->ops->setup = TaoSetUp_BQNK; 239 240 bnk = (TAO_BNK *)tao->data; 241 bnk->computehessian = TaoBQNKComputeHessian; 242 bnk->computestep = TaoBQNKComputeStep; 243 bnk->init_type = BNK_INIT_DIRECTION; 244 245 ierr = PetscNewLog(tao,&bqnk);CHKERRQ(ierr); 246 bnk->ctx = (void*)bqnk; 247 bqnk->is_spd = PETSC_TRUE; 248 249 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);CHKERRQ(ierr); 250 ierr = PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);CHKERRQ(ierr); 251 ierr = MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");CHKERRQ(ierr); 252 ierr = MatSetType(bqnk->B, MATLMVMSR1);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 /*@ 257 TaoGetLMVMMatrix - Returns a pointer to the internal LMVM matrix. Valid 258 only for quasi-Newton family of methods. 259 260 Input Parameters: 261 . tao - Tao solver context 262 263 Output Parameters: 264 . B - LMVM matrix 265 266 Level: advanced 267 268 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoSetLMVMMatrix() 269 @*/ 270 PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B) 271 { 272 TAO_BNK *bnk = (TAO_BNK*)tao->data; 273 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 274 PetscErrorCode ierr; 275 PetscBool flg = PETSC_FALSE; 276 277 PetscFunctionBegin; 278 ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr); 279 if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms"); 280 *B = bqnk->B; 281 PetscFunctionReturn(0); 282 } 283 284 /*@ 285 TaoSetLMVMMatrix - Sets an external LMVM matrix into the Tao solver. Valid 286 only for quasi-Newton family of methods. 287 288 QN family of methods create their own LMVM matrices and users who wish to 289 manipulate this matrix should use TaoGetLMVMMatrix() instead. 290 291 Input Parameters: 292 + tao - Tao solver context 293 - B - LMVM matrix 294 295 Level: advanced 296 297 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoGetLMVMMatrix() 298 @*/ 299 PetscErrorCode TaoSetLMVMMatrix(Tao tao, Mat B) 300 { 301 TAO_BNK *bnk = (TAO_BNK*)tao->data; 302 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 303 PetscErrorCode ierr; 304 PetscBool flg = PETSC_FALSE; 305 306 PetscFunctionBegin; 307 ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr); 308 if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms"); 309 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &flg);CHKERRQ(ierr); 310 if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given matrix is not an LMVM matrix"); 311 if (bqnk->B) { 312 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 313 } 314 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 315 bqnk->B = B; 316 PetscFunctionReturn(0); 317 } 318