autodiff.cpp 17.4 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
	Automatizator::Automatizator(const Automatizator & b) : name(b.name+"_copy")
	{
59
		std::vector<INMOST_DATA_ENUM_TYPE> regs = b.ListRegisteredEntries();
Kirill Terekhov's avatar
Kirill Terekhov committed
60
		for(std::vector<INMOST_DATA_ENUM_TYPE>::iterator kt = regs.begin(); kt != regs.end(); ++kt)
61 62
			RegisterEntry(b.GetEntry(*kt));
		if( b.last_num != 0 ) EnumerateEntries();
Kirill Terekhov's avatar
Kirill Terekhov committed
63 64 65 66 67 68
	}
	Automatizator & Automatizator::operator =(Automatizator const & b)
	{
		if( &b != this )
		{
			name = b.name+"_copy";
69 70 71 72 73 74
			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
Kirill Terekhov committed
75
			for(std::vector<INMOST_DATA_ENUM_TYPE>::iterator kt = regs.begin(); kt != regs.end(); ++kt)
76 77
				RegisterEntry(b.GetEntry(*kt));
			if( b.last_num != 0 ) EnumerateEntries();
Kirill Terekhov's avatar
Kirill Terekhov committed
78 79 80 81
		}
		return *this;
	}
	Automatizator::Automatizator(std::string _name) :name(_name), first_num(0), last_num(0) {}
Kirill Terekhov's avatar
Kirill Terekhov committed
82 83
	Automatizator::~Automatizator()
	{
84 85 86 87 88
		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
89
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
90
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterTag(Tag t, ElementType typemask, MarkerType domain_mask)
Kirill Terekhov's avatar
Kirill Terekhov committed
91
	{
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
		if( t.GetSize() == ENUMUNDEF )
			return RegisterEntry(VectorEntry(typemask,domain_mask,t));
		else if( t.GetSize() == 1 )
			return RegisterEntry(SingleEntry(typemask,domain_mask,t,0));
		else
		{
			BlockEntry b(typemask,domain_mask);
			for(int k = 0; k < t.GetSize(); ++k)
				b.AddTag(t,k);
			return RegisterEntry(b);
		}
	}
	INMOST_DATA_ENUM_TYPE Automatizator::RegisterEntry(const AbstractEntry & b)
	{
		Mesh * m = b.GetMeshLink();
		
		ElementType sparse = b.GetElementType();
		for (ElementType q = NODE; q <= MESH; q = NextElementType(q)) if (q & b.GetElementType())
Kirill Terekhov's avatar
Kirill Terekhov committed
110
		{
111 112
			for(unsigned unk = 0; unk < b.Size(); ++unk)
				sparse &= b.GetValueTag(unk).isSparse(q);
Kirill Terekhov's avatar
Kirill Terekhov committed
113
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
114
		
115 116
		INMOST_DATA_ENUM_TYPE ret = ENUMUNDEF;
		if( del_blocks.empty() )
Kirill Terekhov's avatar
Kirill Terekhov committed
117
		{
118 119 120
			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
Kirill Terekhov committed
121 122 123
		}
		else
		{
124 125 126 127 128 129 130 131 132 133 134 135 136
			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;
		
		{
			std::stringstream tag_name;
			tag_name << name << "_BLK_" << ret << "_Offset";
			reg_blocks[ret]->SetOffsetTag(m->CreateTag(tag_name.str(),DATA_INTEGER,b.GetElementType(),sparse,1));
Kirill Terekhov's avatar
Kirill Terekhov committed
137
		}
138
					
Kirill Terekhov's avatar
Kirill Terekhov committed
139 140
		return ret;
	}
141 142
	
	void Automatizator::UnregisterEntry(INMOST_DATA_ENUM_TYPE ind)
Kirill Terekhov's avatar
Kirill Terekhov committed
143
	{
144 145 146 147 148 149 150
		assert(act_blocks[ind]);
		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
Kirill Terekhov committed
151
	}
152
	void Automatizator::EnumerateEntries()
Kirill Terekhov's avatar
Kirill Terekhov committed
153 154 155
	{
		first_num = last_num = 0;
		const ElementType paralleltypes = NODE | EDGE | FACE | CELL;
156 157
		
		for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
Kirill Terekhov's avatar
Kirill Terekhov committed
158
		{
159 160 161
			AbstractEntry & b = *reg_blocks[it];
			TagInteger offset_tag = b.GetOffsetTag();
			Mesh * m = offset_tag.GetMeshLink();
Kirill Terekhov's avatar
Kirill Terekhov committed
162
			for (ElementType etype = NODE; etype <= MESH; etype = etype << 1)
163
				if (offset_tag.isDefined(etype) && offset_tag.isSparse(etype))
Kirill Terekhov's avatar
Kirill Terekhov committed
164
				{
165 166
					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
167 168
				}
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
169

170
		for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
Kirill Terekhov's avatar
Kirill Terekhov committed
171
		{
172 173 174 175
			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
176
			{
177
				for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) )
Kirill Terekhov's avatar
Kirill Terekhov committed
178
				{
179 180
					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
181
					{
182 183
						offset_tag[jt] = last_num;
						last_num += b.Size(jt);
Kirill Terekhov's avatar
Kirill Terekhov committed
184 185 186 187
					}
				}
			}
		}
188
		std::set<INMOST_DATA_ENUM_TYPE> Pre, Post; //Nonlocal indices
Kirill Terekhov's avatar
Kirill Terekhov committed
189
#if defined(USE_MPI)
Kirill Terekhov's avatar
Kirill Terekhov committed
190 191 192
		int size;
		MPI_Comm_size(MPI_COMM_WORLD,&size);
		if (size > 1)
Kirill Terekhov's avatar
Kirill Terekhov committed
193
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
194
			MPI_Scan(&last_num, &first_num, 1, INMOST_MPI_DATA_ENUM_TYPE, MPI_SUM, MPI_COMM_WORLD);
Kirill Terekhov's avatar
Kirill Terekhov committed
195
			first_num -= last_num;
196
			last_num += first_num;
Kirill Terekhov's avatar
Kirill Terekhov committed
197
			ElementType exch_mask = NONE;
198
			for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
Kirill Terekhov's avatar
Kirill Terekhov committed
199
			{
200 201 202 203
				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
204
				{
205 206
					exch_mask |= etype;
					for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) )
Kirill Terekhov's avatar
Kirill Terekhov committed
207
					{
208 209 210
						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
211 212 213
					}
				}
			}
214
			//synchronize indices
Kirill Terekhov's avatar
Kirill Terekhov committed
215
			{
Kirill Terekhov's avatar
Kirill Terekhov committed
216
				std::map<Mesh *,std::vector<Tag> > exch_tags;
217 218
				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
219 220
				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
221
			}
222
			//compute out-of-bounds indices
223
			for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
224
			{
225 226 227 228
				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 )
229
				{
230
					for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) )
231
					{
232 233
						Element jt = m->ElementByLocalID(etype,kt);
						if ((!(etype & paralleltypes) || (jt.GetStatus() == Element::Ghost)) && b.isValid(jt) && b.Size(jt))
234
						{
Kirill Terekhov's avatar
Kirill Terekhov committed
235 236 237 238 239 240 241 242 243
							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);
								}
							}
244
						}
245
					}
246 247 248
				} //etype
			} //it
			// after cycle
Kirill Terekhov's avatar
Kirill Terekhov committed
249 250
		}
#endif
251 252
		//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
253 254
#if defined(USE_OMP)
#pragma omp parallel
255
#endif //USE_OMP
256
		{
257
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
258
#pragma omp single
259
#endif //USE_OMP
260
			{
261
				merger.resize(MAX_THREADS);
262
			}
263
			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);
264
		}
Kirill Terekhov's avatar
Kirill Terekhov committed
265
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
266
	
267
	std::vector<INMOST_DATA_ENUM_TYPE> Automatizator::ListRegisteredEntries() const
Kirill Terekhov's avatar
Kirill Terekhov committed
268 269
	{
		std::vector<INMOST_DATA_ENUM_TYPE> ret;
270
		for(blk_enum::size_type it = 0; it < reg_blocks.size(); ++it) if( isRegisteredEntry(it) )
Kirill Terekhov's avatar
Kirill Terekhov committed
271 272 273 274
			ret.push_back(static_cast<INMOST_DATA_ENUM_TYPE>(it));
		return ret;
	}
	
275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 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 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 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 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421
	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
	{
		unsigned pos = 0, k = 0;
		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)
	{
		unsigned pos = 0, k = 0;
		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
	{
		unsigned pos = 0, k = 0;
		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
	{
		unsigned pos = 0, k = 0;
		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;
	}
	
	rMatrix MultiEntry::Value(const Storage & e) const
	{
		vMatrix ret(MatrixSize(e),1);
		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;
	}
	
	iMatrix MultiEntry::Index(const Storage & e) const
	{
		iMatrix ret(MatrixSize(e),1);
		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;
	}
	
	uMatrix MultiEntry::operator [](const Storage & e) const
	{
		uMatrix ret(MatrixSize(e),1);
		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;
		while( pos + entries[k]->Size() < unk ) pos += entries[k++]->Size();
		assert(k < entries.size());
		return entries[k]->GetValueComp(unk-pos);
	}
	
	TagRealArray MultiEntry::GetValueTag(INMOST_DATA_ENUM_TYPE unk) const
	{
		unsigned pos = 0, k = 0;
		while( pos + entries[k]->Size() < unk ) pos += entries[k++]->Size();
		assert(k < entries.size());
		return entries[k]->GetValueTag(unk-pos);
	}
	
	AbstractEntry * MultiEntry::Copy() const
	{
		MultiEntry * ret = new MultiEntry(GetElementType(),GetMask());
		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
Kirill Terekhov committed
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
		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) 
			if( status_tbl[stat][k] ) ret++;
		return ret;
	}
	
	AbstractEntry * StatusBlockEntry::Copy() const 
	{
		StatusBlockEntry * ret = new StatusBlockEntry(GetElementType(),GetMask(),status_tag,status_tbl); 
		for(unsigned k = 0; k < Size(); ++k) 
			ret->AddTag(unknown_tags[k],unknown_comp[k]); 
		return ret; 
Kirill Terekhov's avatar
Kirill Terekhov committed
466
	}
Kirill Terekhov's avatar
Kirill Terekhov committed
467
#endif //USE_MESH
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487

#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();
	}
488 489 490 491
	void Residual::ClearHessian()
	{
		for(Sparse::HessianMatrix::iterator it = hessian.Begin(); it != hessian.End(); ++it) it->Clear();
	}
492 493 494 495 496 497 498 499
	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;
500 501
			if( !jacobian.Empty() ) jacobian[k].Clear();
			if( !hessian.Empty() ) hessian[k].Clear();
502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519
		}
	}
	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)
	{
520
		hessian = other.hessian;
521 522 523 524
		jacobian = other.jacobian;
		residual = other.residual;
		return *this;
	}
525
	void Residual::Rescale(INMOST_DATA_ENUM_TYPE p)
526 527 528 529 530 531 532
	{
#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;
533 534 535 536 537 538 539 540 541 542 543 544
			if( p == ENUMUNDEF ) //infinite norm
			{
				for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
					if( norm < fabs(jacobian[k].GetValue(q)) )
						norm = fabs(jacobian[k].GetValue(q));
			}
			else //p-norm
			{
				for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q)
					norm += pow(fabs(jacobian[k].GetValue(q)),p);
				norm = pow(norm,1.0/p);
			}
545 546 547 548 549 550
			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;
551 552 553
				if( !hessian.Empty() )
					for(INMOST_DATA_ENUM_TYPE q = 0; q < hessian[k].Size(); ++q)
						hessian[k].GetValue(q) *= norm;
554 555 556 557
			}
		}
	}
	Residual::Residual(std::string name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm)
558
	: jacobian(name,start,end,_comm),residual(name,start,end,_comm), hessian(name,0,0,_comm)
559 560 561
	{
	}
	Residual::Residual(const Residual & other) 
562
	: jacobian(other.jacobian), residual(other.residual), hessian(other.hessian)
563 564
	{
	}
565 566 567 568 569 570 571 572 573

	Matrix<multivar_expression_reference> Residual::operator [](const AbstractMatrix<INMOST_DATA_INTEGER_TYPE> & rows)	
	{
		Matrix<multivar_expression_reference> ret(rows.Rows(),rows.Cols());
		for(INMOST_DATA_ENUM_TYPE i = 0; i < rows.Rows(); ++i)
			for(INMOST_DATA_ENUM_TYPE j = 0; j < rows.Cols(); ++j)
				ret(i,j) = multivar_expression_reference(residual[rows(i,j)],&jacobian[rows(i,j)]);
		return ret;
	}
574
#endif //USE_SOLVER
Kirill Terekhov's avatar
Kirill Terekhov committed
575 576 577
};

#endif //USE_AUTODIFF