autodiff.cpp 35.3 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 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
  Automatizator * Automatizator::CurrentAutomatizator = NULL;
  bool print_ad_ctor = false;
  bool GetAutodiffPrint() {return print_ad_ctor;}
  void SetAutodiffPrint(bool set) {print_ad_ctor = set;}
  bool CheckCurrentAutomatizator() {return Automatizator::HaveCurrent();}

  void multivar_expression::FromBasicExpression(const basic_expression & expr)
  {
    Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
    expr.GetDerivative(1.0,merger);
    merger.RetriveRow(entries);
    merger.Clear();
  }

  void multivar_expression::AddBasicExpression(INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr)
  {
    Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
    merger.PushRow(multme,entries);
    expr.GetDerivative(multit,merger);
    merger.RetriveRow(entries);
    merger.Clear();
  }

  void multivar_expression::FromGetDerivative(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
  {
    Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
    GetDerivative(mult,merger);
    merger.AddRow(1.0,r);
    merger.RetriveRow(r);
    merger.Clear();
  }

Kirill Terekhov's avatar
Kirill Terekhov committed
48
#if defined(NEW_VERSION)
Kirill Terekhov's avatar
Kirill Terekhov committed
49
  //! returns offset from the end of precomputed z
Kirill Terekhov's avatar
Kirill Terekhov committed
50
	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
51 52
	{
		INMOST_DATA_ENUM_TYPE voffset = var.values_offset(element), doffset = var.derivatives_offset(element);
Kirill Terekhov's avatar
Kirill Terekhov committed
53
		INMOST_DATA_ENUM_TYPE k = static_cast<INMOST_DATA_ENUM_TYPE>(var.data.size()-1);
Kirill Terekhov's avatar
fixes  
Kirill Terekhov committed
54
		Storage e = Storage(m,var.current_stencil[element].first);
Kirill Terekhov's avatar
Kirill Terekhov committed
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
		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
133
					next.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
134 135 136 137 138 139 140 141 142
					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
143
					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
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166
				}
				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
167
		Storage e = Storage(m,var.current_stencil[element].first);
Kirill Terekhov's avatar
Kirill Terekhov committed
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
		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
183
					next.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201
					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
202
			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
203 204
			case AD_VAL: var.values[offset+k] = var.values[offset+it->left.i]; break;
			default:
Kirill Terekhov's avatar
Kirill Terekhov committed
205 206
				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
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
				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
225
	INMOST_DATA_REAL_TYPE Automatizator::Evaluate(expr & var, const Storage & e, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
226 227
	{
		var.current_stencil.resize(1);
Kirill Terekhov's avatar
fixes  
Kirill Terekhov committed
228
		var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
229 230 231 232
		var.resize_for_stencil();
		return EvaluateSub(var,0,ENUMUNDEF,user_data);
	}
	
Kirill Terekhov's avatar
Kirill Terekhov committed
233
	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
234
	{
235
    Sparse::RowMerger & m = GetMerger();
Kirill Terekhov's avatar
Kirill Terekhov committed
236 237
		INMOST_DATA_REAL_TYPE ret;
		var.current_stencil.resize(1);
Kirill Terekhov's avatar
fixes  
Kirill Terekhov committed
238
		var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
Kirill Terekhov's avatar
Kirill Terekhov committed
239 240
		var.resize_for_stencil();
		ret = EvaluateSub(var,0,ENUMUNDEF,user_data);
241 242 243 244
    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
245 246 247
		return ret*multiply;
	}
#else
Kirill Terekhov's avatar
Kirill Terekhov committed
248

Kirill Terekhov's avatar
Kirill Terekhov committed
249
	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
250 251 252 253 254
	{
		assert(var.op != AD_NONE);
		INMOST_DATA_REAL_TYPE lval, rval, ret = 0.0;
		switch (var.op)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
255
		case AD_COND_TYPE:
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
256
			lval = (e->GetElementType() & (*(ElementType*)&var.left))? 1.0 : -1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
257 258
			return lval;
		case AD_COND_MARK:
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
259
			lval = e->GetMarker(*(MarkerType *)&var.left)? 1.0 : -1.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
260
			return lval;
Kirill Terekhov's avatar
Kirill Terekhov committed
261 262 263 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 293
		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
294 295 296 297 298
		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
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
		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
328
			m->GetGeometricData(e->GetHandle(), MEASURE, &ret);
Kirill Terekhov's avatar
Kirill Terekhov committed
329
			return ret*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
330 331
		case AD_VAL:
			lval = DerivativePrecompute(*var.left, e, values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
332
			rval = Evaluate(*var.right,e,user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
333 334
			values.push_back(rval);
			return lval*var.coef;
Kirill Terekhov's avatar
Kirill Terekhov committed
335 336 337
		}
		if (var.op >= AD_FUNC)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
338
			ret = reg_funcs[var.op-AD_FUNC].func(e, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
339 340 341 342 343 344
			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
345
			ret = reg_tables[var.op-AD_TABLE]->get_value(lval);
Kirill Terekhov's avatar
Kirill Terekhov committed
346 347 348 349
			return ret*var.coef;
		}
		else if (var.op >= AD_STNCL)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
350
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
351 352 353 354 355 356
			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
357
				for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
358 359 360 361 362 363 364 365 366
				{
					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
367
				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
368
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
369
					lval = DerivativePrecompute(*var.left, Storage(m,get_st[k].first), values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390
					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
391
	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
392 393 394 395 396
	{
		assert(var.op != AD_NONE);
		INMOST_DATA_REAL_TYPE lval, rval, ret;
		switch (var.op)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
397 398 399
		case AD_COND_MARK:
		case AD_COND_TYPE:
			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 426 427 428 429 430
		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
431 432 433 434
		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
435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460
		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
461 462 463 464
		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
465 466 467 468 469 470 471 472
		}
		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
473
			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
474 475 476 477
			return;
		}
		else if (var.op >= AD_STNCL)
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
478
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
479 480 481 482 483 484
			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
485
				for (INMOST_DATA_ENUM_TYPE k = elems.size(); k > 0; --k) if( elems.at(k-1) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
486 487 488 489 490 491
					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
492
				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
493
					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
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509
			}
			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
510
	INMOST_DATA_REAL_TYPE Automatizator::Evaluate(const expr & var, const Storage & e, void * user_data)
Kirill Terekhov's avatar
Kirill Terekhov committed
511
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
512 513
		assert(var.op != AD_NONE);
		switch (var.op)
Kirill Terekhov's avatar
Kirill Terekhov committed
514
		{
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
515 516
		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
517 518 519 520 521 522 523 524 525 526 527 528 529 530
		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
531
		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
532 533
		case AD_VAL:   return Evaluate(*var.left,e,user_data)*var.coef;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
534 535
		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
536 537 538
		if (var.op >= AD_STNCL)
		{
			INMOST_DATA_REAL_TYPE ret = 0.0;
Kirill Terekhov's avatar
Kirill Terekhov committed
539
			stencil_kind_domain & st = reg_stencils[var.op-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
540 541
			assert(st.domainmask == 0 || e->GetMarker(st.domainmask));
			if (st.kind == 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
542
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
543 544 545
				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
546
				for (INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) if( elems.at(k) != InvalidHandle() )
Kirill Terekhov's avatar
Kirill Terekhov committed
547
					ret += var.coef * Evaluate(*var.left, elems[k], user_data) * coefs[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
548
			}
Kirill Terekhov's avatar
Kirill Terekhov committed
549 550 551 552
			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
553
				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
554
					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
555 556
			}
			return ret;
Kirill Terekhov's avatar
Kirill Terekhov committed
557
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
558 559 560 561
		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
562 563
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
564
	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
565
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
566
    Sparse::RowMerger & m = GetMerger();
Kirill Terekhov's avatar
Kirill Terekhov committed
567 568 569
		INMOST_DATA_REAL_TYPE ret;
		precomp_values_t values;
		ret = DerivativePrecompute(var, e, values, user_data);
Kirill Terekhov's avatar
Kirill Terekhov committed
570 571 572 573
    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
574
		return ret*multiply;
Kirill Terekhov's avatar
Kirill Terekhov committed
575
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
576
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
577 578 579 580 581 582 583
	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
584 585 586
			delete[] (*it)->args;
			delete[] (*it)->vals;
			delete (*it);
Kirill Terekhov's avatar
Kirill Terekhov committed
587 588
		}
		for (stencil_type::iterator it = reg_stencils.begin(); it != reg_stencils.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
589 590
		if (it->kind == 0)
			delete static_cast<stencil_tag *>(it->link);
Kirill Terekhov's avatar
Kirill Terekhov committed
591 592 593
	}
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterFunc(std::string name, func_callback func)
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
594
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_funcs.size()) + AD_FUNC;
Kirill Terekhov's avatar
Kirill Terekhov committed
595 596 597
		func_name_callback v;
		v.name = name;
		v.func = func;
Kirill Terekhov's avatar
Kirill Terekhov committed
598
		reg_funcs.push_back(v);
Kirill Terekhov's avatar
Kirill Terekhov committed
599 600 601
		return ret;
	}
	//register stencil that can be got from tags
Kirill Terekhov's avatar
Kirill Terekhov committed
602
	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
603
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
604
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
Kirill Terekhov's avatar
Kirill Terekhov committed
605 606 607 608 609 610 611 612
		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
613
		reg_stencils.push_back(st);
Kirill Terekhov's avatar
Kirill Terekhov committed
614 615 616
		return ret;
	}
	//register stencil that can be got from function
Kirill Terekhov's avatar
Kirill Terekhov committed
617
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterStencil(std::string name, stencil_callback func, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
618
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
619
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_stencils.size()) + AD_STNCL;
Kirill Terekhov's avatar
Kirill Terekhov committed
620 621 622 623 624
		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
625
		reg_stencils.push_back(st);
Kirill Terekhov's avatar
Kirill Terekhov committed
626 627 628 629
		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
630
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_tables.size()) + AD_TABLE;
Kirill Terekhov's avatar
Kirill Terekhov committed
631 632 633 634 635 636 637
		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
638
		reg_tables.push_back(t);
Kirill Terekhov's avatar
Kirill Terekhov committed
639 640 641 642
		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
643
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterDynamicTag(Tag t, ElementType typemask, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
644 645 646 647 648 649 650 651 652 653 654
	{
		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
655 656
		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
657 658 659 660 661 662 663 664 665 666 667
		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
668
				if (it->indices.isDefined(etype) && it->indices.isSparse(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
669
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
670
					for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
Kirill Terekhov's avatar
Kirill Terekhov committed
671
						jt->DelData(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
672 673
				}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
674 675 676


		for (index_enum::iterator it = index_tags.begin(); it != index_tags.end(); ++it)
Kirill Terekhov's avatar
Kirill Terekhov committed
677
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
678
			for (ElementType etype = NODE; etype <= MESH; etype = etype << 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
679 680 681 682 683 684 685
			{
				if (it->indices.isDefined(etype))
				{
					if (it->indices.GetSize() == ENUMUNDEF)
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
686 687 688 689 690 691 692 693 694 695
							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
696 697 698
						}
						else
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
699 700 701 702 703 704 705 706 707 708
							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
709 710 711 712 713 714
						}
					}
					else
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
715 716 717 718 719 720 721 722 723
							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
724 725 726
						}
						else
						{
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)) && (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
736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757
						}
					}
				}
			}
		}
#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
758 759 760 761 762 763 764 765 766
								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
767 768 769
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
770 771 772 773 774 775 776 777 778 779
								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
780 781 782 783 784 785
							}
						}
						else
						{
							if (!it->indices.isSparse(etype))
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
786 787 788 789 790 791 792 793 794
								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
795 796 797
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
798 799 800 801 802 803 804 805 806
								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
807 808 809 810 811 812 813 814 815
							}
						}
					}
				}
			}
			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
816
				m->ExchangeData(exch_tags, exch_mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
817 818 819
			}
		}
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835
    // 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
836 837 838
	}
	/// register tag, data for which don't change through iterations
	/// don't register tag twice
Kirill Terekhov's avatar
Kirill Terekhov committed
839
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterStaticTag(Tag t, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
840
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
841
		INMOST_DATA_ENUM_TYPE ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_ctags.size()) + AD_CTAG;
Kirill Terekhov's avatar
Kirill Terekhov committed
842 843 844
		tagdomain d;
		d.t = t;
		d.domain_mask = domain_mask;
Kirill Terekhov's avatar
Kirill Terekhov committed
845
		reg_ctags.push_back(d);
Kirill Terekhov's avatar
Kirill Terekhov committed
846 847 848
		return ret;
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
849
	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
850
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
851
		stencil_kind_domain & st = reg_stencils[stnclind-AD_STNCL];
Kirill Terekhov's avatar
Kirill Terekhov committed
852 853
		assert(st.domainmask == 0 || elem->GetMarker(st.domainmask));
		if (st.kind == 0)
Kirill Terekhov's avatar
Kirill Terekhov committed
854
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
855 856 857
			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
858 859
			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
860
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
861 862
		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
863 864
	}

Kirill Terekhov's avatar
Kirill Terekhov committed
865
	INMOST_DATA_ENUM_TYPE Automatizator::table::binary_search(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
866
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
867 868
		int l = 0, r = static_cast<int>(size)-1, mid = 0;
		while (r >= l)
Kirill Terekhov's avatar
Kirill Terekhov committed
869
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
870 871 872 873
			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
874
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
875 876 877
		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
878
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
879
	INMOST_DATA_REAL_TYPE Automatizator::table::get_value(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
880
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
881 882 883
		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
884
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
885
	INMOST_DATA_REAL_TYPE Automatizator::table::get_derivative(INMOST_DATA_REAL_TYPE arg)
Kirill Terekhov's avatar
Kirill Terekhov committed
886
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
887 888 889
		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
890
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
891
	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
892
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
893 894 895 896
		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
897 898 899 900
	}
};

#endif //USE_AUTODIFF