File: | snes/interface/snesj.c |
Warning: | line 79, column 10 3rd function call argument is an uninitialized value |
[?] Use j/k keys for keyboard navigation
1 | ||||
2 | #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ | |||
3 | #include <petscdm.h> | |||
4 | ||||
5 | /*@C | |||
6 | SNESComputeJacobianDefault - Computes the Jacobian using finite differences. | |||
7 | ||||
8 | Collective on SNES | |||
9 | ||||
10 | Input Parameters: | |||
11 | + x1 - compute Jacobian at this point | |||
12 | - ctx - application's function context, as set with SNESSetFunction() | |||
13 | ||||
14 | Output Parameters: | |||
15 | + J - Jacobian matrix (not altered in this routine) | |||
16 | - B - newly computed Jacobian matrix to use with preconditioner (generally the same as J) | |||
17 | ||||
18 | Options Database Key: | |||
19 | + -snes_fd - Activates SNESComputeJacobianDefault() | |||
20 | . -snes_test_err - Square root of function error tolerance, default square root of machine | |||
21 | epsilon (1.e-8 in double, 3.e-4 in single) | |||
22 | - -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS) | |||
23 | ||||
24 | Notes: | |||
25 | This routine is slow and expensive, and is not currently optimized | |||
26 | to take advantage of sparsity in the problem. Although | |||
27 | SNESComputeJacobianDefault() is not recommended for general use | |||
28 | in large-scale applications, It can be useful in checking the | |||
29 | correctness of a user-provided Jacobian. | |||
30 | ||||
31 | An alternative routine that uses coloring to exploit matrix sparsity is | |||
32 | SNESComputeJacobianDefaultColor(). | |||
33 | ||||
34 | Level: intermediate | |||
35 | ||||
36 | .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF() | |||
37 | @*/ | |||
38 | PetscErrorCode SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx) | |||
39 | { | |||
40 | Vec j1a,j2a,x2; | |||
| ||||
41 | PetscErrorCode ierr; | |||
42 | PetscInt i,N,start,end,j,value,root; | |||
43 | PetscScalar dx,*y,wscale; | |||
44 | const PetscScalar *xx; | |||
45 | PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON1.490116119384766e-08; | |||
46 | PetscReal dx_min = 1.e-16,dx_par = 1.e-1,unorm; | |||
47 | MPI_Comm comm; | |||
48 | PetscBool assembled,use_wp = PETSC_TRUE,flg; | |||
49 | const char *list[2] = {"ds","wp"}; | |||
50 | PetscMPIInt size; | |||
51 | const PetscInt *ranges; | |||
52 | ||||
53 | PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize < 64)) { petscstack->function[petscstack->currentsize ] = __func__; petscstack->file[petscstack->currentsize] = "/sandbox/petsc/petsc.next/src/snes/interface/snesj.c"; petscstack ->line[petscstack->currentsize] = 53; petscstack->petscroutine [petscstack->currentsize] = PETSC_TRUE; petscstack->currentsize ++; } if (petscstack) { petscstack->hotdepth += (PETSC_FALSE || petscstack->hotdepth); } ; } while (0); ; } while (0); | |||
54 | /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */ | |||
55 | ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),55,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
56 | ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),56,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
57 | ||||
58 | ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),58,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
59 | ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),59,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
60 | ierr = MatAssembled(B,&assembled);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),60,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
61 | if (assembled) { | |||
62 | ierr = MatZeroEntries(B);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),62,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
63 | } | |||
64 | if (!snes->nvwork) { | |||
65 | if (snes->dm) { | |||
66 | ierr = DMGetGlobalVector(snes->dm,&j1a);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),66,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
67 | ierr = DMGetGlobalVector(snes->dm,&j2a);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),67,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
68 | ierr = DMGetGlobalVector(snes->dm,&x2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),68,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
69 | } else { | |||
70 | snes->nvwork = 3; | |||
71 | ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),71,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
72 | ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork)0;do{int _i; for (_i=0; _i<(snes->nvwork); _i++) {ierr = PetscLogObjectParent((PetscObject)(snes),(PetscObject)(snes-> vwork)[_i]);do {if (__builtin_expect(!!(ierr),0)) return PetscError (((MPI_Comm)0x44000001),72,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0);}}while(0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),72,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
73 | j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; | |||
74 | } | |||
75 | } | |||
76 | ||||
77 | ierr = VecGetSize(x1,&N);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),77,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
78 | ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),78,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
79 | ierr = SNESComputeFunction(snes,x1,j1a);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),79,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
| ||||
80 | ||||
81 | ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES")0; do { PetscOptionItems PetscOptionsObjectBase; PetscOptionItems *PetscOptionsObject = &PetscOptionsObjectBase; PetscMemzero (PetscOptionsObject,sizeof(PetscOptionItems)); for (PetscOptionsObject ->count=(PetscOptionsPublish?-1:1); PetscOptionsObject-> count<2; PetscOptionsObject->count++) { PetscErrorCode _5_ierr = PetscOptionsBegin_Private(PetscOptionsObject,PetscObjectComm ((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options" ,"SNES");do {if (__builtin_expect(!!(_5_ierr),0)) return PetscError (((MPI_Comm)0x44000001),81,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),81,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
82 | ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg)PetscOptionsEList_Private(PetscOptionsObject,"-mat_fd_type","Algorithm to compute difference parameter" ,"SNESComputeJacobianDefault",list,2,"wp",&value,&flg );CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),82,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
83 | ierr = PetscOptionsEnd()_5_ierr = PetscOptionsEnd_Private(PetscOptionsObject);do {if ( __builtin_expect(!!(_5_ierr),0)) return PetscError(((MPI_Comm )0x44000001),83,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);}} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),83,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
84 | if (flg && !value) use_wp = PETSC_FALSE; | |||
85 | ||||
86 | if (use_wp) { | |||
87 | ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),87,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
88 | } | |||
89 | /* Compute Jacobian approximation, 1 column at a time. | |||
90 | x1 = current iterate, j1a = F(x1) | |||
91 | x2 = perturbed iterate, j2a = F(x2) | |||
92 | */ | |||
93 | for (i=0; i<N; i++) { | |||
94 | ierr = VecCopy(x1,x2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),94,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
95 | if (i>= start && i<end) { | |||
96 | ierr = VecGetArrayRead(x1,&xx);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),96,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
97 | if (use_wp) dx = PetscSqrtReal(1.0 + unorm)sqrt(1.0 + unorm); | |||
98 | else dx = xx[i-start]; | |||
99 | ierr = VecRestoreArrayRead(x1,&xx);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),99,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
100 | if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx)(dx) < 0. ? -1. : 1.) * dx_par; | |||
101 | dx *= epsilon; | |||
102 | wscale = 1.0/dx; | |||
103 | ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),103,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
104 | } else { | |||
105 | wscale = 0.0; | |||
106 | } | |||
107 | ierr = VecAssemblyBegin(x2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),107,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
108 | ierr = VecAssemblyEnd(x2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),108,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
109 | ierr = SNESComputeFunction(snes,x2,j2a);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),109,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
110 | ierr = VecAXPY(j2a,-1.0,j1a);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),110,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
111 | /* Communicate scale=1/dx_i to all processors */ | |||
112 | ierr = VecGetOwnershipRanges(x1,&ranges);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),112,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
113 | root = size; | |||
114 | for (j=size-1; j>-1; j--) { | |||
115 | root--; | |||
116 | if (i>=ranges[j]) break; | |||
117 | } | |||
118 | ierr = MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm)((petsc_allreduce_ct += PetscMPIParallelComm((comm)),0) || MPI_Bcast ((&wscale),(1),(((MPI_Datatype)0x4c00080b)),(root),(comm) ));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),118,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
119 | ||||
120 | ierr = VecScale(j2a,wscale);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),120,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
121 | ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),121,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); amax *= 1.e-14; | |||
122 | ierr = VecGetArray(j2a,&y);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),122,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
123 | for (j=start; j<end; j++) { | |||
124 | if (PetscAbsScalar(y[j-start]) > amax || j == i) { | |||
125 | ierr = MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),125,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
126 | } | |||
127 | } | |||
128 | ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),128,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
129 | } | |||
130 | if (snes->dm) { | |||
131 | ierr = DMRestoreGlobalVector(snes->dm,&j1a);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),131,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
132 | ierr = DMRestoreGlobalVector(snes->dm,&j2a);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),132,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
133 | ierr = DMRestoreGlobalVector(snes->dm,&x2);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),133,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
134 | } | |||
135 | ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),135,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
136 | ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),136,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
137 | if (B != J) { | |||
138 | ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),138,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
139 | ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm )0x44000001),139,__func__,"/sandbox/petsc/petsc.next/src/snes/interface/snesj.c" ,ierr,PETSC_ERROR_REPEAT," ");} while (0); | |||
140 | } | |||
141 | PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize > 0) { petscstack->currentsize--; petscstack->function [petscstack->currentsize] = 0; petscstack->file[petscstack ->currentsize] = 0; petscstack->line[petscstack->currentsize ] = 0; petscstack->petscroutine[petscstack->currentsize ] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth = (((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack-> hotdepth-1)); } ; } while (0); return(0);} while (0); | |||
142 | } | |||
143 | ||||
144 |