1 #include <../src/tao/bound/impls/bqnk/bqnk.h> 2 #include <petscksp.h> 3 4 /* Routine for computing the approximate quasi-Newton Hessian */ 5 6 static PetscErrorCode TaoBQNKComputeHessian(Tao tao) 7 { 8 TAO_BNK *bnk = (TAO_BNK *)tao->data; 9 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 /* Alias the LMVM matrix into the TAO hessian */ 14 if (tao->hessian) { 15 ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr); 16 } 17 if (tao->hessian_pre) { 18 ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr); 19 } 20 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 21 tao->hessian = bqnk->B; 22 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 23 tao->hessian_pre = bqnk->B; 24 /* Update the Hessian with the latest solution */ 25 ierr = MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 26 ierr = MatLMVMResetShift(tao->hessian);CHKERRQ(ierr); 27 /* Prepare the reduced sub-matrices for the inactive set */ 28 if (bnk->active_idx) { 29 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 30 ierr = MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);CHKERRQ(ierr); 31 } else { 32 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 33 ierr = PetscObjectReference((PetscObject)tao->hessian);CHKERRQ(ierr); 34 bnk->H_inactive = tao->hessian; 35 } 36 bnk->Hpre_inactive = bnk->H_inactive; 37 PetscFunctionReturn(0); 38 } 39 40 static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason) 41 { 42 TAO_BNK *bnk = (TAO_BNK *)tao->data; 43 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = TaoBNKComputeStep(tao, shift, ksp_reason);CHKERRQ(ierr); 48 if (*ksp_reason < 0) { 49 /* Krylov solver failed to converge so reset the LMVM matrix */ 50 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 51 } 52 PetscFunctionReturn(0); 53 } 54 55 static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject, Tao tao) 56 { 57 TAO_BNK *bnk = (TAO_BNK*)tao->data; 58 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 59 PetscErrorCode ierr; 60 PC pc; 61 PetscBool is_lmvm, is_sym; 62 63 PetscFunctionBegin; 64 ierr = TaoSetFromOptions_BNK(PetscOptionsObject, tao);CHKERRQ(ierr); 65 /* Make sure the interpolation method is disabled for trust region initialization */ 66 if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION; 67 /* Make sure the KSP solver is a CG method and the preconditioner is disabled */ 68 if (!bnk->is_nash && !bnk->is_stcg && !bnk->is_gltr) SETERRQ(PETSC_COMM_SELF,1,"Must use a trust-region CG method for KSP (KSPNASH, KSPSTCG, KSPGLTR)"); 69 ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr); 70 ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr); 71 /* Read mat options and check type and properties */ 72 ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr); 73 ierr = PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);CHKERRQ(ierr); 74 if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type"); 75 ierr = MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);CHKERRQ(ierr); 76 if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric"); 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer) 81 { 82 TAO_BNK *bnk = (TAO_BNK*)tao->data; 83 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 84 PetscErrorCode ierr; 85 PetscBool isascii; 86 87 PetscFunctionBegin; 88 ierr = TaoView_BNK(tao, viewer);CHKERRQ(ierr); 89 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 90 if (isascii) { 91 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 92 ierr = MatView(bqnk->B, viewer);CHKERRQ(ierr); 93 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 94 } 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode TaoDestroy_BQNK(Tao tao) 99 { 100 TAO_BNK *bnk = (TAO_BNK*)tao->data; 101 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 102 PetscErrorCode ierr; 103 104 PetscFunctionBegin; 105 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 106 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 107 ierr = PetscFree(bnk->ctx);CHKERRQ(ierr); 108 ierr = TaoDestroy_BNK(tao);CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } 111 112 PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao) 113 { 114 TAO_BNK *bnk; 115 TAO_BQNK *bqnk; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 120 tao->ops->setfromoptions = TaoSetFromOptions_BQNK; 121 tao->ops->destroy = TaoDestroy_BQNK; 122 tao->ops->view = TaoView_BQNK; 123 124 bnk = (TAO_BNK *)tao->data; 125 bnk->computehessian = TaoBQNKComputeHessian; 126 bnk->computestep = TaoBQNKComputeStep; 127 bnk->init_type = BNK_INIT_DIRECTION; 128 129 ierr = PetscNewLog(tao,&bqnk);CHKERRQ(ierr); 130 bnk->ctx = (void*)bqnk; 131 132 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);CHKERRQ(ierr); 133 ierr = PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);CHKERRQ(ierr); 134 ierr = MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");CHKERRQ(ierr); 135 ierr = MatSetType(bqnk->B, MATLMVMSR1);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 }