1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
<center><a href="snesj.c">Actual source code: snesj.c</a></center><br>
<html>
<head>
<title></title>
<meta name="generator" content="c2html 0.9.5">
<meta name="date" content="2010-04-08T19:42:47+00:00">
</head>
<body bgcolor="#FFFFFF">
<pre width="80"><a name="line1"> 1: </a><strong><font color="#228B22">#define PETSCSNES_DLL</font></strong>
<a name="line3"> 3: </a> #include <A href="../../../include/private/snesimpl.h.html">private/snesimpl.h</A>
<a name="line7"> 7: </a><font color="#B22222">/*@C</font>
<a name="line8"> 8: </a><font color="#B22222"> <A href="../../../docs/manualpages/SNES/SNESDefaultComputeJacobian.html#SNESDefaultComputeJacobian">SNESDefaultComputeJacobian</A> - Computes the Jacobian using finite differences. </font>
<a name="line10"> 10: </a><font color="#B22222"> Collective on <A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A></font>
<a name="line12"> 12: </a><font color="#B22222"> Input Parameters:</font>
<a name="line13"> 13: </a><font color="#B22222">+ x1 - compute Jacobian at this point</font>
<a name="line14"> 14: </a><font color="#B22222">- ctx - application's function context, as set with <A href="../../../docs/manualpages/SNES/SNESSetFunction.html#SNESSetFunction">SNESSetFunction</A>()</font>
<a name="line16"> 16: </a><font color="#B22222"> Output Parameters:</font>
<a name="line17"> 17: </a><font color="#B22222">+ J - Jacobian matrix (not altered in this routine)</font>
<a name="line18"> 18: </a><font color="#B22222">. B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)</font>
<a name="line19"> 19: </a><font color="#B22222">- flag - flag indicating whether the matrix sparsity structure has changed</font>
<a name="line21"> 21: </a><font color="#B22222"> Options Database Key:</font>
<a name="line22"> 22: </a><font color="#B22222">+ -snes_fd - Activates <A href="../../../docs/manualpages/SNES/SNESDefaultComputeJacobian.html#SNESDefaultComputeJacobian">SNESDefaultComputeJacobian</A>()</font>
<a name="line23"> 23: </a><font color="#B22222">. -snes_test_err - Square root of function error tolerance, default square root of machine</font>
<a name="line24"> 24: </a><font color="#B22222"> epsilon (1.e-8 in double, 3.e-4 in single)</font>
<a name="line25"> 25: </a><font color="#B22222">- -mat_fd_type - Either wp or ds (see <A href="../../../docs/manualpages/Mat/MATMFFD_WP.html#MATMFFD_WP">MATMFFD_WP</A> or <A href="../../../docs/manualpages/Mat/MATMFFD_DS.html#MATMFFD_DS">MATMFFD_DS</A>)</font>
<a name="line27"> 27: </a><font color="#B22222"> Notes:</font>
<a name="line28"> 28: </a><font color="#B22222"> This routine is slow and expensive, and is not currently optimized</font>
<a name="line29"> 29: </a><font color="#B22222"> to take advantage of sparsity in the problem. Although</font>
<a name="line30"> 30: </a><font color="#B22222"> <A href="../../../docs/manualpages/SNES/SNESDefaultComputeJacobian.html#SNESDefaultComputeJacobian">SNESDefaultComputeJacobian</A>() is not recommended for general use</font>
<a name="line31"> 31: </a><font color="#B22222"> in large-scale applications, It can be useful in checking the</font>
<a name="line32"> 32: </a><font color="#B22222"> correctness of a user-provided Jacobian.</font>
<a name="line34"> 34: </a><font color="#B22222"> An alternative routine that uses coloring to exploit matrix sparsity is</font>
<a name="line35"> 35: </a><font color="#B22222"> <A href="../../../docs/manualpages/SNES/SNESDefaultComputeJacobianColor.html#SNESDefaultComputeJacobianColor">SNESDefaultComputeJacobianColor</A>().</font>
<a name="line37"> 37: </a><font color="#B22222"> Level: intermediate</font>
<a name="line39"> 39: </a><font color="#B22222">.keywords: <A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>, finite differences, Jacobian</font>
<a name="line41"> 41: </a><font color="#B22222">.seealso: <A href="../../../docs/manualpages/SNES/SNESSetJacobian.html#SNESSetJacobian">SNESSetJacobian</A>(), <A href="../../../docs/manualpages/SNES/SNESDefaultComputeJacobianColor.html#SNESDefaultComputeJacobianColor">SNESDefaultComputeJacobianColor</A>(), <A href="../../../docs/manualpages/SNES/MatCreateSNESMF.html#MatCreateSNESMF">MatCreateSNESMF</A>()</font>
<a name="line42"> 42: </a><font color="#B22222">@*/</font>
<a name="line43"> 43: </a><strong><font color="#4169E1"><a name="SNESDefaultComputeJacobian"></a><A href="../../../docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode">PetscErrorCode</A> <A href="../../../docs/manualpages/SNES/SNESDefaultComputeJacobian.html#SNESDefaultComputeJacobian">SNESDefaultComputeJacobian</A>(<A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> snes,<A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> x1,<A href="../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A> *J,<A href="../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A> *B,<A href="../../../docs/manualpages/Mat/MatStructure.html#MatStructure">MatStructure</A> *flag,void *ctx)</font></strong>
<a name="line44"> 44: </a>{
<a name="line45"> 45: </a> <A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> j1a,j2a,x2;
<a name="line47"> 47: </a> <A href="../../../docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</A> i,N,start,end,j,value,root;
<a name="line48"> 48: </a> <A href="../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A> dx,*y,*xx,wscale;
<a name="line49"> 49: </a> <A href="../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A> amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
<a name="line50"> 50: </a> <A href="../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A> dx_min = 1.e-16,dx_par = 1.e-1,unorm;
<a name="line51"> 51: </a> <A href="../../../docs/manualpages/Sys/MPI_Comm.html#MPI_Comm">MPI_Comm</A> comm;
<a name="line52"> 52: </a> <A href="../../../docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode">PetscErrorCode</A> (*eval_fct)(<A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>,<A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>,<A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>)=0;
<a name="line53"> 53: </a> <A href="../../../docs/manualpages/Sys/PetscTruth.html#PetscTruth">PetscTruth</A> assembled,use_wp = <A href="../../../docs/manualpages/Sys/PETSC_TRUE.html#PETSC_TRUE">PETSC_TRUE</A>,flg;
<a name="line54"> 54: </a> const char *list[2] = {<font color="#666666">"ds"</font>,<font color="#666666">"wp"</font>};
<a name="line55"> 55: </a> <A href="../../../docs/manualpages/Sys/PetscMPIInt.html#PetscMPIInt">PetscMPIInt</A> size;
<a name="line56"> 56: </a> const <A href="../../../docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</A> *ranges;
<a name="line59"> 59: </a> <A href="../../../docs/manualpages/Sys/PetscOptionsGetReal.html#PetscOptionsGetReal">PetscOptionsGetReal</A>(((<A href="../../../docs/manualpages/Sys/PetscObject.html#PetscObject">PetscObject</A>)snes)->prefix,<font color="#666666">"-snes_test_err"</font>,&epsilon,0);
<a name="line60"> 60: </a> eval_fct = <A href="../../../docs/manualpages/SNES/SNESComputeFunction.html#SNESComputeFunction">SNESComputeFunction</A>;
<a name="line62"> 62: </a> <A href="../../../docs/manualpages/Sys/PetscObjectGetComm.html#PetscObjectGetComm">PetscObjectGetComm</A>((<A href="../../../docs/manualpages/Sys/PetscObject.html#PetscObject">PetscObject</A>)x1,&comm);
<a name="line63"> 63: </a> <A href="http://www.mcs.anl.gov/mpi/www/www3/MPI_Comm_size.html#MPI_Comm_size">MPI_Comm_size</A>(comm,&size);
<a name="line64"> 64: </a> <A href="../../../docs/manualpages/Mat/MatAssembled.html#MatAssembled">MatAssembled</A>(*B,&assembled);
<a name="line65"> 65: </a> <font color="#4169E1">if</font> (assembled) {
<a name="line66"> 66: </a> <A href="../../../docs/manualpages/Mat/MatZeroEntries.html#MatZeroEntries">MatZeroEntries</A>(*B);
<a name="line67"> 67: </a> }
<a name="line68"> 68: </a> <font color="#4169E1">if</font> (!snes->nvwork) {
<a name="line69"> 69: </a> snes->nvwork = 3;
<a name="line70"> 70: </a> <A href="../../../docs/manualpages/Vec/VecDuplicateVecs.html#VecDuplicateVecs">VecDuplicateVecs</A>(x1,snes->nvwork,&snes->vwork);
<a name="line71"> 71: </a> PetscLogObjectParents(snes,snes->nvwork,snes->vwork);
<a name="line72"> 72: </a> }
<a name="line73"> 73: </a> j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2];
<a name="line75"> 75: </a> <A href="../../../docs/manualpages/Vec/VecGetSize.html#VecGetSize">VecGetSize</A>(x1,&N);
<a name="line76"> 76: </a> <A href="../../../docs/manualpages/Vec/VecGetOwnershipRange.html#VecGetOwnershipRange">VecGetOwnershipRange</A>(x1,&start,&end);
<a name="line77"> 77: </a> (*eval_fct)(snes,x1,j1a);
<a name="line79"> 79: </a> <A href="../../../docs/manualpages/Sys/PetscOptionsEList.html#PetscOptionsEList">PetscOptionsEList</A>(<font color="#666666">"-mat_fd_type"</font>,<font color="#666666">"Algorithm to compute difference parameter"</font>,<font color="#666666">"<A href="../../../docs/manualpages/SNES/SNESDefaultComputeJacobian.html#SNESDefaultComputeJacobian">SNESDefaultComputeJacobian</A>"</font>,list,2,<font color="#666666">"wp"</font>,&value,&flg);
<a name="line80"> 80: </a> <font color="#4169E1">if</font> (flg && !value) {
<a name="line81"> 81: </a> use_wp = <A href="../../../docs/manualpages/Sys/PETSC_FALSE.html#PETSC_FALSE">PETSC_FALSE</A>;
<a name="line82"> 82: </a> }
<a name="line83"> 83: </a> <font color="#4169E1">if</font> (use_wp) {
<a name="line84"> 84: </a> <A href="../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(x1,<A href="../../../docs/manualpages/Vec/NORM_2.html#NORM_2">NORM_2</A>,&unorm);
<a name="line85"> 85: </a> }
<a name="line86"> 86: </a> <font color="#B22222">/* Compute Jacobian approximation, 1 column at a time. </font>
<a name="line87"> 87: </a><font color="#B22222"> x1 = current iterate, j1a = F(x1)</font>
<a name="line88"> 88: </a><font color="#B22222"> x2 = perturbed iterate, j2a = F(x2)</font>
<a name="line89"> 89: </a><font color="#B22222"> */</font>
<a name="line90"> 90: </a> <font color="#4169E1">for</font> (i=0; i<N; i++) {
<a name="line91"> 91: </a> <A href="../../../docs/manualpages/Vec/VecCopy.html#VecCopy">VecCopy</A>(x1,x2);
<a name="line92"> 92: </a> <font color="#4169E1">if</font> (i>= start && i<end) {
<a name="line93"> 93: </a> <A href="../../../docs/manualpages/Vec/VecGetArray.html#VecGetArray">VecGetArray</A>(x1,&xx);
<a name="line94"> 94: </a> <font color="#4169E1">if</font> (use_wp) {
<a name="line95"> 95: </a> dx = 1.0 + unorm;
<a name="line96"> 96: </a> } <font color="#4169E1">else</font> {
<a name="line97"> 97: </a> dx = xx[i-start];
<a name="line98"> 98: </a> }
<a name="line99"> 99: </a> <A href="../../../docs/manualpages/Vec/VecRestoreArray.html#VecRestoreArray">VecRestoreArray</A>(x1,&xx);
<a name="line100">100: </a><font color="#A020F0">#if !defined(PETSC_USE_COMPLEX)</font>
<a name="line101">101: </a> <font color="#4169E1">if</font> (dx < dx_min && dx >= 0.0) dx = dx_par;
<a name="line102">102: </a> <font color="#4169E1">else</font> <font color="#4169E1">if</font> (dx < 0.0 && dx > -dx_min) dx = -dx_par;
<a name="line103">103: </a><font color="#A020F0">#else</font>
<a name="line104">104: </a> <font color="#4169E1">if</font> (PetscAbsScalar(dx) < dx_min && PetscRealPart(dx) >= 0.0) dx = dx_par;
<a name="line105">105: </a> <font color="#4169E1">else</font> <font color="#4169E1">if</font> (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < dx_min) dx = -dx_par;
<a name="line106">106: </a><font color="#A020F0">#endif</font>
<a name="line107">107: </a> dx *= epsilon;
<a name="line108">108: </a> wscale = 1.0/dx;
<a name="line109">109: </a> <A href="../../../docs/manualpages/Vec/VecSetValues.html#VecSetValues">VecSetValues</A>(x2,1,&i,&dx,<A href="../../../docs/manualpages/Sys/ADD_VALUES.html#ADD_VALUES">ADD_VALUES</A>);
<a name="line110">110: </a> } <font color="#4169E1">else</font> {
<a name="line111">111: </a> wscale = 0.0;
<a name="line112">112: </a> }
<a name="line113">113: </a> (*eval_fct)(snes,x2,j2a);
<a name="line114">114: </a> <A href="../../../docs/manualpages/Vec/VecAXPY.html#VecAXPY">VecAXPY</A>(j2a,-1.0,j1a);
<a name="line115">115: </a> <font color="#B22222">/* Communicate scale=1/dx_i to all processors */</font>
<a name="line116">116: </a> <A href="../../../docs/manualpages/Vec/VecGetOwnershipRanges.html#VecGetOwnershipRanges">VecGetOwnershipRanges</A>(x1,&ranges);
<a name="line117">117: </a> root = size;
<a name="line118">118: </a> <font color="#4169E1">for</font> (j=size-1; j>-1; j--){
<a name="line119">119: </a> root--;
<a name="line120">120: </a> <font color="#4169E1">if</font> (i>=ranges[j]) <font color="#4169E1">break</font>;
<a name="line121">121: </a> }
<a name="line122">122: </a> <A href="http://www.mcs.anl.gov/mpi/www/www3/MPI_Bcast.html#MPI_Bcast">MPI_Bcast</A>(&wscale,1,<A href="../../../docs/manualpages/Sys/MPIU_SCALAR.html#MPIU_SCALAR">MPIU_SCALAR</A>,root,comm);
<a name="line124">124: </a> <A href="../../../docs/manualpages/Vec/VecScale.html#VecScale">VecScale</A>(j2a,wscale);
<a name="line125">125: </a> <A href="../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(j2a,<A href="../../../docs/manualpages/Vec/NORM_INFINITY.html#NORM_INFINITY">NORM_INFINITY</A>,&amax); amax *= 1.e-14;
<a name="line126">126: </a> <A href="../../../docs/manualpages/Vec/VecGetArray.html#VecGetArray">VecGetArray</A>(j2a,&y);
<a name="line127">127: </a> <font color="#4169E1">for</font> (j=start; j<end; j++) {
<a name="line128">128: </a> <font color="#4169E1">if</font> (PetscAbsScalar(y[j-start]) > amax || j == i) {
<a name="line129">129: </a> <A href="../../../docs/manualpages/Mat/MatSetValues.html#MatSetValues">MatSetValues</A>(*B,1,&j,1,&i,y+j-start,<A href="../../../docs/manualpages/Sys/INSERT_VALUES.html#INSERT_VALUES">INSERT_VALUES</A>);
<a name="line130">130: </a> }
<a name="line131">131: </a> }
<a name="line132">132: </a> <A href="../../../docs/manualpages/Vec/VecRestoreArray.html#VecRestoreArray">VecRestoreArray</A>(j2a,&y);
<a name="line133">133: </a> }
<a name="line134">134: </a> <A href="../../../docs/manualpages/Mat/MatAssemblyBegin.html#MatAssemblyBegin">MatAssemblyBegin</A>(*B,MAT_FINAL_ASSEMBLY);
<a name="line135">135: </a> <A href="../../../docs/manualpages/Mat/MatAssemblyEnd.html#MatAssemblyEnd">MatAssemblyEnd</A>(*B,MAT_FINAL_ASSEMBLY);
<a name="line136">136: </a> <font color="#4169E1">if</font> (*B != *J) {
<a name="line137">137: </a> <A href="../../../docs/manualpages/Mat/MatAssemblyBegin.html#MatAssemblyBegin">MatAssemblyBegin</A>(*J,MAT_FINAL_ASSEMBLY);
<a name="line138">138: </a> <A href="../../../docs/manualpages/Mat/MatAssemblyEnd.html#MatAssemblyEnd">MatAssemblyEnd</A>(*J,MAT_FINAL_ASSEMBLY);
<a name="line139">139: </a> }
<a name="line140">140: </a> *flag = DIFFERENT_NONZERO_PATTERN;
<a name="line141">141: </a> <font color="#4169E1">return</font>(0);
<a name="line142">142: </a>}
</pre>
</body>
</html>
|