autodiff.cpp 33.1 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include "inmost.h"

#if defined(USE_AUTODIFF)

#if defined(USE_AUTODIFF_ASMJIT)
#include "asmjit.h"
#endif

#if defined(USE_AUTODIFF_OPENCL)
#include <CL/cl.h>
#endif


namespace INMOST
{
Kirill Terekhov's avatar
Kirill Terekhov committed
16
17
18
19
20
21
#if defined(NEW_VERSION)
	
	//! returns offset from the end of precomputed z
	void Automatizator::DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Solver::Row & entries, INMOST_DATA_REAL_TYPE multval, void * user_data)
	{
		INMOST_DATA_ENUM_TYPE voffset = var.values_offset(element), doffset = var.derivatives_offset(element);
Kirill Terekhov's avatar
Kirill Terekhov committed
22
		INMOST_DATA_ENUM_TYPE k = static_cast<INMOST_DATA_ENUM_TYPE>(var.data.size()-1);
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
23
		Storage e = Storage(m,var.current_stencil[element].first);
Kirill Terekhov's avatar
Kirill Terekhov committed
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
		INMOST_DATA_REAL_TYPE lval, rval, ret;
		var.values[doffset+k] = multval;
		expr::expr_data * arr = &var.data[0], *it;
		//for (expr::data_type::reverse_iterator it = var.data.rbegin(); it != var.data.rend(); ++it)
		do
		{
			it = arr+k;
			switch (it->op)
			{
			case AD_EXT: 
				assert(parent != ENUMUNDEF); 
				it->left.e->values[it->left.e->derivatives_offset(parent)+it->right.i] += var.values[doffset+k]; 
				break;
			case AD_COND:
				{
					expr & next = var.values[voffset+it->left.i] > 0.0 ? *it->right.q->left.e : *it->right.q->right.e;
					DerivativeFill(next, 0, element, entries,var.values[doffset+k], user_data);
				}
				break;
			case AD_PLUS:  
				var.values[doffset+it->left.i] += var.values[doffset+k];
				var.values[doffset+it->right.i] += var.values[doffset+k];
				break;
			case AD_MINUS:
				var.values[doffset+it->left.i] += var.values[doffset+k];
				var.values[doffset+it->right.i] -= var.values[doffset+k];
				break;
			case AD_MULT:
				rval = var.values[voffset+it->right.i];
				lval = var.values[voffset+it->left.i];
				var.values[doffset+it->left.i] += var.values[doffset+k]*rval;
				var.values[doffset+it->right.i] += var.values[doffset+k]*lval;
				break;
			case AD_DIV:
				rval = var.values[voffset+it->right.i];
				lval = var.values[voffset+it->left.i];
				var.values[doffset+it->left.i] += var.values[doffset+k]/rval;
				var.values[doffset+it->right.i] -= var.values[doffset+k]*lval/(rval*rval);
				break;
			case AD_POW:
				ret = var.values[voffset+k];
				rval = var.values[voffset+it->right.i];
				lval = var.values[voffset+it->left.i];
				var.values[doffset+it->right.i] += var.values[doffset+k] * ::log(lval);
				var.values[doffset+it->left.i] += var.values[doffset+k] * ret * rval / lval;
				break;
			case AD_SQRT:
				ret = var.values[voffset+k];
				var.values[doffset+it->left.i] += 0.5 * var.values[doffset+k] / ret;
				break;
			case AD_ABS:
				lval = var.values[voffset+it->left.i];
				var.values[doffset+it->left.i] += var.values[doffset+k] * (lval > 0.0 ? 1.0 : -1.0);
				break;
			case AD_EXP:
				ret = var.values[voffset+k];
				var.values[doffset+it->left.i] += var.values[doffset+k] * ret;
				break;
			case AD_LOG:   
				lval = var.values[voffset+it->left.i];
				var.values[doffset+it->left.i] += var.values[doffset+k] / lval;
				break;
			case AD_SIN:   
				lval = var.values[voffset+it->left.i];
				var.values[doffset+it->left.i] += var.values[doffset+k] * ::cos(lval);
				break;
			case AD_COS:   
				lval = var.values[voffset+it->left.i];
				var.values[doffset+it->left.i] -= var.values[doffset+k] * ::sin(lval);
				break;
			case AD_COND_MARK: break;
			case AD_COND_TYPE: break;
			case AD_CONST: break;
			case AD_MES: break;
			case AD_VAL: 
				{
					expr & next = *it->right.e;
					next.current_stencil.resize(1);
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
102
					next.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
103
104
105
106
107
108
109
110
111
					next.resize_for_stencil();
					var.values[doffset+it->left.i] += var.values[doffset+k] * EvaluateSub(next,0,element,user_data); 
					break;
				}
			default:
				if (it->op >= AD_FUNC) {}
				else if (it->op >= AD_TABLE) 
				{
					lval = var.values[voffset+it->left.i];
Kirill Terekhov's avatar
Kirill Terekhov committed
112
					var.values[doffset+it->left.i] += var.values[doffset+k] *  reg_tables[it->op-AD_TABLE]->get_derivative(lval);
Kirill Terekhov's avatar
Kirill Terekhov committed
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
				}
				else if (it->op >= AD_STNCL)
				{
					for (INMOST_DATA_ENUM_TYPE j = 0; j < it->left.e->current_stencil.size(); ++j) if(  it->left.e->current_stencil[j].first != NULL )
					{
						DerivativeFill(*it->left.e,  j, ENUMUNDEF,  entries, var.values[doffset+k] *  it->left.e->current_stencil[j].second, user_data);
					}
				}
				else if (it->op >= AD_CTAG) {}
				else if (it->op >= AD_TAG)
				{
					if (isDynamicValid(e, it->op))
					{
						entries[GetDynamicIndex(e,it->op,it->left.i)] += var.values[doffset+k];
					}
				}
				else assert(false);
			}
		} while(k-- != 0);
	}
	INMOST_DATA_REAL_TYPE Automatizator::EvaluateSub(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, void * user_data)
	{
		INMOST_DATA_ENUM_TYPE k = 0, offset = var.values_offset(element);
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
136
		Storage e = Storage(m,var.current_stencil[element].first);
Kirill Terekhov's avatar
Kirill Terekhov committed
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
		expr::expr_data * arr = &var.data[0], *it;
		//for (expr::data_type::iterator it = var.data.begin(); it != var.data.end(); ++it)
		do
		{
			it = arr+k;
			switch (it->op)
			{
			case AD_EXT: 
				assert(parent != ENUMUNDEF); 
				var.values[offset+k] = it->left.e->values[it->left.e->values_offset(parent)+it->right.i]; 
				break;
			case AD_COND:
				{
					expr & next = var.values[offset+it->left.i] > 0.0 ? *it->right.q->left.e : *it->right.q->right.e;
					next.current_stencil.resize(1);
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
152
					next.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
					next.resize_for_stencil();
					var.values[offset+k] = EvaluateSub(next, 0,element, user_data);
				}
				break;
			case AD_PLUS:  var.values[offset+k] = var.values[offset+it->left.i] + var.values[offset+it->right.i]; break;
			case AD_MINUS: var.values[offset+k] = var.values[offset+it->left.i] - var.values[offset+it->right.i]; break;
			case AD_MULT:  var.values[offset+k] = var.values[offset+it->left.i] * var.values[offset+it->right.i]; break;
			case AD_DIV:   var.values[offset+k] = var.values[offset+it->left.i] / var.values[offset+it->right.i]; break;
			case AD_POW:   var.values[offset+k] = ::pow(var.values[offset+it->left.i], var.values[offset+it->right.i]); break;
			case AD_SQRT:  var.values[offset+k] = ::sqrt(var.values[offset+it->left.i]); break;
			case AD_ABS:   var.values[offset+k] = ::fabs(var.values[offset+it->left.i]); break;
			case AD_EXP:   var.values[offset+k] = ::exp(var.values[offset+it->left.i]); break;
			case AD_LOG:   var.values[offset+k] = ::log(var.values[offset+it->left.i]); break;
			case AD_SIN:   var.values[offset+k] = ::sin(var.values[offset+it->left.i]); break;
			case AD_COS:   var.values[offset+k] = ::cos(var.values[offset+it->left.i]); break;
			case AD_CONST: var.values[offset+k] = it->left.r; break;
			case AD_COND_TYPE: var.values[offset+k] = ((e->GetElementType() & it->left.i)? 1.0 : -1.0); break;
			case AD_COND_MARK: var.values[offset+k] = e->GetMarker(it->left.i) ? 1.0 : -1.0; break;
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
171
			case AD_MES: assert(!(e->GetElementType() & (ESET | MESH))); m->GetGeometricData(e->GetHandle(), MEASURE, &var.values[offset+k]); break;
Kirill Terekhov's avatar
Kirill Terekhov committed
172
173
			case AD_VAL: var.values[offset+k] = var.values[offset+it->left.i]; break;
			default:
Kirill Terekhov's avatar
Kirill Terekhov committed
174
175
				if (it->op >= AD_FUNC) var.values[offset+k] = reg_funcs[it->op-AD_FUNC].func(e, user_data);
				else if (it->op >= AD_TABLE) var.values[offset+k] = reg_tables[it->op-AD_TABLE]->get_value(var.values[offset+it->left.i]);
Kirill Terekhov's avatar
Kirill Terekhov committed
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
				else if (it->op >= AD_STNCL)
				{
					it->left.e->current_stencil.clear();
					GetStencil(it->op,e,user_data,it->left.e->current_stencil);
					it->left.e->resize_for_stencil();
					var.values[offset+k] = 0.0;
					for (INMOST_DATA_ENUM_TYPE j = 0; j < it->left.e->current_stencil.size(); ++j) if( it->left.e->current_stencil[j].first != NULL )
						var.values[offset+k] += EvaluateSub(*it->left.e,j,ENUMUNDEF,user_data) * it->left.e->current_stencil[j].second;
				}
				else if (it->op >= AD_CTAG) var.values[offset+k] = GetStaticValue(e, it->op, it->left.i);
				else if (it->op >= AD_TAG) var.values[offset+k] = GetDynamicValue(e, it->op, it->left.i);
				else assert(false);
			}
			//k++;
		} while(++k != var.data.size());
		return var.values[offset+var.data.size()-1];
	}

Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
194
	INMOST_DATA_REAL_TYPE Automatizator::Evaluate(expr & var, const Storage & e, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
195
196
	{
		var.current_stencil.resize(1);
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
197
		var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
198
199
200
201
		var.resize_for_stencil();
		return EvaluateSub(var,0,ENUMUNDEF,user_data);
	}
	
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
202
	INMOST_DATA_REAL_TYPE Automatizator::Derivative(expr & var, const Storage & e, Solver::Row & out, Storage::real multiply, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
203
204
205
	{
		INMOST_DATA_REAL_TYPE ret;
		var.current_stencil.resize(1);
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
206
		var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
207
208
209
210
211
212
		var.resize_for_stencil();
		ret = EvaluateSub(var,0,ENUMUNDEF,user_data);
		DerivativeFill(var, 0, ENUMUNDEF, out, multiply, user_data);
		return ret*multiply;
	}
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
213

Kirill Terekhov's avatar
Kirill Terekhov committed
214
	INMOST_DATA_REAL_TYPE Automatizator::DerivativePrecompute(const expr & var, const Storage & e, precomp_values_t & values, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
215
216
217
218
219
	{
		assert(var.op != AD_NONE);
		INMOST_DATA_REAL_TYPE lval, rval, ret = 0.0;
		switch (var.op)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
220
		case AD_COND_TYPE:
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
221
			lval = (e->GetElementType() & (*(ElementType*)&var.left))? 1.0 : -1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
222
223
			return lval;
		case AD_COND_MARK:
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
224
			lval = e->GetMarker(*(MarkerType *)&var.left)? 1.0 : -1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
225
			return lval;
Kirill Terekhov's avatar
Kirill Terekhov committed
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
		case AD_COND:
			lval = Evaluate(*var.left, e, user_data);
			rval = DerivativePrecompute(*(lval > 0.0 ? var.right->left : var.right->right), e, values, user_data);
			values.push_back(lval);
			return rval*var.coef;
		case AD_PLUS:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			rval = DerivativePrecompute(*var.right, e, values, user_data);
			return (lval + rval)*var.coef;
		case AD_MINUS:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			rval = DerivativePrecompute(*var.right, e, values, user_data);
			return (lval - rval)*var.coef;
		case AD_MULT:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			rval = DerivativePrecompute(*var.right, e, values, user_data);
			values.push_back(lval);
			values.push_back(rval);
			return (lval * rval)*var.coef;
		case AD_DIV:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			rval = DerivativePrecompute(*var.right, e, values, user_data);
			values.push_back(lval);
			values.push_back(rval);
			return (lval / rval)*var.coef;
		case AD_POW:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			rval = DerivativePrecompute(*var.right, e, values, user_data);
			values.push_back(lval);
			values.push_back(rval);
			ret = ::pow(lval, rval);
			values.push_back(ret);
			return ret*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
259
260
261
262
263
		case AD_SQRT:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			ret = ::sqrt(lval);
			values.push_back(ret);
			return ret*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
		case AD_INV:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			values.push_back(lval);
			return var.coef / lval;
		case AD_ABS:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			values.push_back(lval);
			return ::fabs(lval)*var.coef;
		case AD_EXP:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			ret = ::exp(lval);
			values.push_back(ret);
			return ret*var.coef;
		case AD_LOG:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			values.push_back(lval);
			return ::log(lval)*var.coef;
		case AD_SIN:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			values.push_back(lval);
			return ::sin(lval)*var.coef;
		case AD_COS:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			values.push_back(lval);
			return ::cos(lval)*var.coef;
		case AD_CONST:
			return var.coef;
		case AD_MES:
			assert(!(e->GetElementType() & (ESET | MESH)));
Kirill Terekhov's avatar
Kirill Terekhov committed
293
			m->GetGeometricData(e->GetHandle(), MEASURE, &ret);
Kirill Terekhov's avatar
Kirill Terekhov committed
294
			return ret*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
295
296
297
298
299
		case AD_VAL:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			rval = Evaluate(*var.left,e,user_data);
			values.push_back(rval);
			return lval*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
300
301
302
		}
		if (var.op >= AD_FUNC)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
303
			ret = reg_funcs[var.op-AD_FUNC].func(e, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
304
305
306
307
308
309
			return ret*var.coef;
		}
		else if (var.op >= AD_TABLE)
		{
			lval = DerivativePrecompute(*var.left, e, values, user_data);
			values.push_back(lval);
Kirill Terekhov's avatar
Kirill Terekhov committed
310
			ret = reg_tables[var.op-AD_TABLE]->get_value(lval);
Kirill Terekhov's avatar
Kirill Terekhov committed
311
312
313
314
			return ret*var.coef;
		}
		else if (var.op >= AD_STNCL)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
315
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
316
317
318
319
320
321
			assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
			if (st.kind == 0)
			{
				Storage::reference_array elems = e->ReferenceArray(static_cast<stencil_tag *>(st.link)->elements);
				Storage::real_array coefs = e->RealArray(static_cast<stencil_tag *>(st.link)->coefs);
				assert(elems.size() == coefs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
322
				for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
323
324
325
326
327
328
329
330
331
				{
					lval = DerivativePrecompute(*var.left, elems[k], values, user_data);
					ret += lval * coefs[k];
				}
			}
			else if (st.kind == 1)
			{
				stencil_pairs get_st;
				reinterpret_cast<stencil_callback>(st.link)(e, get_st, user_data);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
332
				for (INMOST_DATA_ENUM_TYPE k = 0; k < get_st.size(); ++k) if( get_st[k].first != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
333
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
334
					lval = DerivativePrecompute(*var.left, Storage(m,get_st[k].first), values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
					ret += lval * get_st[k].second;
				}
			}
			return ret*var.coef;
		}
		else if (var.op >= AD_CTAG)
		{
			INMOST_DATA_ENUM_TYPE comp = *(INMOST_DATA_ENUM_TYPE *)(&var.left);
			ret = GetStaticValue(e, var.op, comp);
			return ret * var.coef;
		}
		else if (var.op >= AD_TAG)
		{
			INMOST_DATA_ENUM_TYPE comp = *(INMOST_DATA_ENUM_TYPE *)(&var.left);
			ret = GetDynamicValue(e, var.op, comp);
			return ret*var.coef;
		}
		assert(false);
		return 0.0;
	}
	//! returns offset from the end of precomputed values
Kirill Terekhov's avatar
Kirill Terekhov committed
356
	void Automatizator::DerivativeFill(const expr & var, const Storage & e, Solver::Row & entries, precomp_values_t & values, INMOST_DATA_REAL_TYPE multval, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
357
358
359
360
361
	{
		assert(var.op != AD_NONE);
		INMOST_DATA_REAL_TYPE lval, rval, ret;
		switch (var.op)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
362
363
364
		case AD_COND_MARK:
		case AD_COND_TYPE:
			return;
Kirill Terekhov's avatar
Kirill Terekhov committed
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
		case AD_COND:
			lval = values.back(); values.pop_back();
			DerivativeFill(*(lval > 0.0 ? var.right->left : var.right->right), e, entries, values, multval*var.coef, user_data);
			return;
		case AD_PLUS:
			DerivativeFill(*var.right, e, entries, values, multval*var.coef, user_data);
			DerivativeFill(*var.left, e, entries, values, multval*var.coef, user_data);
			return;
		case AD_MINUS:
			DerivativeFill(*var.right, e, entries, values, -multval*var.coef, user_data);
			DerivativeFill(*var.left, e, entries, values, multval*var.coef, user_data);
			return;
		case AD_MULT:
			rval = values.back(); values.pop_back();
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.right, e, entries, values, lval*multval*var.coef, user_data);
			DerivativeFill(*var.left, e, entries, values, rval*multval*var.coef, user_data);
			return;
		case AD_DIV:
			rval = values.back(); values.pop_back();
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.right, e, entries, values, -multval * lval / (rval*rval) * var.coef, user_data);
			DerivativeFill(*var.left, e, entries, values, multval / rval*var.coef, user_data);
			return;
		case AD_POW:
			ret = values.back(); values.pop_back();
			rval = values.back(); values.pop_back();
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.right, e, entries, values, multval * ret * ::log(lval) * var.coef, user_data);
			DerivativeFill(*var.left, e, entries, values, multval * ret * rval / lval * var.coef, user_data);
			return;
Kirill Terekhov's avatar
Kirill Terekhov committed
396
397
398
399
		case AD_SQRT:
			ret = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, 0.5 * multval / ret * var.coef, user_data);
			return;
Kirill Terekhov's avatar
Kirill Terekhov committed
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
		case AD_INV:
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, -multval / (lval * lval) * var.coef, user_data);
			return;
		case AD_ABS:
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, multval * var.coef * (lval > 0 ? 1.0 : -1.0), user_data);
			return;
		case AD_EXP:
			ret = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, multval * ret * var.coef, user_data);
			return;
		case AD_LOG:
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, multval * var.coef / lval, user_data);
			return;
		case AD_SIN:
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, multval * var.coef * ::cos(lval), user_data);
			return;
		case AD_COS:
			lval = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, -multval * var.coef * ::sin(lval), user_data);
			return;
		case AD_CONST: return;
		case AD_MES: return;
Kirill Terekhov's avatar
Kirill Terekhov committed
426
427
428
429
		case AD_VAL: 
			rval = values.back(); values.pop_back();
			DerivativeFill(*var.left, e, entries, values, multval * var.coef * rval, user_data);
			return;
Kirill Terekhov's avatar
Kirill Terekhov committed
430
431
432
433
434
435
436
437
		}
		if (var.op >= AD_FUNC)
		{
			return;
		}
		else if (var.op >= AD_TABLE)
		{
			lval = values.back(); values.pop_back();
Kirill Terekhov's avatar
Kirill Terekhov committed
438
			DerivativeFill(*var.left, e, entries, values, multval * var.coef * reg_tables[var.op-AD_TABLE]->get_derivative(lval), user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
439
440
441
442
			return;
		}
		else if (var.op >= AD_STNCL)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
443
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
444
445
446
447
448
449
			assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
			if (st.kind == 0)
			{
				Storage::reference_array elems = e->ReferenceArray(static_cast<stencil_tag *>(st.link)->elements);
				Storage::real_array coefs = e->RealArray(static_cast<stencil_tag *>(st.link)->coefs);
				assert(elems.size() == coefs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
450
				for (INMOST_DATA_ENUM_TYPE k = elems.size(); k > 0; --k) if( elems.at(k-1) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
451
452
453
454
455
456
					DerivativeFill(*var.left, elems[k - 1], entries, values, var.coef * coefs[k - 1] * multval, user_data);
			}
			else if (st.kind == 1)
			{
				stencil_pairs get_st;
				reinterpret_cast<stencil_callback>(st.link)(e, get_st, user_data);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
457
				for (INMOST_DATA_ENUM_TYPE k = static_cast<INMOST_DATA_ENUM_TYPE>(get_st.size()); k > 0; --k) if( get_st[k-1].first != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
458
					DerivativeFill(*var.left, Storage(m,get_st[k - 1].first), entries, values, var.coef * get_st[k - 1].second*multval, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
			}
			return;
		}
		else if (var.op >= AD_CTAG) return;
		else if (var.op >= AD_TAG)
		{
			if (isDynamicValid(e, var.op))
			{
				INMOST_DATA_ENUM_TYPE ind = GetDynamicIndex(e, var.op, *(INMOST_DATA_ENUM_TYPE *)(&var.left));
				entries[ind] += multval * var.coef;
			}
			return;
		}
		assert(false);
		return;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
475
	INMOST_DATA_REAL_TYPE Automatizator::Evaluate(const expr & var, const Storage & e, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
476
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
477
478
		assert(var.op != AD_NONE);
		switch (var.op)
Kirill Terekhov's avatar
Kirill Terekhov committed
479
		{
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
480
481
		case AD_COND_MARK: return e->GetMarker(*(MarkerType *)&var.left) ? 1.0 : -1.0;
		case AD_COND_TYPE: return (e->GetElementType() & (*(ElementType *)&var.left)) ? 1.0 : -1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
482
483
484
485
486
487
488
489
490
491
492
493
494
495
		case AD_COND:  return Evaluate(*(Evaluate(*var.left, e, user_data) > 0.0 ? var.right->left : var.right->right), e, user_data)*var.coef;
		case AD_PLUS:  return (Evaluate(*var.left, e, user_data) + Evaluate(*var.right, e, user_data))*var.coef;
		case AD_MINUS: return (Evaluate(*var.left, e, user_data) - Evaluate(*var.right, e, user_data))*var.coef;
		case AD_MULT:  return (Evaluate(*var.left, e, user_data) * Evaluate(*var.right, e, user_data))*var.coef;
		case AD_DIV:   return (Evaluate(*var.left, e, user_data) / Evaluate(*var.right, e, user_data))*var.coef;
		case AD_INV:   return var.coef / Evaluate(*var.left, e, user_data);
		case AD_POW:   return ::pow(Evaluate(*var.left, e, user_data), Evaluate(*var.right, e, user_data))*var.coef;
		case AD_SQRT:  return ::sqrt(Evaluate(*var.left, e, user_data))*var.coef;
		case AD_ABS:   return ::fabs(Evaluate(*var.left, e, user_data))*var.coef;
		case AD_EXP:   return ::exp(Evaluate(*var.left, e, user_data))*var.coef;
		case AD_LOG:   return ::log(Evaluate(*var.left, e, user_data))*var.coef;
		case AD_SIN:   return ::sin(Evaluate(*var.left, e, user_data))*var.coef;
		case AD_COS:   return ::cos(Evaluate(*var.left, e, user_data))*var.coef;
		case AD_CONST: return var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
496
		case AD_MES: assert(!(e->GetElementType() & (ESET | MESH))); Storage::real ret; m->GetGeometricData(e->GetHandle(), MEASURE, &ret); return ret*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
497
498
		case AD_VAL:   return Evaluate(*var.left,e,user_data)*var.coef;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
499
500
		if (var.op >= AD_FUNC) return reg_funcs[var.op-AD_FUNC].func(e, user_data);
		if (var.op >= AD_TABLE) return reg_tables[var.op-AD_TABLE]->get_value(Evaluate(*var.left, e, user_data))*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
501
502
503
		if (var.op >= AD_STNCL)
		{
			INMOST_DATA_REAL_TYPE ret = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
504
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
505
506
			assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
			if (st.kind == 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
507
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
508
509
510
				Storage::reference_array elems = e->ReferenceArray(static_cast<stencil_tag *>(st.link)->elements);
				Storage::real_array coefs = e->RealArray(static_cast<stencil_tag *>(st.link)->coefs);
				assert(elems.size() == coefs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
511
				for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
512
					ret += var.coef * Evaluate(*var.left, elems[k], user_data) * coefs[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
513
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
514
515
516
517
			else if (st.kind == 1)
			{
				stencil_pairs get_st;
				reinterpret_cast<stencil_callback>(st.link)(e, get_st, user_data);
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
518
				for (INMOST_DATA_ENUM_TYPE k = 0; k < get_st.size(); ++k) if ( get_st[k].first != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
519
					ret += var.coef * Evaluate(*var.left, Storage(m,get_st[k].first), user_data) * get_st[k].second;
Kirill Terekhov's avatar
Kirill Terekhov committed
520
521
			}
			return ret;
Kirill Terekhov's avatar
Kirill Terekhov committed
522
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
523
524
525
526
		if (var.op >= AD_CTAG) return GetStaticValue(e, var.op, *(INMOST_DATA_ENUM_TYPE *)(&var.left))*var.coef;
		if (var.op >= AD_TAG) return  GetDynamicValue(e, var.op, *(INMOST_DATA_ENUM_TYPE *)(&var.left))*var.coef;
		assert(false);
		return 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
527
528
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
529
	INMOST_DATA_REAL_TYPE Automatizator::Derivative(const expr & var, const Storage & e, Solver::Row & out, Storage::real multiply, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
530
531
532
533
	{
		INMOST_DATA_REAL_TYPE ret;
		precomp_values_t values;
		ret = DerivativePrecompute(var, e, values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
534
535
		DerivativeFill(var, e, out, values, multiply, user_data);
		return ret*multiply;
Kirill Terekhov's avatar
Kirill Terekhov committed
536
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
537
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
538
539
540
541
542
543
544
	Automatizator::Automatizator(Mesh * m) :first_num(0), last_num(0), m(m) {}
	Automatizator::~Automatizator()
	{
		for (unsigned k = 0; k < index_tags.size(); k++)
			index_tags[k].indices = m->DeleteTag(index_tags[k].indices);
		for (table_type::iterator it = reg_tables.begin(); it != reg_tables.end(); ++it)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
545
546
547
			delete[] (*it)->args;
			delete[] (*it)->vals;
			delete (*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
548
549
		}
		for (stencil_type::iterator it = reg_stencils.begin(); it != reg_stencils.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
550
551
		if (it->kind == 0)
			delete static_cast<stencil_tag *>(it->link);
Kirill Terekhov's avatar
Kirill Terekhov committed
552
553
554
	}
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterFunc(std::string name, func_callback func)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
555
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_funcs.size()) + AD_FUNC;
Kirill Terekhov's avatar
Kirill Terekhov committed
556
557
558
		func_name_callback v;
		v.name = name;
		v.func = func;
Kirill Terekhov's avatar
Kirill Terekhov committed
559
		reg_funcs.push_back(v);
Kirill Terekhov's avatar
Kirill Terekhov committed
560
561
562
		return ret;
	}
	//register stencil that can be got from tags
Kirill Terekhov's avatar
Kirill Terekhov committed
563
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterStencil(std::string name, Tag elements_tag, Tag coefs_tag, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
564
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
565
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
Kirill Terekhov's avatar
Kirill Terekhov committed
566
567
568
569
570
571
572
573
		stencil_kind_domain st;
		stencil_tag * save = new stencil_tag;
		st.name = name;
		save->coefs = coefs_tag;
		save->elements = elements_tag;
		st.kind = 0;
		st.link = static_cast<void *>(save);
		st.domainmask = domain_mask;
Kirill Terekhov's avatar
Kirill Terekhov committed
574
		reg_stencils.push_back(st);
Kirill Terekhov's avatar
Kirill Terekhov committed
575
576
577
		return ret;
	}
	//register stencil that can be got from function
Kirill Terekhov's avatar
Kirill Terekhov committed
578
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterStencil(std::string name, stencil_callback func, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
579
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
580
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
Kirill Terekhov's avatar
Kirill Terekhov committed
581
582
583
584
585
		stencil_kind_domain st;
		st.name = name;
		st.kind = 1;
		st.link = reinterpret_cast<void *>(func);
		st.domainmask = domain_mask;
Kirill Terekhov's avatar
Kirill Terekhov committed
586
		reg_stencils.push_back(st);
Kirill Terekhov's avatar
Kirill Terekhov committed
587
588
589
590
		return ret;
	}
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterTable(std::string name, INMOST_DATA_REAL_TYPE * Arguments, INMOST_DATA_REAL_TYPE * Values, INMOST_DATA_ENUM_TYPE size)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
591
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_tables.size()) + AD_TABLE;
Kirill Terekhov's avatar
Kirill Terekhov committed
592
593
594
595
596
597
598
		table_ptr t = new table;
		t->name = name;
		t->args = new INMOST_DATA_REAL_TYPE[size];
		memcpy(t->args, Arguments, sizeof(INMOST_DATA_REAL_TYPE)*size);
		t->vals = new INMOST_DATA_REAL_TYPE[size];
		memcpy(t->vals, Values, sizeof(INMOST_DATA_REAL_TYPE)*size);
		t->size = size;
Kirill Terekhov's avatar
Kirill Terekhov committed
599
		reg_tables.push_back(t);
Kirill Terekhov's avatar
Kirill Terekhov committed
600
601
602
603
		return ret;
	}
	/// set data of tag t defined on domain_mask to be dynamic data
	/// don't register tag twice
Kirill Terekhov's avatar
Kirill Terekhov committed
604
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterDynamicTag(Tag t, ElementType typemask, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
605
606
607
608
609
610
611
612
613
614
615
	{
		tagpair p;
		p.d.domain_mask = domain_mask;
		p.d.t = t;
		ElementType def = NONE, sparse = NONE;
		for (ElementType q = NODE; q <= MESH; q = q << 1) if (q & typemask)
		{
			if (t.isDefined(q)) def |= q;
			if (t.isSparse(q)) sparse |= q;
		}
		p.indices = m->CreateTag(t.GetTagName() + "_index", DATA_INTEGER, def, sparse, t.GetSize());
Kirill Terekhov's avatar
Kirill Terekhov committed
616
617
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_tags.size()) + AD_TAG;
		reg_tags.push_back(p);
Kirill Terekhov's avatar
Kirill Terekhov committed
618
619
620
621
622
623
624
625
626
627
628
		index_tags.push_back(p);
		return ret;
	}
	/// set index for every data entry of dynamic tag
	void Automatizator::EnumerateDynamicTags()
	{
		first_num = last_num = 0;
		const ElementType paralleltypes = NODE | EDGE | FACE | CELL;
		for (index_enum::iterator it = index_tags.begin(); it != index_tags.end(); ++it)
		{
			for (ElementType etype = NODE; etype <= MESH; etype = etype << 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
629
				if (it->indices.isDefined(etype) && it->indices.isSparse(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
630
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
631
					for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
Kirill Terekhov's avatar
Kirill Terekhov committed
632
						jt->DelData(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
633
634
				}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
635
636
637


		for (index_enum::iterator it = index_tags.begin(); it != index_tags.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
638
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
639
			for (ElementType etype = NODE; etype <= MESH; etype = etype << 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
640
641
642
643
644
645
646
			{
				if (it->indices.isDefined(etype))
				{
					if (it->indices.GetSize() == ENUMUNDEF)
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
647
648
							for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
							if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
Kirill Terekhov's avatar
Kirill Terekhov committed
649
650
							{
								Storage::integer_array indarr = jt->IntegerArray(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
651
								indarr.resize(jt->RealArray(it->d.t).size());
Kirill Terekhov's avatar
Kirill Terekhov committed
652
								for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
Kirill Terekhov's avatar
Kirill Terekhov committed
653
									*qt = last_num++;
Kirill Terekhov's avatar
Kirill Terekhov committed
654
655
656
657
							}
						}
						else
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
658
659
							for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
							if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
Kirill Terekhov's avatar
Kirill Terekhov committed
660
661
662
663
							{
								Storage::integer_array indarr = jt->IntegerArray(it->indices);
								indarr.resize(jt->RealArray(it->d.t).size());
								for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
Kirill Terekhov's avatar
Kirill Terekhov committed
664
									*qt = last_num++;
Kirill Terekhov's avatar
Kirill Terekhov committed
665
666
667
668
669
670
671
							}
						}
					}
					else
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
672
673
							for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
							if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
Kirill Terekhov's avatar
Kirill Terekhov committed
674
675
676
							{
								Storage::integer_array indarr = jt->IntegerArray(it->indices);
								for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
Kirill Terekhov's avatar
Kirill Terekhov committed
677
									*qt = last_num++;
Kirill Terekhov's avatar
Kirill Terekhov committed
678
679
680
681
							}
						}
						else
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
682
683
							for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
							if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
Kirill Terekhov's avatar
Kirill Terekhov committed
684
685
686
							{
								Storage::integer_array indarr = jt->IntegerArray(it->indices);
								for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
Kirill Terekhov's avatar
Kirill Terekhov committed
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
									*qt = last_num++;
							}
						}
					}
				}
			}
		}
#if defined(USE_MPI)
		if (m->GetProcessorsNumber() > 1)
		{
			MPI_Scan(&last_num, &first_num, 1, INMOST_MPI_DATA_ENUM_TYPE, MPI_SUM, m->GetCommunicator());
			first_num -= last_num;
			ElementType exch_mask = NONE;
			for (index_enum::iterator it = index_tags.begin(); it != index_tags.end(); ++it)
			{
				for (ElementType etype = NODE; etype <= MESH; etype = etype << 1)
				{
					if (it->indices.isDefined(etype))
					{
						exch_mask |= etype;
						if (it->indices.GetSize() == ENUMUNDEF)
						{
							if (!it->indices.isSparse(etype))
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
711
712
								for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
Kirill Terekhov's avatar
Kirill Terekhov committed
713
714
715
716
717
718
719
720
								{
									Storage::integer_array indarr = jt->IntegerArray(it->indices);
									for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										*qt += first_num;
								}
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
721
722
								for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
Kirill Terekhov's avatar
Kirill Terekhov committed
723
724
725
726
727
728
729
730
731
732
733
734
								{
									Storage::integer_array indarr = jt->IntegerArray(it->indices);
									indarr.resize(jt->RealArray(it->d.t).size());
									for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										*qt += first_num;
								}
							}
						}
						else
						{
							if (!it->indices.isSparse(etype))
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
735
736
								for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
Kirill Terekhov's avatar
Kirill Terekhov committed
737
738
739
740
741
742
743
744
								{
									Storage::integer_array indarr = jt->IntegerArray(it->indices);
									for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										*qt += first_num;
								}
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
745
746
								for (Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt)
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
Kirill Terekhov's avatar
Kirill Terekhov committed
747
748
749
750
751
								{
									Storage::integer_array indarr = jt->IntegerArray(it->indices);
									for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										*qt += first_num;
								}
Kirill Terekhov's avatar
Kirill Terekhov committed
752
753
754
755
756
757
758
759
760
							}
						}
					}
				}
			}
			last_num += first_num;
			{
				std::vector<Tag> exch_tags;
				for (index_enum::iterator it = index_tags.begin(); it != index_tags.end(); ++it) exch_tags.push_back(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
761
				m->ExchangeData(exch_tags, exch_mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
762
763
764
765
766
767
768
			}
		}
#endif

	}
	/// register tag, data for which don't change through iterations
	/// don't register tag twice
Kirill Terekhov's avatar
Kirill Terekhov committed
769
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterStaticTag(Tag t, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
770
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
771
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_ctags.size()) + AD_CTAG;
Kirill Terekhov's avatar
Kirill Terekhov committed
772
773
774
		tagdomain d;
		d.t = t;
		d.domain_mask = domain_mask;
Kirill Terekhov's avatar
Kirill Terekhov committed
775
		reg_ctags.push_back(d);
Kirill Terekhov's avatar
Kirill Terekhov committed
776
777
778
		return ret;
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
779
	INMOST_DATA_ENUM_TYPE Automatizator::GetStencil(INMOST_DATA_ENUM_TYPE stnclind, const Storage & elem, void * user_data, stencil_pairs & ret)
Kirill Terekhov's avatar
Kirill Terekhov committed
780
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
781
		stencil_kind_domain & st = reg_stencils[stnclind-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
782
783
		assert(st.domainmask == 0 || elem->GetMarker(st.domainmask));
		if (st.kind == 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
784
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
785
786
787
			Storage::reference_array elems = elem->ReferenceArray(static_cast<stencil_tag *>(st.link)->elements);
			Storage::real_array coefs = elem->RealArray(static_cast<stencil_tag *>(st.link)->coefs);
			assert(elems.size() == coefs.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
788
789
			for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) 
				ret.push_back(std::make_pair(elems.at(k), coefs[k]));
Kirill Terekhov's avatar
Kirill Terekhov committed
790
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
791
792
		else if (st.kind == 1) reinterpret_cast<stencil_callback>(st.link)(elem, ret, user_data);
		return static_cast<INMOST_DATA_ENUM_TYPE>(ret.size());
Kirill Terekhov's avatar
Kirill Terekhov committed
793
794
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
795
	INMOST_DATA_ENUM_TYPE Automatizator::table::binary_search(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
796
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
797
798
		int l = 0, r = static_cast<int>(size)-1, mid = 0;
		while (r >= l)
Kirill Terekhov's avatar
Kirill Terekhov committed
799
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
800
801
802
803
			mid = (l + r) / 2;
			if (args[mid] > arg) r = mid - 1;
			else if (args[mid] < arg) l = mid + 1;
			else return mid;
Kirill Terekhov's avatar
Kirill Terekhov committed
804
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
805
806
807
		mid = (l + r) / 2;
		if (mid > static_cast<int>(size - 2)) mid = static_cast<int>(size - 2);
		return static_cast<INMOST_DATA_ENUM_TYPE>(mid);
Kirill Terekhov's avatar
Kirill Terekhov committed
808
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
809
	INMOST_DATA_REAL_TYPE Automatizator::table::get_value(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
810
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
811
812
813
		if (arg < args[0]) return vals[0];
		INMOST_DATA_ENUM_TYPE i = binary_search(arg);
		return vals[i] + (vals[i + 1] - vals[i]) * (arg - args[i]) / (args[i + 1] - args[i]);
Kirill Terekhov's avatar
Kirill Terekhov committed
814
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
815
	INMOST_DATA_REAL_TYPE Automatizator::table::get_derivative(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
816
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
817
818
819
		if (arg < args[0]) return 0.0;
		INMOST_DATA_ENUM_TYPE i = binary_search(arg);
		return (vals[i + 1] - vals[i]) / (args[i + 1] - args[i]);
Kirill Terekhov's avatar
Kirill Terekhov committed
820
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
821
	std::pair<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> Automatizator::table::get_both(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
822
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
823
824
825
826
		if (arg < args[0]) return std::make_pair(vals[0], 0.0);
		INMOST_DATA_ENUM_TYPE i = binary_search(arg);
		INMOST_DATA_REAL_TYPE der = (vals[i + 1] - vals[i]) / (args[i + 1] - args[i]);
		return std::make_pair(vals[i] + der * (arg - args[i]), der);
Kirill Terekhov's avatar
Kirill Terekhov committed
827
828
829
830
	}
};

#endif //USE_AUTODIFF