1 #include <../src/tao/bound/impls/bqnk/bqnk.h> 2 #include <petscksp.h> 3 4 static PetscErrorCode TaoBQNKComputeHessian(Tao tao) 5 { 6 TAO_BNK *bnk = (TAO_BNK *)tao->data; 7 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 /* Alias the LMVM matrix into the TAO hessian */ 12 if (tao->hessian) { 13 ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr); 14 } 15 if (tao->hessian_pre) { 16 ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr); 17 } 18 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 19 tao->hessian = bqnk->B; 20 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 21 tao->hessian_pre = bqnk->B; 22 /* Update the Hessian with the latest solution */ 23 ierr = MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 24 ierr = MatLMVMResetShift(tao->hessian);CHKERRQ(ierr); 25 /* Prepare the reduced sub-matrices for the inactive set */ 26 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 27 if (bnk->active_idx) { 28 ierr = MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);CHKERRQ(ierr); 29 if (bqnk->pc) { 30 ierr = PCLMVMSetIS(bqnk->pc, bnk->inactive_idx);CHKERRQ(ierr); 31 } 32 } else { 33 ierr = PetscObjectReference((PetscObject)tao->hessian);CHKERRQ(ierr); 34 bnk->H_inactive = tao->hessian; 35 if (bqnk->pc) { 36 ierr = PCLMVMClearIS(bqnk->pc);CHKERRQ(ierr); 37 } 38 } 39 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 40 ierr = PetscObjectReference((PetscObject)bnk->H_inactive);CHKERRQ(ierr); 41 bnk->Hpre_inactive = bnk->H_inactive; 42 PetscFunctionReturn(0); 43 } 44 45 static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason) 46 { 47 TAO_BNK *bnk = (TAO_BNK *)tao->data; 48 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 ierr = TaoBNKComputeStep(tao, shift, ksp_reason);CHKERRQ(ierr); 53 if (*ksp_reason < 0) { 54 /* Krylov solver failed to converge so reset the LMVM matrix */ 55 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 56 } 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode TaoSetUp_BQNK(Tao tao) 61 { 62 TAO_BNK *bnk = (TAO_BNK *)tao->data; 63 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 64 PetscErrorCode ierr; 65 PetscInt n, N; 66 PetscBool is_lmvm, is_sym, is_sr1; 67 68 PetscFunctionBegin; 69 ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr); 70 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 71 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 72 ierr = MatSetSizes(bqnk->B, n, n, N, N);CHKERRQ(ierr); 73 ierr = MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);CHKERRQ(ierr); 74 ierr = PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);CHKERRQ(ierr); 75 if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type"); 76 ierr = MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);CHKERRQ(ierr); 77 if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric"); 78 ierr = PetscObjectTypeCompare((PetscObject)bqnk->B, MATLMVMSR1, &is_sr1);CHKERRQ(ierr); 79 ierr = KSPGetPC(tao->ksp, &bqnk->pc);CHKERRQ(ierr); 80 if (is_sr1) { 81 ierr = PCSetType(bqnk->pc, PCNONE);CHKERRQ(ierr); 82 bqnk->pc = NULL; 83 bqnk->no_scale = PETSC_TRUE; 84 } else { 85 ierr = PCSetType(bqnk->pc, PCLMVM);CHKERRQ(ierr); 86 ierr = PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);CHKERRQ(ierr); 87 } 88 if (!bqnk->no_scale) { 89 if (!bqnk->Bscale) { 90 ierr = MatCreateLMVMDiagBrdn(PetscObjectComm((PetscObject)bqnk->B), n, N, &bqnk->Bscale);CHKERRQ(ierr); 91 ierr = MatSetOptionsPrefix(bqnk->Bscale, "tao_bqnk_scale_");CHKERRQ(ierr); 92 ierr = MatLMVMAllocate(bqnk->Bscale, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 93 } 94 ierr = MatLMVMSetJ0(bqnk->B, bqnk->Bscale);CHKERRQ(ierr); 95 } 96 PetscFunctionReturn(0); 97 } 98 99 static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject, Tao tao) 100 { 101 TAO_BNK *bnk = (TAO_BNK*)tao->data; 102 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 103 PetscErrorCode ierr; 104 105 PetscFunctionBegin; 106 ierr = PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");CHKERRQ(ierr); 107 ierr = PetscOptionsBool("-tao_bqnk_no_scale","(developer) disable the diagonal Broyden scaling of the BFGS approximation","",bqnk->no_scale,&bqnk->no_scale,NULL);CHKERRQ(ierr); 108 ierr = PetscOptionsTail();CHKERRQ(ierr); 109 ierr = TaoSetFromOptions_BNK(PetscOptionsObject, tao);CHKERRQ(ierr); 110 if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION; 111 ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr); 112 PetscFunctionReturn(0); 113 } 114 115 static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer) 116 { 117 TAO_BNK *bnk = (TAO_BNK*)tao->data; 118 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 119 PetscErrorCode ierr; 120 PetscBool isascii; 121 122 PetscFunctionBegin; 123 ierr = TaoView_BNK(tao, viewer);CHKERRQ(ierr); 124 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 125 if (isascii) { 126 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 127 ierr = MatView(bqnk->B, viewer);CHKERRQ(ierr); 128 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 129 } 130 PetscFunctionReturn(0); 131 } 132 133 static PetscErrorCode TaoDestroy_BQNK(Tao tao) 134 { 135 TAO_BNK *bnk = (TAO_BNK*)tao->data; 136 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 137 PetscErrorCode ierr; 138 139 PetscFunctionBegin; 140 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 141 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 142 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 143 if (!bqnk->no_scale) { 144 ierr = MatDestroy(&bqnk->Bscale);CHKERRQ(ierr); 145 } 146 ierr = PetscFree(bnk->ctx);CHKERRQ(ierr); 147 ierr = TaoDestroy_BNK(tao);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao) 152 { 153 TAO_BNK *bnk; 154 TAO_BQNK *bqnk; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 159 tao->ops->setfromoptions = TaoSetFromOptions_BQNK; 160 tao->ops->destroy = TaoDestroy_BQNK; 161 tao->ops->view = TaoView_BQNK; 162 tao->ops->setup = TaoSetUp_BQNK; 163 164 bnk = (TAO_BNK *)tao->data; 165 bnk->computehessian = TaoBQNKComputeHessian; 166 bnk->computestep = TaoBQNKComputeStep; 167 bnk->init_type = BNK_INIT_DIRECTION; 168 169 ierr = PetscNewLog(tao,&bqnk);CHKERRQ(ierr); 170 bnk->ctx = (void*)bqnk; 171 bqnk->no_scale = PETSC_FALSE; 172 173 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);CHKERRQ(ierr); 174 ierr = PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);CHKERRQ(ierr); 175 ierr = MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");CHKERRQ(ierr); 176 ierr = MatSetType(bqnk->B, MATLMVMSR1);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 }