autodiff.cpp 34.2 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
#if defined(NEW_VERSION)
	
	//! returns offset from the end of precomputed z
Kirill Terekhov's avatar
Kirill Terekhov committed
19
	void Automatizator::DerivativeFill(expr & var, INMOST_DATA_ENUM_TYPE element, INMOST_DATA_ENUM_TYPE parent, Sparse::RowMerger & entries, INMOST_DATA_REAL_TYPE multval, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
20
21
	{
		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
Kirill Terekhov committed
202
	INMOST_DATA_REAL_TYPE Automatizator::Derivative(expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
203
	{
204
    Sparse::RowMerger & m = GetMerger();
Kirill Terekhov's avatar
Kirill Terekhov committed
205
206
		INMOST_DATA_REAL_TYPE ret;
		var.current_stencil.resize(1);
Kirill Terekhov's avatar
fixes    
Kirill Terekhov committed
207
		var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
208
209
		var.resize_for_stencil();
		ret = EvaluateSub(var,0,ENUMUNDEF,user_data);
210
211
212
213
    m.PushRow(1.0,out);
		DerivativeFill(var, 0, ENUMUNDEF, m, multiply, user_data);
    m.RetriveRow(out);
    m.Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
214
215
216
		return ret*multiply;
	}
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
217

Kirill Terekhov's avatar
Kirill Terekhov committed
218
	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
219
220
221
222
223
	{
		assert(var.op != AD_NONE);
		INMOST_DATA_REAL_TYPE lval, rval, ret = 0.0;
		switch (var.op)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
224
		case AD_COND_TYPE:
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
225
			lval = (e->GetElementType() & (*(ElementType*)&var.left))? 1.0 : -1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
226
227
			return lval;
		case AD_COND_MARK:
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
228
			lval = e->GetMarker(*(MarkerType *)&var.left)? 1.0 : -1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
229
			return lval;
Kirill Terekhov's avatar
Kirill Terekhov committed
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
259
260
261
262
		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
263
264
265
266
267
		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
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
293
294
295
296
		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
297
			m->GetGeometricData(e->GetHandle(), MEASURE, &ret);
Kirill Terekhov's avatar
Kirill Terekhov committed
298
			return ret*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
299
300
		case AD_VAL:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
301
			rval = Evaluate(*var.right,e,user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
302
303
			values.push_back(rval);
			return lval*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
304
305
306
		}
		if (var.op >= AD_FUNC)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
307
			ret = reg_funcs[var.op-AD_FUNC].func(e, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
308
309
310
311
312
313
			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
314
			ret = reg_tables[var.op-AD_TABLE]->get_value(lval);
Kirill Terekhov's avatar
Kirill Terekhov committed
315
316
317
318
			return ret*var.coef;
		}
		else if (var.op >= AD_STNCL)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
319
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
320
321
322
323
324
325
			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
326
				for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
327
328
329
330
331
332
333
334
335
				{
					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
336
				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
337
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
338
					lval = DerivativePrecompute(*var.left, Storage(m,get_st[k].first), values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
					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
360
	void Automatizator::DerivativeFill(const expr & var, const Storage & e, Sparse::RowMerger & entries, precomp_values_t & values, INMOST_DATA_REAL_TYPE multval, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
361
362
363
364
365
	{
		assert(var.op != AD_NONE);
		INMOST_DATA_REAL_TYPE lval, rval, ret;
		switch (var.op)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
366
367
368
		case AD_COND_MARK:
		case AD_COND_TYPE:
			return;
Kirill Terekhov's avatar
Kirill Terekhov committed
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
396
397
398
399
		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
400
401
402
403
		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
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
		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
430
431
432
433
		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
434
435
436
437
438
439
440
441
		}
		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
442
			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
443
444
445
446
			return;
		}
		else if (var.op >= AD_STNCL)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
447
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
448
449
450
451
452
453
			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
454
				for (INMOST_DATA_ENUM_TYPE k = elems.size(); k > 0; --k) if( elems.at(k-1) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
455
456
457
458
459
460
					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
461
				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
462
					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
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
			}
			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
479
	INMOST_DATA_REAL_TYPE Automatizator::Evaluate(const expr & var, const Storage & e, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
480
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
481
482
		assert(var.op != AD_NONE);
		switch (var.op)
Kirill Terekhov's avatar
Kirill Terekhov committed
483
		{
Kirill Terekhov's avatar
Fixes    
Kirill Terekhov committed
484
485
		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
486
487
488
489
490
491
492
493
494
495
496
497
498
499
		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
500
		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
501
502
		case AD_VAL:   return Evaluate(*var.left,e,user_data)*var.coef;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
503
504
		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
505
506
507
		if (var.op >= AD_STNCL)
		{
			INMOST_DATA_REAL_TYPE ret = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
508
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
509
510
			assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
			if (st.kind == 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
511
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
512
513
514
				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
515
				for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
516
					ret += var.coef * Evaluate(*var.left, elems[k], user_data) * coefs[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
517
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
518
519
520
521
			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
522
				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
523
					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
524
525
			}
			return ret;
Kirill Terekhov's avatar
Kirill Terekhov committed
526
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
527
528
529
530
		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
531
532
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
533
	INMOST_DATA_REAL_TYPE Automatizator::Derivative(const expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
534
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
535
    Sparse::RowMerger & m = GetMerger();
Kirill Terekhov's avatar
Kirill Terekhov committed
536
537
538
		INMOST_DATA_REAL_TYPE ret;
		precomp_values_t values;
		ret = DerivativePrecompute(var, e, values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
539
540
541
542
    m.PushRow(1.0,out);
		DerivativeFill(var, e, m, values, multiply, user_data);
    m.RetriveRow(out);
    m.Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
543
		return ret*multiply;
Kirill Terekhov's avatar
Kirill Terekhov committed
544
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
545
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
546
547
548
549
550
551
552
	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
553
554
555
			delete[] (*it)->args;
			delete[] (*it)->vals;
			delete (*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
556
557
		}
		for (stencil_type::iterator it = reg_stencils.begin(); it != reg_stencils.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
558
559
		if (it->kind == 0)
			delete static_cast<stencil_tag *>(it->link);
Kirill Terekhov's avatar
Kirill Terekhov committed
560
561
562
	}
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterFunc(std::string name, func_callback func)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
563
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_funcs.size()) + AD_FUNC;
Kirill Terekhov's avatar
Kirill Terekhov committed
564
565
566
		func_name_callback v;
		v.name = name;
		v.func = func;
Kirill Terekhov's avatar
Kirill Terekhov committed
567
		reg_funcs.push_back(v);
Kirill Terekhov's avatar
Kirill Terekhov committed
568
569
570
		return ret;
	}
	//register stencil that can be got from tags
Kirill Terekhov's avatar
Kirill Terekhov committed
571
	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
572
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
573
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
Kirill Terekhov's avatar
Kirill Terekhov committed
574
575
576
577
578
579
580
581
		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
582
		reg_stencils.push_back(st);
Kirill Terekhov's avatar
Kirill Terekhov committed
583
584
585
		return ret;
	}
	//register stencil that can be got from function
Kirill Terekhov's avatar
Kirill Terekhov committed
586
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterStencil(std::string name, stencil_callback func, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
587
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
588
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
Kirill Terekhov's avatar
Kirill Terekhov committed
589
590
591
592
593
		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
594
		reg_stencils.push_back(st);
Kirill Terekhov's avatar
Kirill Terekhov committed
595
596
597
598
		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
599
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_tables.size()) + AD_TABLE;
Kirill Terekhov's avatar
Kirill Terekhov committed
600
601
602
603
604
605
606
		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
607
		reg_tables.push_back(t);
Kirill Terekhov's avatar
Kirill Terekhov committed
608
609
610
611
		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
612
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterDynamicTag(Tag t, ElementType typemask, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
613
614
615
616
617
618
619
620
621
622
623
	{
		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
624
625
		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
626
627
628
629
630
631
632
633
634
635
636
		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
637
				if (it->indices.isDefined(etype) && it->indices.isSparse(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
638
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
639
					for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
Kirill Terekhov's avatar
Kirill Terekhov committed
640
						jt->DelData(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
641
642
				}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
643
644
645


		for (index_enum::iterator it = index_tags.begin(); it != index_tags.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
646
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
647
			for (ElementType etype = NODE; etype <= MESH; etype = etype << 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
648
649
650
651
652
653
654
			{
				if (it->indices.isDefined(etype))
				{
					if (it->indices.GetSize() == ENUMUNDEF)
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
655
656
657
658
659
660
661
662
663
664
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
              {
							  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
							  {
								  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 = last_num++;
							  }
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
665
666
667
						}
						else
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
668
669
670
671
672
673
674
675
676
677
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
              {
							  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
							  {
								  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 = last_num++;
							  }
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
678
679
680
681
682
683
						}
					}
					else
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
684
685
686
687
688
689
690
691
692
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
              {
							  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
							  {
								  Storage::integer_array indarr = jt->IntegerArray(it->indices);
								  for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
									  *qt = last_num++;
							  }
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
693
694
695
						}
						else
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
696
697
698
699
700
701
702
703
704
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
              {
							  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
							  {
								  Storage::integer_array indarr = jt->IntegerArray(it->indices);
								  for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
									  *qt = last_num++;
							  }
              }
Kirill Terekhov's avatar
Kirill Terekhov committed
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
						}
					}
				}
			}
		}
#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
727
728
729
730
731
732
733
734
735
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
                {
								  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
								  {
									  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
736
737
738
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
739
740
741
742
743
744
745
746
747
748
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
                {
                  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
								  {
									  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;
								  }
                }
Kirill Terekhov's avatar
Kirill Terekhov committed
749
750
751
752
753
754
							}
						}
						else
						{
							if (!it->indices.isSparse(etype))
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
755
756
757
758
759
760
761
762
763
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
                {
								  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask)))
								  {
									  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
764
765
766
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
767
768
769
770
771
772
773
774
775
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
                {
								  if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->d.t) && (it->d.domain_mask == 0 || jt->GetMarker(it->d.domain_mask))))
								  {
									  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
776
777
778
779
780
781
782
783
784
							}
						}
					}
				}
			}
			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
785
				m->ExchangeData(exch_tags, exch_mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
786
787
788
			}
		}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
    // this version will fail in parallel
    //merger.Resize(first_num,last_num,false);
    // use this version until there is a way to define multiple intervals in RowMerger
    INMOST_DATA_INTEGER_TYPE max_unknowns = m->AggregateMax(static_cast<INMOST_DATA_INTEGER_TYPE>(last_num));
#if defined(USE_OMP)
#pragma omp parallel
    {
#pragma omp single
      {
        merger.resize(omp_get_num_procs());
      }
      merger[omp_get_thread_num()].Resize(0,max_unknowns,false);
    }
#else
    merger.Resize(0,max_unknowns,false);
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
805
806
807
	}
	/// register tag, data for which don't change through iterations
	/// don't register tag twice
Kirill Terekhov's avatar
Kirill Terekhov committed
808
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterStaticTag(Tag t, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
809
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
810
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_ctags.size()) + AD_CTAG;
Kirill Terekhov's avatar
Kirill Terekhov committed
811
812
813
		tagdomain d;
		d.t = t;
		d.domain_mask = domain_mask;
Kirill Terekhov's avatar
Kirill Terekhov committed
814
		reg_ctags.push_back(d);
Kirill Terekhov's avatar
Kirill Terekhov committed
815
816
817
		return ret;
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
818
	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
819
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
820
		stencil_kind_domain & st = reg_stencils[stnclind-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
821
822
		assert(st.domainmask == 0 || elem->GetMarker(st.domainmask));
		if (st.kind == 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
823
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
824
825
826
			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
827
828
			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
829
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
830
831
		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
832
833
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
834
	INMOST_DATA_ENUM_TYPE Automatizator::table::binary_search(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
835
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
836
837
		int l = 0, r = static_cast<int>(size)-1, mid = 0;
		while (r >= l)
Kirill Terekhov's avatar
Kirill Terekhov committed
838
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
839
840
841
842
			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
843
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
844
845
846
		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
847
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
848
	INMOST_DATA_REAL_TYPE Automatizator::table::get_value(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
849
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
850
851
852
		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
853
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
854
	INMOST_DATA_REAL_TYPE Automatizator::table::get_derivative(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
855
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
856
857
858
		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
859
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
860
	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
861
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
862
863
864
865
		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
866
867
868
869
	}
};

#endif //USE_AUTODIFF