1 #include <../src/tao/bound/impls/bqnk/bqnk.h> /*I "petsctao.h" I*/ /*I "petscmat.h" I*/ 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 PetscReal gnorm2, delta; 10 11 PetscFunctionBegin; 12 /* Alias the LMVM matrix into the TAO hessian */ 13 if (tao->hessian) { 14 ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr); 15 } 16 if (tao->hessian_pre) { 17 ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr); 18 } 19 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 20 tao->hessian = bqnk->B; 21 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 22 tao->hessian_pre = bqnk->B; 23 /* Update the Hessian with the latest solution */ 24 if (bqnk->is_spd) { 25 gnorm2 = bnk->gnorm*bnk->gnorm; 26 if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON; 27 if (bnk->f == 0.0) { 28 delta = 2.0 / gnorm2; 29 } else { 30 delta = 2.0 * PetscAbsScalar(bnk->f) / gnorm2; 31 } 32 ierr = MatLMVMSymBroydenSetDelta(bqnk->B, delta);CHKERRQ(ierr); 33 } 34 ierr = MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 35 ierr = MatLMVMResetShift(tao->hessian);CHKERRQ(ierr); 36 /* Prepare the reduced sub-matrices for the inactive set */ 37 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 38 if (bnk->active_idx) { 39 ierr = MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);CHKERRQ(ierr); 40 ierr = PCLMVMSetIS(bqnk->pc, bnk->inactive_idx);CHKERRQ(ierr); 41 } else { 42 ierr = PetscObjectReference((PetscObject)tao->hessian);CHKERRQ(ierr); 43 bnk->H_inactive = tao->hessian; 44 ierr = PCLMVMClearIS(bqnk->pc);CHKERRQ(ierr); 45 } 46 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 47 ierr = PetscObjectReference((PetscObject)bnk->H_inactive);CHKERRQ(ierr); 48 bnk->Hpre_inactive = bnk->H_inactive; 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type) 53 { 54 TAO_BNK *bnk = (TAO_BNK *)tao->data; 55 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = TaoBNKComputeStep(tao, shift, ksp_reason, step_type);CHKERRQ(ierr); 60 if (*ksp_reason < 0) { 61 /* Krylov solver failed to converge so reset the LMVM matrix */ 62 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 63 ierr = MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 64 } 65 PetscFunctionReturn(0); 66 } 67 68 PetscErrorCode TaoSolve_BQNK(Tao tao) 69 { 70 TAO_BNK *bnk = (TAO_BNK *)tao->data; 71 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 72 Mat_LMVM *lmvm = (Mat_LMVM*)bqnk->B->data; 73 Mat_LMVM *J0; 74 Mat_SymBrdn *diag_ctx; 75 PetscBool flg = PETSC_FALSE; 76 PetscErrorCode ierr; 77 78 PetscFunctionBegin; 79 if (!tao->recycle) { 80 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 81 lmvm->nresets = 0; 82 if (lmvm->J0) { 83 ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &flg);CHKERRQ(ierr); 84 if (flg) { 85 J0 = (Mat_LMVM*)lmvm->J0->data; 86 J0->nresets = 0; 87 } 88 } 89 flg = PETSC_FALSE; 90 ierr = PetscObjectTypeCompareAny((PetscObject)bqnk->B, &flg, MATLMVMSYMBROYDEN, MATLMVMSYMBADBROYDEN, MATLMVMBFGS, MATLMVMDFP, "");CHKERRQ(ierr); 91 if (flg) { 92 diag_ctx = (Mat_SymBrdn*)lmvm->ctx; 93 J0 = (Mat_LMVM*)diag_ctx->D->data; 94 J0->nresets = 0; 95 } 96 } 97 ierr = (*bqnk->solve)(tao);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode TaoSetUp_BQNK(Tao tao) 102 { 103 TAO_BNK *bnk = (TAO_BNK *)tao->data; 104 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 105 PetscErrorCode ierr; 106 PetscInt n, N; 107 PetscBool is_lmvm, is_sym, is_spd; 108 109 PetscFunctionBegin; 110 ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr); 111 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 112 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 113 ierr = MatSetSizes(bqnk->B, n, n, N, N);CHKERRQ(ierr); 114 ierr = MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);CHKERRQ(ierr); 115 ierr = PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);CHKERRQ(ierr); 116 if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type"); 117 ierr = MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);CHKERRQ(ierr); 118 if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric"); 119 ierr = MatGetOption(bqnk->B, MAT_SPD, &is_spd);CHKERRQ(ierr); 120 ierr = KSPGetPC(tao->ksp, &bqnk->pc);CHKERRQ(ierr); 121 ierr = PCSetType(bqnk->pc, PCLMVM);CHKERRQ(ierr); 122 ierr = PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject,Tao tao) 127 { 128 TAO_BNK *bnk = (TAO_BNK *)tao->data; 129 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = TaoSetFromOptions_BNK(PetscOptionsObject,tao);CHKERRQ(ierr); 134 ierr = TaoSetFromOptions(bnk->bncg);CHKERRQ(ierr); 135 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 136 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 137 if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION; 138 ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr); 139 ierr = MatGetOption(bqnk->B, MAT_SPD, &bqnk->is_spd);CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer) 144 { 145 TAO_BNK *bnk = (TAO_BNK*)tao->data; 146 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 147 PetscErrorCode ierr; 148 PetscBool isascii; 149 150 PetscFunctionBegin; 151 ierr = TaoView_BNK(tao, viewer);CHKERRQ(ierr); 152 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 153 if (isascii) { 154 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 155 ierr = MatView(bqnk->B, viewer);CHKERRQ(ierr); 156 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 157 } 158 PetscFunctionReturn(0); 159 } 160 161 static PetscErrorCode TaoDestroy_BQNK(Tao tao) 162 { 163 TAO_BNK *bnk = (TAO_BNK*)tao->data; 164 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 165 PetscErrorCode ierr; 166 167 PetscFunctionBegin; 168 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 169 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 170 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 171 ierr = PetscFree(bnk->ctx);CHKERRQ(ierr); 172 ierr = TaoDestroy_BNK(tao);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao) 177 { 178 TAO_BNK *bnk; 179 TAO_BQNK *bqnk; 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 184 ierr = KSPSetOptionsPrefix(tao->ksp, "tao_bqnk_");CHKERRQ(ierr); 185 tao->ops->solve = TaoSolve_BQNK; 186 tao->ops->setfromoptions = TaoSetFromOptions_BQNK; 187 tao->ops->destroy = TaoDestroy_BQNK; 188 tao->ops->view = TaoView_BQNK; 189 tao->ops->setup = TaoSetUp_BQNK; 190 191 bnk = (TAO_BNK *)tao->data; 192 bnk->computehessian = TaoBQNKComputeHessian; 193 bnk->computestep = TaoBQNKComputeStep; 194 bnk->init_type = BNK_INIT_DIRECTION; 195 196 ierr = PetscNewLog(tao,&bqnk);CHKERRQ(ierr); 197 bnk->ctx = (void*)bqnk; 198 bqnk->is_spd = PETSC_TRUE; 199 200 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);CHKERRQ(ierr); 201 ierr = PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);CHKERRQ(ierr); 202 ierr = MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");CHKERRQ(ierr); 203 ierr = MatSetType(bqnk->B, MATLMVMSR1);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 /*@ 208 TaoGetLMVMMatrix - Returns a pointer to the internal LMVM matrix. Valid 209 only for quasi-Newton family of methods. 210 211 Input Parameters: 212 . tao - Tao solver context 213 214 Output Parameters: 215 . B - LMVM matrix 216 217 Level: advanced 218 219 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoSetLMVMMatrix() 220 @*/ 221 PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B) 222 { 223 TAO_BNK *bnk = (TAO_BNK*)tao->data; 224 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 225 PetscErrorCode ierr; 226 PetscBool flg = PETSC_FALSE; 227 228 PetscFunctionBegin; 229 ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr); 230 if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms"); 231 *B = bqnk->B; 232 PetscFunctionReturn(0); 233 } 234 235 /*@ 236 TaoSetLMVMMatrix - Sets an external LMVM matrix into the Tao solver. Valid 237 only for quasi-Newton family of methods. 238 239 QN family of methods create their own LMVM matrices and users who wish to 240 manipulate this matrix should use TaoGetLMVMMatrix() instead. 241 242 Input Parameters: 243 + tao - Tao solver context 244 - B - LMVM matrix 245 246 Level: advanced 247 248 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoGetLMVMMatrix() 249 @*/ 250 PetscErrorCode TaoSetLMVMMatrix(Tao tao, Mat B) 251 { 252 TAO_BNK *bnk = (TAO_BNK*)tao->data; 253 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 254 PetscErrorCode ierr; 255 PetscBool flg = PETSC_FALSE; 256 257 PetscFunctionBegin; 258 ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr); 259 if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms"); 260 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &flg);CHKERRQ(ierr); 261 if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given matrix is not an LMVM matrix"); 262 if (bqnk->B) { 263 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 264 } 265 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 266 bqnk->B = B; 267 PetscFunctionReturn(0); 268 } 269