autodiff.cpp 17.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
	
#if defined(USE_MESH)
18 19 20 21 22
	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
23

24
	void FromBasicExpression(Sparse::Row & entries, const basic_expression & expr)
Kirill Terekhov's avatar
Kirill Terekhov committed
25
	{
26 27 28 29
		Sparse::RowMerger & merger = Automatizator::GetCurrent()->GetMerger();
		expr.GetJacobian(1.0,merger);
		merger.RetriveRow(entries);
		merger.Clear();
Kirill Terekhov's avatar
Kirill Terekhov committed
30 31
	}

32
	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
33
	{
34 35 36 37 38
		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
39
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
40

41
	void FromGetJacobian(const basic_expression & expr, INMOST_DATA_REAL_TYPE mult, Sparse::Row & r)
Kirill Terekhov's avatar
Kirill Terekhov committed
42
	{
43 44 45 46 47
		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
48
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
49 50 51 52 53 54 55 56
#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
Kirill Terekhov committed
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
	Automatizator::Automatizator(const Automatizator & b) : name(b.name+"_copy")
	{
		std::vector<INMOST_DATA_ENUM_TYPE> regs = b.ListRegisteredTags();
		for(std::vector<INMOST_DATA_ENUM_TYPE>::iterator kt = regs.begin(); kt != regs.end(); ++kt)
			RegisterTag(b.GetValueTag(*kt),b.GetElementType(*kt),b.GetMask(*kt));
		if( b.last_num != 0 ) EnumerateTags();
	}
	Automatizator & Automatizator::operator =(Automatizator const & b)
	{
		if( &b != this )
		{
			name = b.name+"_copy";
			del_tags.clear();
			reg_tags.clear();
			std::vector<INMOST_DATA_ENUM_TYPE> regs = b.ListRegisteredTags();
			for(std::vector<INMOST_DATA_ENUM_TYPE>::iterator kt = regs.begin(); kt != regs.end(); ++kt)
				RegisterTag(b.GetValueTag(*kt),b.GetElementType(*kt),b.GetMask(*kt));
			if( b.last_num != 0 ) EnumerateTags();
		}
		return *this;
	}
	Automatizator::Automatizator(std::string _name) :name(_name), first_num(0), last_num(0) {}
Kirill Terekhov's avatar
Kirill Terekhov committed
79 80
	Automatizator::~Automatizator()
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
81 82 83
		del_tags.clear();
		for (unsigned k = 0; k < reg_tags.size(); k++) if( reg_tags[k].active )
			reg_tags[k].indices = reg_tags[k].indices.GetMeshLink()->DeleteTag(reg_tags[k].indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
84
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
85
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterTag(Tag t, ElementType typemask, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
86
	{
Kirill Terekhov's avatar
Kirill Terekhov committed
87 88 89
		tagdata p;
		p.domain_mask = domain_mask;
		p.t = t;
Kirill Terekhov's avatar
Kirill Terekhov committed
90 91 92 93 94 95
		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;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
		p.indices = t.GetMeshLink()->CreateTag(t.GetTagName() + "_index_" + name, DATA_INTEGER, def, sparse, t.GetSize());
		p.active = true;
		INMOST_DATA_ENUM_TYPE ret;
		
		if( del_tags.empty() )
		{
			ret = static_cast<INMOST_DATA_ENUM_TYPE>(reg_tags.size());
			reg_tags.push_back(p);
		}
		else
		{
			ret = del_tags.back();
			assert(!reg_tags[ret].active);
			del_tags.pop_back();
			reg_tags[ret] = p;
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
112 113
		return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
114 115 116 117 118 119
	void Automatizator::UnregisterTag(INMOST_DATA_ENUM_TYPE ind)
	{
		assert(reg_tags[ind].active);
		del_tags.push_back(ind);
		reg_tags[ind].active = false;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
120
	void Automatizator::EnumerateTags()
Kirill Terekhov's avatar
Kirill Terekhov committed
121 122 123
	{
		first_num = last_num = 0;
		const ElementType paralleltypes = NODE | EDGE | FACE | CELL;
Kirill Terekhov's avatar
Kirill Terekhov committed
124
		for (tag_enum::iterator it = reg_tags.begin(); it != reg_tags.end(); ++it) if( it->active )
Kirill Terekhov's avatar
Kirill Terekhov committed
125
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
126
			Mesh * m = it->indices.GetMeshLink();
Kirill Terekhov's avatar
Kirill Terekhov committed
127
			for (ElementType etype = NODE; etype <= MESH; etype = etype << 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
128
				if (it->indices.isDefined(etype) && it->indices.isSparse(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
129
				{
Kirill Terekhov's avatar
Kirill Terekhov committed
130
					for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
Kirill Terekhov's avatar
Kirill Terekhov committed
131
						jt->DelData(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
132 133
				}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
134 135


Kirill Terekhov's avatar
Kirill Terekhov committed
136
		for (tag_enum::iterator it = reg_tags.begin(); it != reg_tags.end(); ++it) if( it->active )
Kirill Terekhov's avatar
Kirill Terekhov committed
137
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
138
			Mesh * m = it->indices.GetMeshLink();
Kirill Terekhov's avatar
Kirill Terekhov committed
139
			for (ElementType etype = MESH; etype >= NODE; etype = PrevElementType(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
140 141 142 143 144 145 146
			{
				if (it->indices.isDefined(etype))
				{
					if (it->indices.GetSize() == ENUMUNDEF)
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
147
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
148
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
149
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask)))
150 151
								{
									Storage::integer_array indarr = jt->IntegerArray(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
152
									indarr.resize(jt->RealArray(it->t).size());
153 154 155 156
									for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										*qt = last_num++;
								}
							}
Kirill Terekhov's avatar
Kirill Terekhov committed
157 158 159
						}
						else
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
160
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
161
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
162
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->t) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask))))
163 164
								{
									Storage::integer_array indarr = jt->IntegerArray(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
165
									indarr.resize(jt->RealArray(it->t).size());
166 167 168 169
									for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										*qt = last_num++;
								}
							}
Kirill Terekhov's avatar
Kirill Terekhov committed
170 171 172 173 174 175
						}
					}
					else
					{
						if (!it->indices.isSparse(etype))
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
176
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
177
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
178
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask)))
179 180 181 182 183 184
								{
									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
185 186 187
						}
						else
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
188
							for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
189
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
190
								if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->t) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask))))
191 192 193 194 195 196
								{
									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
197 198 199 200 201
						}
					}
				}
			}
		}
202
		std::set<INMOST_DATA_ENUM_TYPE> Pre, Post; //Nonlocal indices
Kirill Terekhov's avatar
Kirill Terekhov committed
203
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
204 205 206
		int size;
		MPI_Comm_size(MPI_COMM_WORLD,&size);
		if (size > 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
207
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
208
			MPI_Scan(&last_num, &first_num, 1, INMOST_MPI_DATA_ENUM_TYPE, MPI_SUM, MPI_COMM_WORLD);
Kirill Terekhov's avatar
Kirill Terekhov committed
209 210
			first_num -= last_num;
			ElementType exch_mask = NONE;
Kirill Terekhov's avatar
Kirill Terekhov committed
211
			for (tag_enum::iterator it = reg_tags.begin(); it != reg_tags.end(); ++it) if( it->active )
Kirill Terekhov's avatar
Kirill Terekhov committed
212
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
213
				Mesh * m = it->indices.GetMeshLink();
Kirill Terekhov's avatar
Kirill Terekhov committed
214
				for (ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
215 216 217 218 219 220 221 222
				{
					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
223
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
224
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
225
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask)))
226 227 228 229 230 231
									{
										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
232 233 234
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
235
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
236
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
237
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->t) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask))))
238 239
									{
										Storage::integer_array indarr = jt->IntegerArray(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
240
										indarr.resize(jt->RealArray(it->t).size());
241 242 243 244
										for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
											*qt += first_num;
									}
								}
Kirill Terekhov's avatar
Kirill Terekhov committed
245 246 247 248 249 250
							}
						}
						else
						{
							if (!it->indices.isSparse(etype))
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
251
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
252
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
253
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask)))
254 255 256 257 258 259
									{
										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
260 261 262
							}
							else
							{
Kirill Terekhov's avatar
Kirill Terekhov committed
263
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
264
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
265
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() != Element::Ghost)) && (jt->HaveData(it->t) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask))))
266 267 268 269 270 271
									{
										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
272 273 274 275 276
							}
						}
					}
				}
			}
277
			//synchronize indices
Kirill Terekhov's avatar
Kirill Terekhov committed
278 279
			last_num += first_num;
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
280
				std::map<Mesh *,std::vector<Tag> > exch_tags;
Kirill Terekhov's avatar
Kirill Terekhov committed
281
				for (tag_enum::iterator it = reg_tags.begin(); it != reg_tags.end(); ++it) if( it->active )
Kirill Terekhov's avatar
Kirill Terekhov committed
282 283 284
					exch_tags[it->indices.GetMeshLink()].push_back(it->indices);
				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
			//compute out-of-bounds indices
Kirill Terekhov's avatar
Kirill Terekhov committed
287
			for (tag_enum::iterator it = reg_tags.begin(); it != reg_tags.end(); ++it) if( it->active )
288
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
289
				Mesh * m = it->indices.GetMeshLink();
Kirill Terekhov's avatar
Kirill Terekhov committed
290
				for (ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
291 292 293 294 295 296 297 298 299 300
				{
					if (it->indices.isDefined(etype))
					{
						exch_mask |= etype;
						if (it->indices.GetSize() == ENUMUNDEF)
						{
							if (!it->indices.isSparse(etype))
							{
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
301
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() == Element::Ghost)) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask)))
302 303 304 305
									{
										Storage::integer_array indarr = jt->IntegerArray(it->indices);
										for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										{
306 307
											if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
											else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
308 309 310 311 312 313 314 315
										}
									}
								}
							}
							else
							{
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
316
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() == Element::Ghost)) && (jt->HaveData(it->t) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask))))
317 318
									{
										Storage::integer_array indarr = jt->IntegerArray(it->indices);
Kirill Terekhov's avatar
Kirill Terekhov committed
319
										indarr.resize(jt->RealArray(it->t).size());
320 321
										for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										{
322 323
											if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
											else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
324 325 326 327 328 329 330 331 332 333 334
										}
									}
								}
							}
						}
						else //getsize
						{
							if (!it->indices.isSparse(etype))
							{
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
335
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() == Element::Ghost)) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask)))
336 337 338 339
									{
										Storage::integer_array indarr = jt->IntegerArray(it->indices);
										for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										{
340 341
											if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
											else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
342 343 344 345 346 347 348 349
										}
									}
								}
							}
							else
							{
								for (Mesh::iteratorStorage jt = m->Begin(etype); jt != m->End(); ++jt)
								{
Kirill Terekhov's avatar
Kirill Terekhov committed
350
									if ((!(etype & paralleltypes) || ((etype & paralleltypes) && jt->getAsElement()->GetStatus() == Element::Ghost)) && (jt->HaveData(it->t) && (it->domain_mask == 0 || jt->GetMarker(it->domain_mask))))
351 352 353 354
									{
										Storage::integer_array indarr = jt->IntegerArray(it->indices);
										for (Storage::integer_array::iterator qt = indarr.begin(); qt != indarr.end(); ++qt)
										{
355 356
											if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) < first_num ) Pre.insert(*qt);
											else if( static_cast<INMOST_DATA_ENUM_TYPE>(*qt) >= last_num ) Post.insert(*qt);
357 358 359 360 361 362 363 364 365
										}
									}
								}
							}
						} //getsize
					} //isdefined
				} //etype
			} //it
			// after cycle
Kirill Terekhov's avatar
Kirill Terekhov committed
366 367
		}
#endif
368 369 370
		// 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
371 372
		//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
373 374
#if defined(USE_OMP)
#pragma omp parallel
375
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
376
#pragma omp single
377 378 379
			{
				merger.resize(omp_get_num_procs());
			}
380
			merger[omp_get_thread_num()].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);
381
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
382
#else
383
		merger.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);
Kirill Terekhov's avatar
Kirill Terekhov committed
384
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
385
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402
	
	std::vector<INMOST_DATA_ENUM_TYPE> Automatizator::ListRegisteredTags() const
	{
		std::vector<INMOST_DATA_ENUM_TYPE> ret;
		for(tag_enum::size_type it = 0; it < reg_tags.size(); ++it) if( reg_tags[it].active )
			ret.push_back(static_cast<INMOST_DATA_ENUM_TYPE>(it));
		return ret;
	}
	
	ElementType Automatizator::GetElementType(INMOST_DATA_ENUM_TYPE ind) const
	{
		Tag index = GetIndexTag(ind);
		ElementType ret = NONE;
		for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
			ret |= (index.isDefined(etype) ? etype : NONE);
			return ret;
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
403
#endif //USE_MESH
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 431 432 433 434 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 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483

#if defined(USE_SOLVER)
	void Residual::GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const 
	{
		start = residual.GetFirstIndex(); 
		end = residual.GetLastIndex();
	}
	void Residual::SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end)
	{
		jacobian.SetInterval(beg,end);
		residual.SetInterval(beg,end);
	}
	void Residual::ClearResidual()
	{
		for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0;
	}
	void Residual::ClearJacobian()
	{
		for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it) it->Clear();
	}
	void Residual::Clear()
	{
#if defined(USE_OMP)
#pragma omp for
#endif //USE_OMP
		for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) 
		{
			residual[k] = 0.0;
			jacobian[k].Clear();
		}
	}
	INMOST_DATA_REAL_TYPE Residual::Norm()
	{
		INMOST_DATA_REAL_TYPE ret = 0;
#if defined(USE_OMP)
#pragma omp parallel for reduction(+:ret)
#endif //USE_OMP
		for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) 
			ret += residual[k]*residual[k];
#if defined(USE_MPI)
		INMOST_DATA_REAL_TYPE tmp = ret;
		MPI_Allreduce(&tmp, &ret, 1, INMOST_MPI_DATA_REAL_TYPE, MPI_SUM, jacobian.GetCommunicator());
#endif
		return sqrt(ret);
	}
	Residual & Residual::operator =(Residual const & other)
	{
		jacobian = other.jacobian;
		residual = other.residual;
		return *this;
	}
	void Residual::Rescale()
	{
#if defined(USE_OMP)
#pragma omp parallel for
#endif //USE_OMP
		for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k)
		{
			INMOST_DATA_REAL_TYPE norm = 0.0;
			for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
				norm += jacobian[k].GetValue(q)*jacobian[k].GetValue(q);
			norm = sqrt(norm);
			if( norm )
			{
				norm = 1.0/norm;
				residual[k] *= norm;
				for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
					jacobian[k].GetValue(q) *= norm;
			}
		}
	}
	Residual::Residual(std::string name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm)
	: jacobian(name,start,end,_comm),residual(name,start,end,_comm) 
	{
	}
	Residual::Residual(const Residual & other) 
	: jacobian(other.jacobian), residual(other.residual) 
	{
	}
#endif //USE_SOLVER
Kirill Terekhov's avatar
Kirill Terekhov committed
484 485 486
};

#endif //USE_AUTODIFF