← Back to branch summary

~ubuntu-branches/ubuntu/warty/petsc/warty

« back to all changes in this revision

Viewing changes to src/ts/interface/tsfd.c.html

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2004-06-07 13:41:43 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040607134143-92p586zrauvie0le
Tags: 2.2.0-2
http://bugs.debian.org/249617
* Upstream patch level 2.
* New PETSC_BOPT_EXTRA option for different BOPT and lib names, with _c++
  symlinks only for plain and single (closes: #249617).
* New DEBIAN_DIST=contrib option to link with hypre, parmetis (closes:
  #249619).
* Combined petsc-c and petsc-fortran substvars into petsc-compilers.
* Extra quote in -dev prerm eliminates "too many arguments" problem.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/ts/interface/tsfd.c.html
3
3
<html>
4
4
<head>
5
5
<title></title>
6
 
<meta name="generator" content="c2html 0.9.1">
7
 
<meta name="date" content="2002-05-31T16:34:03+00:00">
 
6
<meta name="generator" content="c2html 0.9.4">
 
7
<meta name="date" content="2004-02-27T20:03:29+00:00">
8
8
</head>
9
9
 
10
10
<body bgcolor="#FFFFFF">
13
13
<a name="line3"> 3: </a> #include <A href="../../../src/mat/matimpl.h.html">src/mat/matimpl.h</A>
14
14
<a name="line4"> 4: </a> #include <A href="../../../src/ts/tsimpl.h.html">src/ts/tsimpl.h</A>
15
15
 
16
 
<a name="line6">  6: </a><font color="#B22222">/*@C</font>
17
 
<a name="line7">  7: </a><font color="#B22222">    <A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A> - Computes the Jacobian using</font>
18
 
<a name="line8">  8: </a><font color="#B22222">    finite differences and coloring to exploit matrix sparsity.  </font>
19
 
<a name="line9">  9: </a><font color="#B22222">  </font>
20
 
<a name="line10"> 10: </a><font color="#B22222">    Collective on <A href="../../../docs/manualpages/TS/TS.html#TS">TS</A>, <A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> and <A href="../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A></font>
21
 
 
22
 
<a name="line12"> 12: </a><font color="#B22222">    Input Parameters:</font>
23
 
<a name="line13"> 13: </a><font color="#B22222">+   ts - nonlinear solver object</font>
24
 
<a name="line14"> 14: </a><font color="#B22222">.   t - current time</font>
25
 
<a name="line15"> 15: </a><font color="#B22222">.   x1 - location at which to evaluate Jacobian</font>
26
 
<a name="line16"> 16: </a><font color="#B22222">-   ctx - coloring context, where ctx must have type <A href="../../../docs/manualpages/Mat/MatFDColoring.html#MatFDColoring">MatFDColoring</A>, </font>
27
 
<a name="line17"> 17: </a><font color="#B22222">          as created via <A href="../../../docs/manualpages/MatFD/MatFDColoringCreate.html#MatFDColoringCreate">MatFDColoringCreate</A>()</font>
28
 
 
29
 
<a name="line19"> 19: </a><font color="#B22222">    Output Parameters:</font>
30
 
<a name="line20"> 20: </a><font color="#B22222">+   J - Jacobian matrix (not altered in this routine)</font>
31
 
<a name="line21"> 21: </a><font color="#B22222">.   B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)</font>
32
 
<a name="line22"> 22: </a><font color="#B22222">-   flag - flag indicating whether the matrix sparsity structure has changed</font>
33
 
 
34
 
<a name="line24"> 24: </a><font color="#B22222">   Options Database Keys:</font>
35
 
<a name="line25"> 25: </a><font color="#B22222">$  -mat_fd_coloring_freq &lt;freq&gt; </font>
36
 
 
37
 
<a name="line27"> 27: </a><font color="#B22222">   Level: intermediate</font>
38
 
 
39
 
<a name="line29"> 29: </a><font color="#B22222">.keywords: <A href="../../../docs/manualpages/TS/TS.html#TS">TS</A>, finite differences, Jacobian, coloring, sparse</font>
40
 
 
41
 
<a name="line31"> 31: </a><font color="#B22222">.seealso: TSSetJacobian(), <A href="../../../docs/manualpages/MatFD/MatFDColoringCreate.html#MatFDColoringCreate">MatFDColoringCreate</A>(), <A href="../../../docs/manualpages/MatFD/MatFDColoringSetFunction.html#MatFDColoringSetFunction">MatFDColoringSetFunction</A>()</font>
42
 
<a name="line32"> 32: </a><font color="#B22222">@*/</font>
43
 
<a name="line33"> 33: </a><strong><font color="#4169E1"><a name="TSDefaultComputeJacobianColor"></a>int <A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>(<A href="../../../docs/manualpages/TS/TS.html#TS">TS</A> ts,PetscReal t,<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>
44
 
<a name="line34"> 34: </a>{
45
 
<a name="line35"> 35: </a>  <A href="../../../docs/manualpages/Mat/MatFDColoring.html#MatFDColoring">MatFDColoring</A> color = (<A href="../../../docs/manualpages/Mat/MatFDColoring.html#MatFDColoring">MatFDColoring</A>) ctx;
46
 
<a name="line36"> 36: </a>  <A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>          snes;
47
 
<a name="line37"> 37: </a>  int           ierr,freq,it;
48
 
 
49
 
<a name="line40"> 40: </a>  <font color="#B22222">/*</font>
50
 
<a name="line41"> 41: </a><font color="#B22222">       If we are not using <A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> we have no way to know the current iteration.</font>
51
 
<a name="line42"> 42: </a><font color="#B22222">  */</font>
52
 
<a name="line43"> 43: </a>  <A href="../../../docs/manualpages/TS/TSGetSNES.html#TSGetSNES">TSGetSNES</A>(ts,&amp;snes);
53
 
<a name="line44"> 44: </a>  <font color="#4169E1">if</font> (snes) {
54
 
<a name="line45"> 45: </a>    <A href="../../../docs/manualpages/MatFD/MatFDColoringGetFrequency.html#MatFDColoringGetFrequency">MatFDColoringGetFrequency</A>(color,&amp;freq);
55
 
<a name="line46"> 46: </a>    <A href="../../../docs/manualpages/SNES/SNESGetIterationNumber.html#SNESGetIterationNumber">SNESGetIterationNumber</A>(snes,&amp;it);
56
 
 
57
 
<a name="line48"> 48: </a>    <font color="#4169E1">if</font> ((freq &gt; 1) &amp;&amp; ((it % freq) != 1)) {
58
 
<a name="line49"> 49: </a>      <A href="../../../docs/manualpages/Profiling/PetscLogInfo.html#PetscLogInfo">PetscLogInfo</A>(color,<font color="#666666">"<A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>:Skipping Jacobian, it %d, freq %dn"</font>,it,freq);
59
 
<a name="line50"> 50: </a>      *flag = SAME_PRECONDITIONER;
60
 
<a name="line51"> 51: </a>      <font color="#4169E1">return</font>(0);
61
 
<a name="line52"> 52: </a>    } <font color="#4169E1">else</font> {
62
 
<a name="line53"> 53: </a>      <A href="../../../docs/manualpages/Profiling/PetscLogInfo.html#PetscLogInfo">PetscLogInfo</A>(color,<font color="#666666">"<A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>:Computing Jacobian, it %d, freq %dn"</font>,it,freq);
63
 
<a name="line54"> 54: </a>      *flag = SAME_NONZERO_PATTERN;
64
 
<a name="line55"> 55: </a>    }
65
 
<a name="line56"> 56: </a>  }
66
 
<a name="line57"> 57: </a>  <A href="../../../docs/manualpages/MatFD/MatFDColoringApplyTS.html#MatFDColoringApplyTS">MatFDColoringApplyTS</A>(*B,color,t,x1,flag,ts);
67
 
<a name="line58"> 58: </a>  <font color="#4169E1">return</font>(0);
68
 
<a name="line59"> 59: </a>}
69
 
 
70
 
<a name="line61"> 61: </a><font color="#B22222">/*</font>
71
 
<a name="line62"> 62: </a><font color="#B22222">   TSDefaultComputeJacobian - Computes the Jacobian using finite differences.</font>
72
 
 
73
 
<a name="line64"> 64: </a><font color="#B22222">   Input Parameters:</font>
74
 
<a name="line65"> 65: </a><font color="#B22222">+  ts - <A href="../../../docs/manualpages/TS/TS.html#TS">TS</A> context</font>
75
 
<a name="line66"> 66: </a><font color="#B22222">.  xx1 - compute Jacobian at this point</font>
76
 
<a name="line67"> 67: </a><font color="#B22222">-  ctx - application's function context, as set with <A href="../../../docs/manualpages/SNES/SNESSetFunction.html#SNESSetFunction">SNESSetFunction</A>()</font>
77
 
 
78
 
<a name="line69"> 69: </a><font color="#B22222">   Output Parameters:</font>
79
 
<a name="line70"> 70: </a><font color="#B22222">+  J - Jacobian</font>
80
 
<a name="line71"> 71: </a><font color="#B22222">.  B - preconditioner, same as Jacobian</font>
81
 
<a name="line72"> 72: </a><font color="#B22222">-  flag - matrix flag</font>
82
 
 
83
 
<a name="line74"> 74: </a><font color="#B22222">   Notes:</font>
84
 
<a name="line75"> 75: </a><font color="#B22222">   This routine is slow and expensive, and is not optimized.</font>
85
 
 
86
 
<a name="line77"> 77: </a><font color="#B22222">   Sparse approximations using colorings are also available and</font>
87
 
<a name="line78"> 78: </a><font color="#B22222">   would be a much better alternative!</font>
88
 
 
89
 
<a name="line80"> 80: </a><font color="#B22222">   Level: intermediate</font>
90
 
 
91
 
<a name="line82"> 82: </a><font color="#B22222">.seealso: <A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>()</font>
92
 
<a name="line83"> 83: </a><font color="#B22222">*/</font>
93
 
<a name="line84"> 84: </a><strong><font color="#4169E1"><a name="TSDefaultComputeJacobian"></a>int TSDefaultComputeJacobian(<A href="../../../docs/manualpages/TS/TS.html#TS">TS</A> ts,PetscReal t,<A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> xx1,<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>
94
 
<a name="line85"> 85: </a>{
95
 
<a name="line86"> 86: </a>  <A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>         jj1,jj2,xx2;
96
 
<a name="line87"> 87: </a>  int         i,ierr,N,start,end,j;
97
 
<a name="line88"> 88: </a>  <A href="../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A> dx,mone = -1.0,*y,scale,*xx,wscale;
98
 
<a name="line89"> 89: </a>  PetscReal   amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
99
 
<a name="line90"> 90: </a>  PetscReal   dx_min = 1.e-16,dx_par = 1.e-1;
100
 
<a name="line91"> 91: </a>  MPI_Comm    <A href="../../../docs/manualpages/Sys/comm.html#comm">comm</A>;
101
 
 
102
 
<a name="line94"> 94: </a>  <A href="../../../docs/manualpages/Vec/VecDuplicate.html#VecDuplicate">VecDuplicate</A>(xx1,&amp;jj1);
103
 
<a name="line95"> 95: </a>  <A href="../../../docs/manualpages/Vec/VecDuplicate.html#VecDuplicate">VecDuplicate</A>(xx1,&amp;jj2);
104
 
<a name="line96"> 96: </a>  <A href="../../../docs/manualpages/Vec/VecDuplicate.html#VecDuplicate">VecDuplicate</A>(xx1,&amp;xx2);
105
 
 
106
 
<a name="line98"> 98: </a>  <A href="../../../docs/manualpages/Sys/PetscObjectGetComm.html#PetscObjectGetComm">PetscObjectGetComm</A>((<A href="../../../docs/manualpages/Sys/PetscObject.html#PetscObject">PetscObject</A>)xx1,&amp;<A href="../../../docs/manualpages/Sys/comm.html#comm">comm</A>);
107
 
<a name="line99"> 99: </a>  <A href="../../../docs/manualpages/Mat/MatZeroEntries.html#MatZeroEntries">MatZeroEntries</A>(*J);
108
 
 
109
 
<a name="line101">101: </a>  <A href="../../../docs/manualpages/Vec/VecGetSize.html#VecGetSize">VecGetSize</A>(xx1,&amp;N);
110
 
<a name="line102">102: </a>  <A href="../../../docs/manualpages/Vec/VecGetOwnershipRange.html#VecGetOwnershipRange">VecGetOwnershipRange</A>(xx1,&amp;start,&amp;end);
111
 
<a name="line103">103: </a>  TSComputeRHSFunction(ts,ts-&gt;ptime,xx1,jj1);
112
 
 
113
 
<a name="line105">105: </a>  <font color="#B22222">/* Compute Jacobian approximation, 1 column at a time.</font>
114
 
<a name="line106">106: </a><font color="#B22222">      xx1 = current iterate, jj1 = F(xx1)</font>
115
 
<a name="line107">107: </a><font color="#B22222">      xx2 = perturbed iterate, jj2 = F(xx2)</font>
116
 
<a name="line108">108: </a><font color="#B22222">   */</font>
117
 
<a name="line109">109: </a>  <font color="#4169E1">for</font> (i=0; i&lt;N; i++) {
118
 
<a name="line110">110: </a>    <A href="../../../docs/manualpages/Vec/VecCopy.html#VecCopy">VecCopy</A>(xx1,xx2);
119
 
<a name="line111">111: </a>    <font color="#4169E1">if</font> (i&gt;= start &amp;&amp; i&lt;end) {
120
 
<a name="line112">112: </a>       <A href="../../../docs/manualpages/Vec/VecGetArray.html#VecGetArray">VecGetArray</A>(xx1,&amp;xx);
121
 
<a name="line113">113: </a>      dx   = xx[i-start];
122
 
<a name="line114">114: </a>       <A href="../../../docs/manualpages/Vec/VecRestoreArray.html#VecRestoreArray">VecRestoreArray</A>(xx1,&amp;xx);
123
 
<a name="line115">115: </a><font color="#A020F0">#if !defined(PETSC_USE_COMPLEX)</font>
124
 
<a name="line116">116: </a>      <font color="#4169E1">if</font> (dx &lt; dx_min &amp;&amp; dx &gt;= 0.0) dx = dx_par;
125
 
<a name="line117">117: </a>      <font color="#4169E1">else</font> <font color="#4169E1">if</font> (dx &lt; 0.0 &amp;&amp; dx &gt; -dx_min) dx = -dx_par;
126
 
<a name="line118">118: </a><font color="#A020F0">#else</font>
127
 
<a name="line119">119: </a>      <font color="#4169E1">if</font> (PetscAbsScalar(dx) &lt; dx_min &amp;&amp; PetscRealPart(dx) &gt;= 0.0) dx = dx_par;
128
 
<a name="line120">120: </a>      <font color="#4169E1">else</font> <font color="#4169E1">if</font> (PetscRealPart(dx) &lt; 0.0 &amp;&amp; PetscAbsScalar(dx) &lt; dx_min) dx = -dx_par;
129
 
<a name="line121">121: </a><font color="#A020F0">#endif</font>
130
 
<a name="line122">122: </a>      dx *= epsilon;
131
 
<a name="line123">123: </a>      wscale = 1.0/dx;
132
 
<a name="line124">124: </a>       <A href="../../../docs/manualpages/Vec/VecSetValues.html#VecSetValues">VecSetValues</A>(xx2,1,&amp;i,&amp;dx,ADD_VALUES);
133
 
<a name="line125">125: </a>    } <font color="#4169E1">else</font> {
134
 
<a name="line126">126: </a>      wscale = 0.0;
135
 
<a name="line127">127: </a>    }
136
 
<a name="line128">128: </a>    TSComputeRHSFunction(ts,t,xx2,jj2);
137
 
<a name="line129">129: </a>    <A href="../../../docs/manualpages/Vec/VecAXPY.html#VecAXPY">VecAXPY</A>(&amp;mone,jj1,jj2);
138
 
<a name="line130">130: </a>    <font color="#B22222">/* Communicate scale to all processors */</font>
139
 
<a name="line131">131: </a>    <A href="http://www.mcs.anl.gov/mpi/www/www3/MPI_Allreduce.html#MPI_Allreduce">MPI_Allreduce</A>(&amp;wscale,&amp;scale,1,MPIU_SCALAR,PetscSum_Op,<A href="../../../docs/manualpages/Sys/comm.html#comm">comm</A>);
140
 
<a name="line132">132: </a>    <A href="../../../docs/manualpages/Vec/VecScale.html#VecScale">VecScale</A>(&amp;scale,jj2);
141
 
<a name="line133">133: </a>    <A href="../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(jj2,NORM_INFINITY,&amp;amax); amax *= 1.e-14;
142
 
<a name="line134">134: </a>    <A href="../../../docs/manualpages/Vec/VecGetArray.html#VecGetArray">VecGetArray</A>(jj2,&amp;y);
143
 
<a name="line135">135: </a>    <font color="#4169E1">for</font> (j=start; j&lt;end; j++) {
144
 
<a name="line136">136: </a>      <font color="#4169E1">if</font> (PetscAbsScalar(y[j-start]) &gt; amax) {
145
 
<a name="line137">137: </a>        <A href="../../../docs/manualpages/Mat/MatSetValues.html#MatSetValues">MatSetValues</A>(*J,1,&amp;j,1,&amp;i,y+j-start,INSERT_VALUES);
146
 
<a name="line138">138: </a>      }
147
 
<a name="line139">139: </a>    }
148
 
<a name="line140">140: </a>    <A href="../../../docs/manualpages/Vec/VecRestoreArray.html#VecRestoreArray">VecRestoreArray</A>(jj2,&amp;y);
149
 
<a name="line141">141: </a>  }
150
 
<a name="line142">142: </a>  <A href="../../../docs/manualpages/Mat/MatAssemblyBegin.html#MatAssemblyBegin">MatAssemblyBegin</A>(*J,MAT_FINAL_ASSEMBLY);
151
 
<a name="line143">143: </a>  <A href="../../../docs/manualpages/Mat/MatAssemblyEnd.html#MatAssemblyEnd">MatAssemblyEnd</A>(*J,MAT_FINAL_ASSEMBLY);
152
 
<a name="line144">144: </a>  *flag =  DIFFERENT_NONZERO_PATTERN;
153
 
 
154
 
<a name="line146">146: </a>  <A href="../../../docs/manualpages/Vec/VecDestroy.html#VecDestroy">VecDestroy</A>(jj1);
155
 
<a name="line147">147: </a>  <A href="../../../docs/manualpages/Vec/VecDestroy.html#VecDestroy">VecDestroy</A>(jj2);
156
 
<a name="line148">148: </a>  <A href="../../../docs/manualpages/Vec/VecDestroy.html#VecDestroy">VecDestroy</A>(xx2);
157
 
 
158
 
<a name="line150">150: </a>  <font color="#4169E1">return</font>(0);
159
 
<a name="line151">151: </a>}
 
16
<a name="line8">  8: </a><font color="#B22222">/*@C</font>
 
17
<a name="line9">  9: </a><font color="#B22222">    <A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A> - Computes the Jacobian using</font>
 
18
<a name="line10"> 10: </a><font color="#B22222">    finite differences and coloring to exploit matrix sparsity.  </font>
 
19
<a name="line11"> 11: </a><font color="#B22222">  </font>
 
20
<a name="line12"> 12: </a><font color="#B22222">    Collective on <A href="../../../docs/manualpages/TS/TS.html#TS">TS</A>, <A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> and <A href="../../../docs/manualpages/Mat/Mat.html#Mat">Mat</A></font>
 
21
 
 
22
<a name="line14"> 14: </a><font color="#B22222">    Input Parameters:</font>
 
23
<a name="line15"> 15: </a><font color="#B22222">+   ts - nonlinear solver object</font>
 
24
<a name="line16"> 16: </a><font color="#B22222">.   t - current time</font>
 
25
<a name="line17"> 17: </a><font color="#B22222">.   x1 - location at which to evaluate Jacobian</font>
 
26
<a name="line18"> 18: </a><font color="#B22222">-   ctx - coloring context, where ctx must have type <A href="../../../docs/manualpages/Mat/MatFDColoring.html#MatFDColoring">MatFDColoring</A>, </font>
 
27
<a name="line19"> 19: </a><font color="#B22222">          as created via <A href="../../../docs/manualpages/MatFD/MatFDColoringCreate.html#MatFDColoringCreate">MatFDColoringCreate</A>()</font>
 
28
 
 
29
<a name="line21"> 21: </a><font color="#B22222">    Output Parameters:</font>
 
30
<a name="line22"> 22: </a><font color="#B22222">+   J - Jacobian matrix (not altered in this routine)</font>
 
31
<a name="line23"> 23: </a><font color="#B22222">.   B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)</font>
 
32
<a name="line24"> 24: </a><font color="#B22222">-   flag - flag indicating whether the matrix sparsity structure has changed</font>
 
33
 
 
34
<a name="line26"> 26: </a><font color="#B22222">   Options Database Keys:</font>
 
35
<a name="line27"> 27: </a><font color="#B22222">$  -mat_fd_coloring_freq &lt;freq&gt; </font>
 
36
 
 
37
<a name="line29"> 29: </a><font color="#B22222">   Level: intermediate</font>
 
38
 
 
39
<a name="line31"> 31: </a><font color="#B22222">.keywords: <A href="../../../docs/manualpages/TS/TS.html#TS">TS</A>, finite differences, Jacobian, coloring, sparse</font>
 
40
 
 
41
<a name="line33"> 33: </a><font color="#B22222">.seealso: TSSetJacobian(), <A href="../../../docs/manualpages/MatFD/MatFDColoringCreate.html#MatFDColoringCreate">MatFDColoringCreate</A>(), <A href="../../../docs/manualpages/MatFD/MatFDColoringSetFunction.html#MatFDColoringSetFunction">MatFDColoringSetFunction</A>()</font>
 
42
<a name="line34"> 34: </a><font color="#B22222">@*/</font>
 
43
<a name="line35"> 35: </a><strong><font color="#4169E1"><a name="TSDefaultComputeJacobianColor"></a>int <A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>(<A href="../../../docs/manualpages/TS/TS.html#TS">TS</A> ts,<A href="../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A> t,<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>
 
44
<a name="line36"> 36: </a>{
 
45
<a name="line37"> 37: </a>  <A href="../../../docs/manualpages/Mat/MatFDColoring.html#MatFDColoring">MatFDColoring</A> color = (<A href="../../../docs/manualpages/Mat/MatFDColoring.html#MatFDColoring">MatFDColoring</A>) ctx;
 
46
<a name="line38"> 38: </a>  <A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A>          snes;
 
47
<a name="line39"> 39: </a>  int           ierr,freq,it;
 
48
 
 
49
<a name="line42"> 42: </a>  <font color="#B22222">/*</font>
 
50
<a name="line43"> 43: </a><font color="#B22222">       If we are not using <A href="../../../docs/manualpages/SNES/SNES.html#SNES">SNES</A> we have no way to know the current iteration.</font>
 
51
<a name="line44"> 44: </a><font color="#B22222">  */</font>
 
52
<a name="line45"> 45: </a>  <A href="../../../docs/manualpages/TS/TSGetSNES.html#TSGetSNES">TSGetSNES</A>(ts,&amp;snes);
 
53
<a name="line46"> 46: </a>  <font color="#4169E1">if</font> (snes) {
 
54
<a name="line47"> 47: </a>    <A href="../../../docs/manualpages/MatFD/MatFDColoringGetFrequency.html#MatFDColoringGetFrequency">MatFDColoringGetFrequency</A>(color,&amp;freq);
 
55
<a name="line48"> 48: </a>    <A href="../../../docs/manualpages/SNES/SNESGetIterationNumber.html#SNESGetIterationNumber">SNESGetIterationNumber</A>(snes,&amp;it);
 
56
 
 
57
<a name="line50"> 50: </a>    <font color="#4169E1">if</font> ((freq &gt; 1) &amp;&amp; ((it % freq) != 1)) {
 
58
<a name="line51"> 51: </a>      <A href="../../../docs/manualpages/Profiling/PetscLogInfo.html#PetscLogInfo">PetscLogInfo</A>(color,<font color="#666666">"<A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>:Skipping Jacobian, it %d, freq %d\n"</font>,it,freq);
 
59
<a name="line52"> 52: </a>      *flag = SAME_PRECONDITIONER;
 
60
<a name="line53"> 53: </a>      <font color="#4169E1">return</font>(0);
 
61
<a name="line54"> 54: </a>    } <font color="#4169E1">else</font> {
 
62
<a name="line55"> 55: </a>      <A href="../../../docs/manualpages/Profiling/PetscLogInfo.html#PetscLogInfo">PetscLogInfo</A>(color,<font color="#666666">"<A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>:Computing Jacobian, it %d, freq %d\n"</font>,it,freq);
 
63
<a name="line56"> 56: </a>      *flag = SAME_NONZERO_PATTERN;
 
64
<a name="line57"> 57: </a>    }
 
65
<a name="line58"> 58: </a>  }
 
66
<a name="line59"> 59: </a>  <A href="../../../docs/manualpages/MatFD/MatFDColoringApplyTS.html#MatFDColoringApplyTS">MatFDColoringApplyTS</A>(*B,color,t,x1,flag,ts);
 
67
<a name="line60"> 60: </a>  <font color="#4169E1">return</font>(0);
 
68
<a name="line61"> 61: </a>}
 
69
 
 
70
<a name="line65"> 65: </a><font color="#B22222">/*</font>
 
71
<a name="line66"> 66: </a><font color="#B22222">   TSDefaultComputeJacobian - Computes the Jacobian using finite differences.</font>
 
72
 
 
73
<a name="line68"> 68: </a><font color="#B22222">   Input Parameters:</font>
 
74
<a name="line69"> 69: </a><font color="#B22222">+  ts - <A href="../../../docs/manualpages/TS/TS.html#TS">TS</A> context</font>
 
75
<a name="line70"> 70: </a><font color="#B22222">.  xx1 - compute Jacobian at this point</font>
 
76
<a name="line71"> 71: </a><font color="#B22222">-  ctx - application's function context, as set with <A href="../../../docs/manualpages/SNES/SNESSetFunction.html#SNESSetFunction">SNESSetFunction</A>()</font>
 
77
 
 
78
<a name="line73"> 73: </a><font color="#B22222">   Output Parameters:</font>
 
79
<a name="line74"> 74: </a><font color="#B22222">+  J - Jacobian</font>
 
80
<a name="line75"> 75: </a><font color="#B22222">.  B - preconditioner, same as Jacobian</font>
 
81
<a name="line76"> 76: </a><font color="#B22222">-  flag - matrix flag</font>
 
82
 
 
83
<a name="line78"> 78: </a><font color="#B22222">   Notes:</font>
 
84
<a name="line79"> 79: </a><font color="#B22222">   This routine is slow and expensive, and is not optimized.</font>
 
85
 
 
86
<a name="line81"> 81: </a><font color="#B22222">   Sparse approximations using colorings are also available and</font>
 
87
<a name="line82"> 82: </a><font color="#B22222">   would be a much better alternative!</font>
 
88
 
 
89
<a name="line84"> 84: </a><font color="#B22222">   Level: intermediate</font>
 
90
 
 
91
<a name="line86"> 86: </a><font color="#B22222">.seealso: <A href="../../../docs/manualpages/TS/TSDefaultComputeJacobianColor.html#TSDefaultComputeJacobianColor">TSDefaultComputeJacobianColor</A>()</font>
 
92
<a name="line87"> 87: </a><font color="#B22222">*/</font>
 
93
<a name="line88"> 88: </a><strong><font color="#4169E1"><a name="TSDefaultComputeJacobian"></a>int TSDefaultComputeJacobian(<A href="../../../docs/manualpages/TS/TS.html#TS">TS</A> ts,<A href="../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A> t,<A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A> xx1,<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>
 
94
<a name="line89"> 89: </a>{
 
95
<a name="line90"> 90: </a>  <A href="../../../docs/manualpages/Vec/Vec.html#Vec">Vec</A>         jj1,jj2,xx2;
 
96
<a name="line91"> 91: </a>  int         i,ierr,N,start,end,j;
 
97
<a name="line92"> 92: </a>  <A href="../../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A> dx,mone = -1.0,*y,scale,*xx,wscale;
 
98
<a name="line93"> 93: </a>  <A href="../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>   amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
 
99
<a name="line94"> 94: </a>  <A href="../../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>   dx_min = 1.e-16,dx_par = 1.e-1;
 
100
<a name="line95"> 95: </a>  <A href="../../../docs/manualpages/Sys/MPI_Comm.html#MPI_Comm">MPI_Comm</A>    <A href="../../../docs/manualpages/Sys/comm.html#comm">comm</A>;
 
101
 
 
102
<a name="line98"> 98: </a>  <A href="../../../docs/manualpages/Vec/VecDuplicate.html#VecDuplicate">VecDuplicate</A>(xx1,&amp;jj1);
 
103
<a name="line99"> 99: </a>  <A href="../../../docs/manualpages/Vec/VecDuplicate.html#VecDuplicate">VecDuplicate</A>(xx1,&amp;jj2);
 
104
<a name="line100">100: </a>  <A href="../../../docs/manualpages/Vec/VecDuplicate.html#VecDuplicate">VecDuplicate</A>(xx1,&amp;xx2);
 
105
 
 
106
<a name="line102">102: </a>  <A href="../../../docs/manualpages/Sys/PetscObjectGetComm.html#PetscObjectGetComm">PetscObjectGetComm</A>((<A href="../../../docs/manualpages/Sys/PetscObject.html#PetscObject">PetscObject</A>)xx1,&amp;<A href="../../../docs/manualpages/Sys/comm.html#comm">comm</A>);
 
107
<a name="line103">103: </a>  <A href="../../../docs/manualpages/Mat/MatZeroEntries.html#MatZeroEntries">MatZeroEntries</A>(*J);
 
108
 
 
109
<a name="line105">105: </a>  <A href="../../../docs/manualpages/Vec/VecGetSize.html#VecGetSize">VecGetSize</A>(xx1,&amp;N);
 
110
<a name="line106">106: </a>  <A href="../../../docs/manualpages/Vec/VecGetOwnershipRange.html#VecGetOwnershipRange">VecGetOwnershipRange</A>(xx1,&amp;start,&amp;end);
 
111
<a name="line107">107: </a>  TSComputeRHSFunction(ts,ts-&gt;ptime,xx1,jj1);
 
112
 
 
113
<a name="line109">109: </a>  <font color="#B22222">/* Compute Jacobian approximation, 1 column at a time.</font>
 
114
<a name="line110">110: </a><font color="#B22222">      xx1 = current iterate, jj1 = F(xx1)</font>
 
115
<a name="line111">111: </a><font color="#B22222">      xx2 = perturbed iterate, jj2 = F(xx2)</font>
 
116
<a name="line112">112: </a><font color="#B22222">   */</font>
 
117
<a name="line113">113: </a>  <font color="#4169E1">for</font> (i=0; i&lt;N; i++) {
 
118
<a name="line114">114: </a>    <A href="../../../docs/manualpages/Vec/VecCopy.html#VecCopy">VecCopy</A>(xx1,xx2);
 
119
<a name="line115">115: </a>    <font color="#4169E1">if</font> (i&gt;= start &amp;&amp; i&lt;end) {
 
120
<a name="line116">116: </a>       <A href="../../../docs/manualpages/Vec/VecGetArray.html#VecGetArray">VecGetArray</A>(xx1,&amp;xx);
 
121
<a name="line117">117: </a>      dx   = xx[i-start];
 
122
<a name="line118">118: </a>       <A href="../../../docs/manualpages/Vec/VecRestoreArray.html#VecRestoreArray">VecRestoreArray</A>(xx1,&amp;xx);
 
123
<a name="line119">119: </a><font color="#A020F0">#if !defined(PETSC_USE_COMPLEX)</font>
 
124
<a name="line120">120: </a>      <font color="#4169E1">if</font> (dx &lt; dx_min &amp;&amp; dx &gt;= 0.0) dx = dx_par;
 
125
<a name="line121">121: </a>      <font color="#4169E1">else</font> <font color="#4169E1">if</font> (dx &lt; 0.0 &amp;&amp; dx &gt; -dx_min) dx = -dx_par;
 
126
<a name="line122">122: </a><font color="#A020F0">#else</font>
 
127
<a name="line123">123: </a>      <font color="#4169E1">if</font> (PetscAbsScalar(dx) &lt; dx_min &amp;&amp; PetscRealPart(dx) &gt;= 0.0) dx = dx_par;
 
128
<a name="line124">124: </a>      <font color="#4169E1">else</font> <font color="#4169E1">if</font> (PetscRealPart(dx) &lt; 0.0 &amp;&amp; PetscAbsScalar(dx) &lt; dx_min) dx = -dx_par;
 
129
<a name="line125">125: </a><font color="#A020F0">#endif</font>
 
130
<a name="line126">126: </a>      dx *= epsilon;
 
131
<a name="line127">127: </a>      wscale = 1.0/dx;
 
132
<a name="line128">128: </a>       <A href="../../../docs/manualpages/Vec/VecSetValues.html#VecSetValues">VecSetValues</A>(xx2,1,&amp;i,&amp;dx,ADD_VALUES);
 
133
<a name="line129">129: </a>    } <font color="#4169E1">else</font> {
 
134
<a name="line130">130: </a>      wscale = 0.0;
 
135
<a name="line131">131: </a>    }
 
136
<a name="line132">132: </a>    TSComputeRHSFunction(ts,t,xx2,jj2);
 
137
<a name="line133">133: </a>    <A href="../../../docs/manualpages/Vec/VecAXPY.html#VecAXPY">VecAXPY</A>(&amp;mone,jj1,jj2);
 
138
<a name="line134">134: </a>    <font color="#B22222">/* Communicate scale to all processors */</font>
 
139
<a name="line135">135: </a>    <A href="http://www.mcs.anl.gov/mpi/www/www3/MPI_Allreduce.html#MPI_Allreduce">MPI_Allreduce</A>(&amp;wscale,&amp;scale,1,MPIU_SCALAR,PetscSum_Op,<A href="../../../docs/manualpages/Sys/comm.html#comm">comm</A>);
 
140
<a name="line136">136: </a>    <A href="../../../docs/manualpages/Vec/VecScale.html#VecScale">VecScale</A>(&amp;scale,jj2);
 
141
<a name="line137">137: </a>    <A href="../../../docs/manualpages/Vec/VecNorm.html#VecNorm">VecNorm</A>(jj2,NORM_INFINITY,&amp;amax); amax *= 1.e-14;
 
142
<a name="line138">138: </a>    <A href="../../../docs/manualpages/Vec/VecGetArray.html#VecGetArray">VecGetArray</A>(jj2,&amp;y);
 
143
<a name="line139">139: </a>    <font color="#4169E1">for</font> (j=start; j&lt;end; j++) {
 
144
<a name="line140">140: </a>      <font color="#4169E1">if</font> (PetscAbsScalar(y[j-start]) &gt; amax) {
 
145
<a name="line141">141: </a>        <A href="../../../docs/manualpages/Mat/MatSetValues.html#MatSetValues">MatSetValues</A>(*J,1,&amp;j,1,&amp;i,y+j-start,INSERT_VALUES);
 
146
<a name="line142">142: </a>      }
 
147
<a name="line143">143: </a>    }
 
148
<a name="line144">144: </a>    <A href="../../../docs/manualpages/Vec/VecRestoreArray.html#VecRestoreArray">VecRestoreArray</A>(jj2,&amp;y);
 
149
<a name="line145">145: </a>  }
 
150
<a name="line146">146: </a>  <A href="../../../docs/manualpages/Mat/MatAssemblyBegin.html#MatAssemblyBegin">MatAssemblyBegin</A>(*J,MAT_FINAL_ASSEMBLY);
 
151
<a name="line147">147: </a>  <A href="../../../docs/manualpages/Mat/MatAssemblyEnd.html#MatAssemblyEnd">MatAssemblyEnd</A>(*J,MAT_FINAL_ASSEMBLY);
 
152
<a name="line148">148: </a>  *flag =  DIFFERENT_NONZERO_PATTERN;
 
153
 
 
154
<a name="line150">150: </a>  <A href="../../../docs/manualpages/Vec/VecDestroy.html#VecDestroy">VecDestroy</A>(jj1);
 
155
<a name="line151">151: </a>  <A href="../../../docs/manualpages/Vec/VecDestroy.html#VecDestroy">VecDestroy</A>(jj2);
 
156
<a name="line152">152: </a>  <A href="../../../docs/manualpages/Vec/VecDestroy.html#VecDestroy">VecDestroy</A>(xx2);
 
157
 
 
158
<a name="line154">154: </a>  <font color="#4169E1">return</font>(0);
 
159
<a name="line155">155: </a>}
160
160
 
161
161
 
162
162
</pre>