autodiff.cpp 19.3 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
1
2
3
4
5
6
7
8
#include "inmost.h"

#if defined(USE_AUTODIFF)

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

9

Kirill Terekhov's avatar
Kirill Terekhov committed
10
11
12
13


namespace INMOST
{
Kirill Terekhov's avatar
Kirill Terekhov committed
14
	
15
16
17
18
19
20
	template<> Demote<INMOST_DATA_REAL_TYPE>::type    AbstractEntry::Access<INMOST_DATA_REAL_TYPE>   (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Value(e,pos);}
	template<> Demote<INMOST_DATA_INTEGER_TYPE>::type AbstractEntry::Access<INMOST_DATA_INTEGER_TYPE>(const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Index(e,pos);}
	template<> Demote<unknown>::type                  AbstractEntry::Access<unknown>                 (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Unknown(e,pos);}
	template<> Demote<variable>::type                 AbstractEntry::Access<variable>                (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Unknown(e,pos);}
	template<> Demote<hessian_variable>::type         AbstractEntry::Access<hessian_variable>        (const Storage& e, INMOST_DATA_ENUM_TYPE pos) const {return Unknown(e,pos);}
	
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
	template<>
	Matrix<Demote<INMOST_DATA_REAL_TYPE>::type, pool_array_t<Demote<INMOST_DATA_REAL_TYPE>::type> >
	AbstractEntry::Access<INMOST_DATA_REAL_TYPE>   (const Storage& e) const {return Value(e);}
	template<>
	Matrix<Demote<INMOST_DATA_INTEGER_TYPE>::type, pool_array_t<Demote<INMOST_DATA_INTEGER_TYPE>::type> >
	AbstractEntry::Access<INMOST_DATA_INTEGER_TYPE>(const Storage& e) const {return Index(e);}
	template<>
	Matrix<Demote<unknown>::type, pool_array_t<Demote<unknown>::type> >
	AbstractEntry::Access<unknown>(const Storage& e) const {return Unknown(e);}
	template<>
	Matrix<Demote<variable>::type, pool_array_t<Demote<variable>::type> >
	AbstractEntry::Access<variable>(const Storage& e) const {return Unknown(e);}
	template<>
	Matrix<Demote<hessian_variable>::type,pool_array_t<Demote<hessian_variable>::type> >
	AbstractEntry::Access<hessian_variable>(const Storage& e) const {return Unknown(e);}
36
	
Kirill Terekhov's avatar
Kirill Terekhov committed
37
	
38
#if defined(USE_MESH) //Automatizator class does not exist without mesh
39
40
41
42
43
	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();}
Kirill Terekhov's avatar
Kirill Terekhov committed
44

45
	void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr)
Kirill Terekhov's avatar
Kirill Terekhov committed
46
	{
47
48
49
50
		Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
		expr.GetJacobian(1.0,merger);
		merger.RetriveRow(entries);
		merger.Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
51
52
	}

53
	void AddBasicExpression(Sparse::Row & entries, INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr)
Kirill Terekhov's avatar
Kirill Terekhov committed
54
	{
55
56
57
58
59
		Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
		merger.PushRow(multme,entries);
		expr.GetJacobian(multit,merger);
		merger.RetriveRow(entries);
		merger.Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
60
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
61

62
	void FromGetJacobian(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r)
Kirill Terekhov's avatar
Kirill Terekhov committed
63
	{
64
65
66
67
68
		Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
		expr.GetJacobian(mult,merger);
		merger.AddRow(1.0,r);
		merger.RetriveRow(r);
		merger.Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
69
	}
70
	Sparse::RowMerger & GetCurrentMerger() {return Automatizator::GetCurrent()->GetMerger();}
Kirill Terekhov's avatar
Kirill Terekhov committed
71
72
73
74
75
76
77
78
#else //USE_MESH
	bool CheckCurrentAutomatizator() {return false;}
	void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr) {}
	void AddBasicExpression(Sparse::Row & entries, INMOST_DATA_REAL_TYPE multme, INMOST_DATA_REAL_TYPE multit, const basic_expression & expr) {}
	void FromGetJacobian(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) {}
#endif //USE_MESH
	
#if defined(USE_MESH)
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
79
80
	Automatizator::Automatizator(const Automatizator & b) : name(b.name+"_copy")
	{
81
		std::vector<INMOST_DATA_ENUM_TYPE> regs = b.ListRegisteredEntries();
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
82
		for(std::vector<INMOST_DATA_ENUM_TYPE>::iterator kt = regs.begin(); kt != regs.end(); ++kt)
83
84
			RegisterEntry(b.GetEntry(*kt));
		if( b.last_num != 0 ) EnumerateEntries();
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
85
86
87
88
89
90
	}
	Automatizator & Automatizator::operator =(Automatizator const & b)
	{
		if( &b != this )
		{
			name = b.name+"_copy";
91
92
93
94
95
96
			for (unsigned k = 0; k < reg_blocks.size(); k++)
				if( isRegisteredEntry(k) ) UnregisterEntry(k);
			del_blocks.clear();
			reg_blocks.clear();
			act_blocks.clear();
			std::vector<INMOST_DATA_ENUM_TYPE> regs = b.ListRegisteredEntries();
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
97
			for(std::vector<INMOST_DATA_ENUM_TYPE>::iterator kt = regs.begin(); kt != regs.end(); ++kt)
98
99
				RegisterEntry(b.GetEntry(*kt));
			if( b.last_num != 0 ) EnumerateEntries();
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
100
101
102
103
		}
		return *this;
	}
	Automatizator::Automatizator(std::string _name) :name(_name), first_num(0), last_num(0) {}
Kirill Terekhov's avatar
Kirill Terekhov committed
104
105
	Automatizator::~Automatizator()
	{
106
107
108
109
110
		for (unsigned k = 0; k < reg_blocks.size(); k++) 
			if( isRegisteredEntry(k) ) UnregisterEntry(k);
		del_blocks.clear();
		act_blocks.clear();
		reg_blocks.clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
111
	}
112
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterTag(Tag t, ElementType typemask, MarkerType domain_mask, bool inverse)
Kirill Terekhov's avatar
Kirill Terekhov committed
113
	{
114
		if( t.GetSize() == ENUMUNDEF )
115
			return RegisterEntry(VectorEntry(typemask,domain_mask,inverse,t));
116
		else if( t.GetSize() == 1 )
117
			return RegisterEntry(SingleEntry(typemask,domain_mask,inverse,t,0));
118
119
		else
		{
120
			BlockEntry b(typemask,domain_mask,inverse);
121
            for(INMOST_DATA_ENUM_TYPE k = 0; k < t.GetSize(); ++k)
122
123
124
125
126
127
128
129
				b.AddTag(t,k);
			return RegisterEntry(b);
		}
	}
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterEntry(const AbstractEntry & b)
	{
		Mesh * m = b.GetMeshLink();
		
130
		ElementType sparse = NONE;// b.GetElementType();
131
		for (ElementType q = NODE; q <= MESH; q = NextElementType(q)) if (q & b.GetElementType())
Kirill Terekhov's avatar
Kirill Terekhov committed
132
		{
133
			for(unsigned unk = 0; unk < b.Size(); ++unk)
ptomin's avatar
ptomin committed
134
				sparse |= b.GetValueTag(unk).isSparse(q) ? q : NONE;
Kirill Terekhov's avatar
Kirill Terekhov committed
135
		}
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
136
		
137
138
		INMOST_DATA_ENUM_TYPE ret = ENUMUNDEF;
		if( del_blocks.empty() )
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
139
		{
140
141
142
			ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_blocks.size());
			reg_blocks.push_back(b.Copy());
			act_blocks.push_back(true);
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
143
144
145
		}
		else
		{
146
147
148
149
150
151
152
153
			ret = del_blocks.back();
			assert(!act_blocks[ret]);
			del_blocks.pop_back();
			reg_blocks[ret] = b.Copy();
			act_blocks[ret] = true;
		}
		
		reg_blocks[ret]->reg_index = ret;
Kirill Terekhov's avatar
Kirill Terekhov committed
154
		//b.reg_index = ret;
155
156
157
158
		
		{
			std::stringstream tag_name;
			tag_name << name << "_BLK_" << ret << "_Offset";
159
			reg_blocks[ret]->SetOffsetTag(m->CreateTag(tag_name.str(),DATA_INTEGER,b.GetElementType() | MESH,sparse,1));
Kirill Terekhov's avatar
Kirill Terekhov committed
160
			//b.SetOffsetTag(reg_blocks[ret]->GetOffsetTag());
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
161
		}
162
					
Kirill Terekhov's avatar
Kirill Terekhov committed
163
164
		return ret;
	}
165
	
Kirill Terekhov's avatar
Kirill Terekhov committed
166
167
168
169
170
171
172
173
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterEntry(AbstractEntry & b)
	{
		INMOST_DATA_ENUM_TYPE ret = RegisterEntry(static_cast<const AbstractEntry &>(b));
		b.reg_index = reg_blocks[ret]->GetRegistrationIndex();
		b.SetOffsetTag(reg_blocks[ret]->GetOffsetTag());
		return ret;
	}
	
174
	void Automatizator::UnregisterEntry(INMOST_DATA_ENUM_TYPE ind)
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
175
	{
176
		assert(reg_blocks[ind]);
177
178
179
180
181
182
		if( reg_blocks[ind]->GetOffsetTag().isValid() ) 
			reg_blocks[ind]->GetOffsetTag().GetMeshLink()->DeleteTag(reg_blocks[ind]->GetOffsetTag());
		delete reg_blocks[ind];
		reg_blocks[ind] = NULL;
		del_blocks.push_back(ind);
		act_blocks[ind] = false;
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
183
	}
184
185
186
187
	
	void Automatizator::DeactivateEntry(INMOST_DATA_ENUM_TYPE ind)
	{
		assert(reg_blocks[ind] != NULL); ///This block was not deleted
188
189
190
191
192
		AbstractEntry & b = GetEntry(ind);
		Mesh * m = b.GetMeshLink();
		b.reg_index = ENUMUNDEF;
		//std::cout << "delete " << reg_blocks[ind]->GetOffsetTag().GetTagName() << " on " << ElementTypeName(reg_blocks[ind]->GetElementType()) << std::endl;
		m->DeleteTag(b.GetOffsetTag(),b.GetElementType());
193
194
195
196
197
198
		act_blocks[ind] = false;
	}
	
	void Automatizator::ActivateEntry(INMOST_DATA_ENUM_TYPE ind)
	{
		assert(reg_blocks[ind] != NULL); ///This block was not deleted
199
200
201
202
203
204
205
206
207
208
209
210
211
212
		AbstractEntry & b = GetEntry(ind);
		Mesh * m = b.GetMeshLink();
		b.reg_index = ind;
		{
			ElementType sparse = NONE;// b.GetElementType();
			for (ElementType q = NODE; q <= MESH; q = NextElementType(q)) if (q & b.GetElementType())
			{
				for(unsigned unk = 0; unk < b.Size(); ++unk)
					sparse |= b.GetValueTag(unk).isSparse(q) ? q : NONE;
			}
			std::stringstream tag_name;
			tag_name << name << "_BLK_" << ind << "_Offset";
			b.SetOffsetTag(m->CreateTag(tag_name.str(),DATA_INTEGER,b.GetElementType(),sparse,1));
		}
213
214
215
		act_blocks[ind] = true;
	}
	
216
	void Automatizator::EnumerateEntries()
Kirill Terekhov's avatar
Kirill Terekhov committed
217
218
219
	{
		first_num = last_num = 0;
		const ElementType paralleltypes = NODE | EDGE | FACE | CELL;
220
221
		
		for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
Kirill Terekhov's avatar
Kirill Terekhov committed
222
		{
223
224
225
			AbstractEntry & b = *reg_blocks[it];
			TagInteger offset_tag = b.GetOffsetTag();
			Mesh * m = offset_tag.GetMeshLink();
226
			for (ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
227
				if (offset_tag.isDefined(etype) && offset_tag.isSparse(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
228
				{
229
230
					for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) )
						m->ElementByLocalID(etype,kt).DelData(offset_tag);
Kirill Terekhov's avatar
Kirill Terekhov committed
231
232
				}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
233

234
		for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
Kirill Terekhov's avatar
Kirill Terekhov committed
235
		{
236
237
238
239
			AbstractEntry & b = *reg_blocks[it];
			TagInteger offset_tag = b.GetOffsetTag();
			Mesh * m = offset_tag.GetMeshLink();
			for (ElementType etype = MESH; etype >= NODE; etype = PrevElementType(etype)) if( b.GetElementType() & etype )
Kirill Terekhov's avatar
Kirill Terekhov committed
240
			{
241
				for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) )
Kirill Terekhov's avatar
Kirill Terekhov committed
242
				{
243
244
					Element jt = m->ElementByLocalID(etype,kt);
					if ((!(etype & paralleltypes) || (jt.GetStatus() != Element::Ghost)) && b.isValid(jt) && b.Size(jt))
Kirill Terekhov's avatar
Kirill Terekhov committed
245
					{
246
247
						offset_tag[jt] = last_num;
						last_num += b.Size(jt);
Kirill Terekhov's avatar
Kirill Terekhov committed
248
249
250
251
					}
				}
			}
		}
252
		//~ std::set<INMOST_DATA_ENUM_TYPE> Pre, Post; //Nonlocal indices
Kirill Terekhov's avatar
Kirill Terekhov committed
253
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
254
255
256
		int size;
		MPI_Comm_size(MPI_COMM_WORLD,&size);
		if (size > 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
257
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
258
			MPI_Scan(&last_num, &first_num, 1, INMOST_MPI_DATA_ENUM_TYPE, MPI_SUM, MPI_COMM_WORLD);
Kirill Terekhov's avatar
Kirill Terekhov committed
259
			first_num -= last_num;
260
			last_num += first_num;
Kirill Terekhov's avatar
Kirill Terekhov committed
261
			ElementType exch_mask = NONE;
262
			for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
Kirill Terekhov's avatar
Kirill Terekhov committed
263
			{
264
265
266
267
				AbstractEntry & b = *reg_blocks[it];
				TagInteger offset_tag = b.GetOffsetTag();
				Mesh * m = offset_tag.GetMeshLink();
				for (ElementType etype = MESH; etype >= NODE; etype = PrevElementType(etype)) if( b.GetElementType() & etype )
Kirill Terekhov's avatar
Kirill Terekhov committed
268
				{
269
270
					exch_mask |= etype;
					for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) )
Kirill Terekhov's avatar
Kirill Terekhov committed
271
					{
272
273
274
						Element jt = m->ElementByLocalID(etype,kt);
						if ((!(etype & paralleltypes) || (jt.GetStatus() != Element::Ghost)) && b.isValid(jt) && b.Size(jt))
							offset_tag[jt] += first_num;
Kirill Terekhov's avatar
Kirill Terekhov committed
275
276
277
					}
				}
			}
278
			//synchronize indices
Kirill Terekhov's avatar
Kirill Terekhov committed
279
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
280
				std::map<Mesh *,std::vector<Tag> > exch_tags;
281
282
				for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
					exch_tags[reg_blocks[it]->GetOffsetTag().GetMeshLink()].push_back(reg_blocks[it]->GetOffsetTag());
Kirill Terekhov's avatar
Kirill Terekhov committed
283
284
				for(std::map<Mesh *,std::vector<Tag> >::iterator it = exch_tags.begin(); it != exch_tags.end(); ++it)
					it->first->ExchangeData(it->second, exch_mask,0);
Kirill Terekhov's avatar
Kirill Terekhov committed
285
			}
286
#if 0
287
			//compute out-of-bounds indices
288
			for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
289
			{
290
291
292
293
				AbstractEntry & b = *reg_blocks[it];
				TagInteger offset_tag = b.GetOffsetTag();
				Mesh * m = offset_tag.GetMeshLink();
				for (ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype)) if( b.GetElementType() & etype )
294
				{
295
					for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) )
296
					{
297
298
						Element jt = m->ElementByLocalID(etype,kt);
						if ((!(etype & paralleltypes) || (jt.GetStatus() == Element::Ghost)) && b.isValid(jt) && b.Size(jt))
299
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
300
301
302
303
304
305
306
307
308
							for(unsigned q = 0; q < b.MatrixSize(jt); ++q) 
							{
								INMOST_DATA_ENUM_TYPE ind =  b.Index(jt,q);
								if( ind != ENUMUNDEF ) 
								{
									if( ind < first_num ) Pre.insert(ind);
									if( ind >= last_num ) Post.insert(ind);
								}
							}
309
						}
310
					}
311
312
				} //etype
			} //it
313
#endif
314
			// after cycle
Kirill Terekhov's avatar
Kirill Terekhov committed
315
316
		}
#endif
317
318
		//INMOST_DATA_INTEGER_TYPE max_unknowns = m->AggregateMax(static_cast<INMOST_DATA_INTEGER_TYPE>(last_num));
		//std::cout << "Proc " << m->GetProcessorRank() << " size " << last_num-first_num <<  " pre " << Pre.size() << " post " << Post.size() << " max " << max_unknowns << std::endl;
Kirill Terekhov's avatar
Kirill Terekhov committed
319
320
#if defined(USE_OMP)
#pragma omp parallel
321
#endif //USE_OMP
322
		{
323
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
324
#pragma omp single
325
#endif //USE_OMP
326
			{
327
				merger.resize(MAX_THREADS);
328
			}
329
330
			//~ merger[OMP_THREAD].Resize(first_num,last_num,std::vector<INMOST_DATA_ENUM_TYPE>(Pre.begin(),Pre.end()),std::vector<INMOST_DATA_ENUM_TYPE>(Post.begin(),Post.end()),false);
			merger[OMP_THREAD].Resize(first_num,last_num,false);
331
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
332
	}
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
333
	
334
	std::vector<INMOST_DATA_ENUM_TYPE> Automatizator::ListRegisteredEntries() const
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
335
336
	{
		std::vector<INMOST_DATA_ENUM_TYPE> ret;
337
338
339
        for(blk_enum::size_type it = 0; it < reg_blocks.size(); ++it)
            if( isRegisteredEntry(static_cast<INMOST_DATA_ENUM_TYPE>(it)) )
                ret.push_back(static_cast<INMOST_DATA_ENUM_TYPE>(it));
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
340
341
342
		return ret;
	}
	
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
	void BlockEntry::AddTag(Tag value, INMOST_DATA_ENUM_TYPE comp)
	{
		assert(unknown_tags.empty() || GetMeshLink() == value.GetMeshLink());
		if( comp == ENUMUNDEF )
		{
			if( value.GetSize() != ENUMUNDEF )
			{
				for(unsigned k = 0; k < value.GetSize(); ++k)
				{
					unknown_tags.push_back(value);
					unknown_comp.push_back(k);
				}
			}
			else throw "Cannot add variable-sized tag to block";
		}
		else
		{
			unknown_tags.push_back(value);
			unknown_comp.push_back(comp);
		}
	}
	
	
	INMOST_DATA_ENUM_TYPE MultiEntry::MatrixSize(const Storage & e) const 
	{
		INMOST_DATA_ENUM_TYPE ret = 0; 
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
			ret += entries[k]->MatrixSize(e); 
		return ret;
	}
	
	INMOST_DATA_REAL_TYPE MultiEntry::Value(const Storage & e, INMOST_DATA_ENUM_TYPE unk) const
	{
376
        unsigned pos = 0;
377
378
379
380
381
382
383
384
385
386
387
388
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
		{
			unsigned s = entries[k]->MatrixSize(e);
			if( pos + s > unk )
				return entries[k]->Value(e,unk-pos);
			else pos += s; 
		}
		throw Impossible;
	}
	
	INMOST_DATA_REAL_TYPE & MultiEntry::Value(const Storage & e, INMOST_DATA_ENUM_TYPE unk)
	{
389
        unsigned pos = 0;
390
391
392
393
394
395
396
397
398
399
400
401
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
		{
			unsigned s = entries[k]->MatrixSize(e);
			if( pos + s > unk )
				return entries[k]->Value(e,unk-pos);
			else pos += s; 
		}
		throw Impossible;
	}
	
	INMOST_DATA_ENUM_TYPE MultiEntry::Index(const Storage & e, INMOST_DATA_ENUM_TYPE unk) const
	{
402
        unsigned pos = 0;
403
404
405
406
407
408
409
410
411
412
413
414
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
		{
			unsigned s = entries[k]->MatrixSize(e);
			if( pos + s > unk )
				return entries[k]->Index(e,unk-pos);
			else pos += s; 
		}
		throw Impossible;
	}
	
	unknown MultiEntry::Unknown(const Storage & e, INMOST_DATA_ENUM_TYPE unk) const
	{
415
        unsigned pos = 0;
416
417
418
419
420
421
422
423
424
425
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
		{
			unsigned s = entries[k]->MatrixSize(e);
			if( pos + s > unk )
				return entries[k]->Unknown(e,unk-pos);
			else pos += s; 
		}
		throw Impossible;
	}
	
426
	rpMatrix MultiEntry::Value(const Storage & e) const
427
	{
428
		rpMatrix ret(MatrixSize(e),1);
429
430
431
432
433
434
435
436
437
438
		unsigned l = 0, r, t;
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
		{
			t = entries[k]->MatrixSize(e);
			for(r = 0; r < t; ++r)
				ret(l++,0) = entries[k]->Value(e,r);
		}
		return ret;
	}
	
439
	ipMatrix MultiEntry::Index(const Storage & e) const
440
	{
441
		ipMatrix ret(MatrixSize(e),1);
442
443
444
445
446
447
448
449
450
451
		unsigned l = 0, r, t;
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
		{
			t = entries[k]->MatrixSize(e);
			for(r = 0; r < t; ++r)
				ret(l++,0) = entries[k]->Index(e,r);
		}
		return ret;
	}
	
452
	upMatrix MultiEntry::operator [](const Storage & e) const
453
	{
454
		upMatrix ret(MatrixSize(e),1);
455
456
457
458
459
460
461
462
463
464
465
466
467
		unsigned l = 0, r, t;
		for(unsigned k = 0; k < entries.size(); ++k) if( entries[k]->isValid(e) )
		{
			t = entries[k]->MatrixSize(e);
			for(r = 0; r < t; ++r)
				ret(l++,0) = entries[k]->Unknown(e,r);
		}
		return ret;
	}

	INMOST_DATA_ENUM_TYPE MultiEntry::GetValueComp(INMOST_DATA_ENUM_TYPE unk) const
	{
		unsigned pos = 0, k = 0;
468
		while( pos + entries[k]->Size() <= unk ) pos += entries[k++]->Size();
469
470
471
472
473
474
475
		assert(k < entries.size());
		return entries[k]->GetValueComp(unk-pos);
	}
	
	TagRealArray MultiEntry::GetValueTag(INMOST_DATA_ENUM_TYPE unk) const
	{
		unsigned pos = 0, k = 0;
476
		while( pos + entries[k]->Size() <= unk ) pos += entries[k++]->Size();
477
478
479
480
481
482
		assert(k < entries.size());
		return entries[k]->GetValueTag(unk-pos);
	}
	
	AbstractEntry * MultiEntry::Copy() const
	{
483
		MultiEntry * ret = new MultiEntry(GetElementType(),GetMask(),GetMaskInverse());
484
485
486
487
488
489
		for(unsigned k = 0; k < entries.size(); ++k)
			ret->entries.push_back(entries[k]->Copy());
		return ret;
	}
	
	void MultiEntry::AddEntry(const AbstractEntry & entry)
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
490
	{
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
		assert(entries.empty() || (GetMeshLink() == entry.GetMeshLink())); 
		SetElementType(GetElementType() | entry.GetElementType());
		entries.push_back(entry.Copy());
	}
	
	
	void StatusBlockEntry::AddTag(Tag value, INMOST_DATA_ENUM_TYPE comp)
	{
		assert(unknown_tags.empty() || GetMeshLink() == value.GetMeshLink());
		if( comp == ENUMUNDEF )
		{
			if( value.GetSize() != ENUMUNDEF )
			{
				for(unsigned k = 0; k < value.GetSize(); ++k)
				{
					unknown_tags.push_back(value);
					unknown_comp.push_back(k);
				}
			}
			else throw "Cannot add variable-sized tag to block";
		}
		else
		{
			unknown_tags.push_back(value);
			unknown_comp.push_back(comp);
		}
	}
	
	INMOST_DATA_ENUM_TYPE StatusBlockEntry::Size(const Storage & e) const 
	{
		INMOST_DATA_ENUM_TYPE ret = 0; 
		int stat = status_tag[e];
		for(unsigned k = 0; k < unknown_tags.size(); ++k) 
524
			if( e.HaveData(unknown_tags[k]) && status_tbl[stat][k] ) ret++;
525
526
527
528
529
		return ret;
	}
	
	AbstractEntry * StatusBlockEntry::Copy() const 
	{
530
		StatusBlockEntry * ret = new StatusBlockEntry(GetElementType(),GetMask(),GetMaskInverse(),status_tag,status_tbl);
531
532
533
		for(unsigned k = 0; k < Size(); ++k) 
			ret->AddTag(unknown_tags[k],unknown_comp[k]); 
		return ret; 
Kirill Terekhov's avatar
Updates    
Kirill Terekhov committed
534
	}
535
536
	
	void AbstractEntry::SynchronizeData()
537
	{
538
		//synchronize indices
539
		{
540
541
542
543
544
545
546
547
548
549
			Mesh * m = NULL;
			std::set<Tag> exch_tags;
			for(unsigned jt = 0; jt < Size(); ++jt)
			{
				if( m == NULL )
					m = GetValueTag(jt).GetMeshLink();
				else assert(m == GetValueTag(jt).GetMeshLink());
				exch_tags.insert(GetValueTag(jt));
			}
			m->ExchangeData(std::vector<Tag>(exch_tags.begin(),exch_tags.end()), GetElementType(),0);
550
551
		}
	}
552
553
554
	
	
	void Automatizator::SynchronizeData()
555
	{
556
557
558
		//TODO:
		// optimize std::map usage
		// don't sort Tag by itself, only by name, otherwise memory location may affect order and result
559
		//synchronize indices
560
		{
561
			std::map<Mesh *,std::map<std::string,Tag> > exch_tags;
562
563
			ElementType exch_mask = NONE;
			for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
564
			{
565
566
				exch_mask |= reg_blocks[it]->GetElementType();
				for(unsigned jt = 0; jt < reg_blocks[it]->Size(); ++jt)
567
					exch_tags[reg_blocks[it]->GetValueTag(jt).GetMeshLink()][reg_blocks[it]->GetValueTag(jt).GetTagName()] = reg_blocks[it]->GetValueTag(jt);
568
			}
569
			for(std::map<Mesh *,std::map<std::string,Tag> >::iterator it = exch_tags.begin(); it != exch_tags.end(); ++it)
570
			{
571
572
573
574
575
576
577
578
579
580
				std::vector<Tag> sync;
				//std::cout << "Synchronize tags: ";
				for(std::map<std::string,Tag>::iterator jt = it->second.begin(); jt != it->second.end(); ++jt)
				{
					sync.push_back(jt->second);
					//std::cout << jt->first << " ";
				}
				//std::cout << std::endl;
				
				it->first->ExchangeData(sync, exch_mask,0);
581
582
583
			}
		}
	}
584
585
#endif //USE_MESH
} //namespace INMOST
Kirill Terekhov's avatar
Kirill Terekhov committed
586
587

#endif //USE_AUTODIFF