1 #include <petsctao.h> /*I "petsctao.h" I*/ 2 #include <petsc/private/taoimpl.h> 3 #include <petscsnes.h> 4 #include <petscdmshell.h> 5 6 /* 7 For finited difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault 8 */ 9 static PetscErrorCode Fsnes(SNES snes,Vec X,Vec G,void* ctx) 10 { 11 PetscErrorCode ierr; 12 Tao tao = (Tao)ctx; 13 14 PetscFunctionBegin; 15 PetscValidHeaderSpecific(tao,TAO_CLASSID,4); 16 ierr = TaoComputeGradient(tao,X,G);CHKERRQ(ierr); 17 PetscFunctionReturn(0); 18 } 19 20 /*@C 21 TaoDefaultComputeGradient - computes the gradient using finite differences. 22 23 Collective on Tao 24 25 Input Parameters: 26 + tao - the Tao context 27 . X - compute gradient at this point 28 - dummy - not used 29 30 Output Parameters: 31 . G - Gradient Vector 32 33 Options Database Key: 34 + -tao_fd_gradient - activates TaoDefaultComputeGradient() 35 - -tao_fd_delta <delta> - change in X used to calculate finite differences 36 37 Level: advanced 38 39 Notes: 40 This routine is slow and expensive, and is not currently optimized 41 to take advantage of sparsity in the problem. Although 42 TaoDefaultComputeGradient is not recommended for general use 43 in large-scale applications, It can be useful in checking the 44 correctness of a user-provided gradient. Use the tao method TAOTEST 45 to get an indication of whether your gradient is correct. 46 This finite difference gradient evaluation can be set using the routine TaoSetGradientRoutine() or by using the command line option -tao_fd_gradient 47 48 .seealso: TaoSetGradientRoutine() 49 @*/ 50 PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec Xin,Vec G,void *dummy) 51 { 52 Vec X; 53 PetscScalar *g; 54 PetscReal f, f2; 55 PetscErrorCode ierr; 56 PetscInt low,high,N,i; 57 PetscBool flg; 58 PetscReal h=.5*PETSC_SQRT_MACHINE_EPSILON; 59 60 PetscFunctionBegin; 61 ierr = PetscOptionsGetReal(((PetscObject)tao)->options,((PetscObject)tao)->prefix,"-tao_fd_delta",&h,&flg);CHKERRQ(ierr); 62 ierr = VecDuplicate(Xin,&X);CHKERRQ(ierr); 63 ierr = VecCopy(Xin,X);CHKERRQ(ierr); 64 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 65 ierr = VecGetOwnershipRange(X,&low,&high);CHKERRQ(ierr); 66 ierr = VecSetOption(X,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 67 ierr = VecGetArray(G,&g);CHKERRQ(ierr); 68 for (i=0;i<N;i++) { 69 ierr = VecSetValue(X,i,-h,ADD_VALUES);CHKERRQ(ierr); 70 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); 71 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); 72 ierr = TaoComputeObjective(tao,X,&f);CHKERRQ(ierr); 73 ierr = VecSetValue(X,i,2.0*h,ADD_VALUES);CHKERRQ(ierr); 74 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); 75 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); 76 ierr = TaoComputeObjective(tao,X,&f2);CHKERRQ(ierr); 77 ierr = VecSetValue(X,i,-h,ADD_VALUES);CHKERRQ(ierr); 78 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); 79 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); 80 if (i>=low && i<high) { 81 g[i-low]=(f2-f)/(2.0*h); 82 } 83 } 84 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr); 85 ierr = VecDestroy(&X);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 /*@C 90 TaoDefaultComputeHessian - Computes the Hessian using finite differences. 91 92 Collective on Tao 93 94 Input Parameters: 95 + tao - the Tao context 96 . V - compute Hessian at this point 97 - dummy - not used 98 99 Output Parameters: 100 + H - Hessian matrix (not altered in this routine) 101 - B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 102 103 Options Database Key: 104 . -tao_fd_hessian - activates TaoDefaultComputeHessian() 105 106 Level: advanced 107 108 Notes: 109 This routine is slow and expensive, and is not currently optimized 110 to take advantage of sparsity in the problem. Although 111 TaoDefaultComputeHessian() is not recommended for general use 112 in large-scale applications, It can be useful in checking the 113 correctness of a user-provided Hessian. 114 115 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradientRoutine(), TaoDefaultComputeGradient() 116 @*/ 117 PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void *dummy) 118 { 119 PetscErrorCode ierr; 120 SNES snes; 121 DM dm; 122 123 PetscFunctionBegin; 124 ierr = PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n");CHKERRQ(ierr); 125 ierr = SNESCreate(PetscObjectComm((PetscObject)H),&snes);CHKERRQ(ierr); 126 ierr = SNESSetFunction(snes,NULL,Fsnes,tao);CHKERRQ(ierr); 127 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 128 ierr = DMShellSetGlobalVector(dm,V);CHKERRQ(ierr); 129 ierr = SNESSetUp(snes);CHKERRQ(ierr); 130 if (H) { 131 PetscInt n,N; 132 133 ierr = VecGetSize(V,&N);CHKERRQ(ierr); 134 ierr = VecGetLocalSize(V,&n);CHKERRQ(ierr); 135 ierr = MatSetSizes(H,n,n,N,N);CHKERRQ(ierr); 136 ierr = MatSetUp(H);CHKERRQ(ierr); 137 } 138 if (B && B != H) { 139 PetscInt n,N; 140 141 ierr = VecGetSize(V,&N);CHKERRQ(ierr); 142 ierr = VecGetLocalSize(V,&n);CHKERRQ(ierr); 143 ierr = MatSetSizes(B,n,n,N,N);CHKERRQ(ierr); 144 ierr = MatSetUp(B);CHKERRQ(ierr); 145 } 146 ierr = SNESComputeJacobianDefault(snes,V,H,B,NULL);CHKERRQ(ierr); 147 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 /*@C 152 TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 153 154 Collective on Tao 155 156 Input Parameters: 157 + tao - the Tao context 158 . V - compute Hessian at this point 159 - ctx - the PetscColoring object (must be of type MatFDColoring) 160 161 Output Parameters: 162 + H - Hessian matrix (not altered in this routine) 163 - B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 164 165 Level: advanced 166 167 168 .seealso: TaoSetHessianRoutine(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradientRoutine() 169 @*/ 170 PetscErrorCode TaoDefaultComputeHessianColor(Tao tao,Vec V,Mat H,Mat B,void *ctx) 171 { 172 PetscErrorCode ierr; 173 MatFDColoring coloring = (MatFDColoring)ctx; 174 175 PetscFunctionBegin; 176 PetscValidHeaderSpecific(ctx,MAT_FDCOLORING_CLASSID,6); 177 ierr = PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n");CHKERRQ(ierr); 178 ierr = MatFDColoringApply(B,coloring,V,ctx);CHKERRQ(ierr); 179 if (H != B) { 180 ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 181 ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 182 } 183 PetscFunctionReturn(0); 184 } 185 186 PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao,Vec X,Mat H,Mat B,void *ctx) 187 { 188 PetscInt n,N; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 if (B && B != H) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Preconditioning Hessian matrix"); 193 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 194 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 195 ierr = MatSetSizes(H,n,n,N,N);CHKERRQ(ierr); 196 ierr = MatSetType(H,MATMFFD);CHKERRQ(ierr); 197 ierr = MatSetUp(H);CHKERRQ(ierr); 198 ierr = MatMFFDSetBase(H,X,NULL);CHKERRQ(ierr); 199 ierr = MatMFFDSetFunction(H,(PetscErrorCode (*)(void*,Vec,Vec))TaoComputeGradient,tao);CHKERRQ(ierr); 200 ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 201 ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204