inmost_autodiff.h 6.62 KB
Newer Older
Kirill Terekhov's avatar
Fixes  
Kirill Terekhov committed
1

Kirill Terekhov's avatar
Kirill Terekhov committed
2 3 4 5
#ifndef INMOST_AUTODIFF_H_INCLUDED
#define INMOST_AUTODIFF_H_INCLUDED
#include "inmost_common.h"
#include "inmost_mesh.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
6
#include "inmost_solver.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
7
#include "inmost_variable.h"
Kirill Terekhov's avatar
Kirill Terekhov committed
8 9
#include <sstream> //for debug

Kirill Terekhov's avatar
Kirill Terekhov committed
10
//#define NEW_VERSION
Kirill Terekhov's avatar
Kirill Terekhov committed
11

Kirill Terekhov's avatar
Kirill Terekhov committed
12 13 14 15
//#if defined(USE_AUTODIFF) && (!defined(USE_MESH))
//#warning "USE_AUTODIFF require USE_MESH"
//#undef USE_AUTODIFF
//#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
16 17 18 19 20 21 22 23 24 25 26 27 28

//#define DPRNT



#if defined(USE_AUTODIFF)
#include <math.h>



namespace INMOST
{
	class Automatizator; //forward declaration
29

Kirill Terekhov's avatar
Kirill Terekhov committed
30
#if defined(USE_SOLVER)
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
	class Residual
	{
		Sparse::Matrix jacobian;
		Sparse::Vector residual;
		Sparse::LockService locks;
	public:
		Residual(std::string name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD)
			: jacobian(name,start,end,_comm),residual(name,start,end,_comm) {}
		Residual(const Residual & other)
			: jacobian(other.jacobian), residual(other.residual)
		{}
		Residual & operator =(Residual const & other)
		{
			jacobian = other.jacobian;
			residual = other.residual;
			return *this;
		}
		Sparse::Matrix & GetJacobian() {return jacobian;}
		const Sparse::Matrix & GetJacobian() const {return jacobian;}
		Sparse::Vector & GetResidual() {return residual;}
		const Sparse::Vector & GetResidual() const {return residual;}
		INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();}
		INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();}
		void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const 
		{
			start = residual.GetFirstIndex(); 
			end = residual.GetLastIndex();
		}
		void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end)
		{
			jacobian.SetInterval(beg,end);
			residual.SetInterval(beg,end);
		}
		multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row)
		{
			return multivar_expression_reference(residual[row],&jacobian[row]);
		}
		void ClearResidual()
		{
			for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0;
		}
		void ClearJacobian()
		{
			for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it)
				it->Clear();
		}
		void Clear()
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
79 80 81
#if defined(USE_OMP)
#pragma omp for
#endif
82 83 84 85 86 87 88 89 90
			for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) 
			{
				residual[k] = 0.0;
				jacobian[k].Clear();
			}
		}
		INMOST_DATA_REAL_TYPE Norm()
		{
			INMOST_DATA_REAL_TYPE ret = 0;
Kirill Terekhov's avatar
Kirill Terekhov committed
91
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
92
#pragma omp parallel for reduction(+:ret)
Kirill Terekhov's avatar
Kirill Terekhov committed
93
#endif
94 95
			for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) 
				ret += residual[k]*residual[k];
Kirill Terekhov's avatar
Kirill Terekhov committed
96
#if defined(USE_MPI)
97 98
			INMOST_DATA_REAL_TYPE tmp = ret;
			MPI_Allreduce(&tmp, &ret, 1, INMOST_MPI_DATA_REAL_TYPE, MPI_SUM, jacobian.GetCommunicator());
Kirill Terekhov's avatar
Kirill Terekhov committed
99
#endif
100 101 102 103 104
			return sqrt(ret);
		}
		/// Normalize entries in jacobian and right hand side
		void Rescale()
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
105
#if defined(USE_OMP)
Kirill Terekhov's avatar
Kirill Terekhov committed
106
#pragma omp parallel for
Kirill Terekhov's avatar
Kirill Terekhov committed
107
#endif
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
			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;
				}
			}
		}
		void InitLocks() {locks.SetInterval(GetFirstIndex(),GetLastIndex());}
		void Lock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.Lock(pos);}
		void UnLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.UnLock(pos);}
		void TestLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.TestLock(pos);}
	};
Kirill Terekhov's avatar
Kirill Terekhov committed
128 129 130
#endif
	
#if defined(USE_MESH)
Kirill Terekhov's avatar
Kirill Terekhov committed
131 132 133
	class Automatizator
	{
	private:
134
		static Automatizator * CurrentAutomatizator;
Kirill Terekhov's avatar
Kirill Terekhov committed
135
#if defined(USE_OMP)
136
		std::vector<Sparse::RowMerger> merger;
Kirill Terekhov's avatar
Kirill Terekhov committed
137
#else
138
		Sparse::RowMerger merger;
Kirill Terekhov's avatar
Kirill Terekhov committed
139
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
140
		typedef struct{ Tag t; MarkerType domain_mask; } tagdomain;
Kirill Terekhov's avatar
Kirill Terekhov committed
141
		typedef struct{ tagdomain d; Tag indices; } tagpair;
Kirill Terekhov's avatar
Kirill Terekhov committed
142
		typedef dynarray<tagpair, 128> tagpairs_type;
Kirill Terekhov's avatar
Kirill Terekhov committed
143 144
		typedef std::vector<tagpair> index_enum;
	private:
Kirill Terekhov's avatar
Kirill Terekhov committed
145 146
		index_enum            index_tags;
		tagpairs_type         reg_tags;
Kirill Terekhov's avatar
Kirill Terekhov committed
147 148 149
		INMOST_DATA_ENUM_TYPE first_num;
		INMOST_DATA_ENUM_TYPE last_num;
	public:
Kirill Terekhov's avatar
Kirill Terekhov committed
150
		Automatizator();
Kirill Terekhov's avatar
Kirill Terekhov committed
151
		~Automatizator();
Kirill Terekhov's avatar
Kirill Terekhov committed
152 153
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetFirstIndex() { return first_num; }
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetLastIndex() { return last_num; }
154 155 156 157 158
		/// Set data of tag t defined on domain_mask to be dynamic data.
		/// \warning
		/// Don't register tag twice.
		/// \todo
		/// Read comments inside, change merger.Resize() behavior.
Kirill Terekhov's avatar
Kirill Terekhov committed
159
		INMOST_DATA_ENUM_TYPE                                            RegisterTag(Tag t, ElementType typemask, MarkerType domain_mask = 0);
160
		/// Set index for every data entry of dynamic tag.
Kirill Terekhov's avatar
Kirill Terekhov committed
161 162 163 164 165 166 167 168
		void                                                             EnumerateTags();
		__INLINE Tag                                                     GetValueTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind].d.t; }
		__INLINE Tag                                                     GetIndexTag(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind].indices; }
		__INLINE MarkerType                                              GetMask(INMOST_DATA_ENUM_TYPE ind) { return reg_tags[ind].d.domain_mask; }
		__INLINE INMOST_DATA_REAL_TYPE                                   GetValue(const Storage & e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->RealArray(GetValueTag(ind))[comp]; }
		__INLINE INMOST_DATA_ENUM_TYPE                                   GetIndex(const Storage & e, INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_ENUM_TYPE comp = 0) { return e->IntegerArray(GetIndexTag(ind))[comp]; }
		__INLINE bool                                                    isValid(const Storage & e, INMOST_DATA_ENUM_TYPE ind) { MarkerType mask = GetMask(ind); return mask == 0 || e->GetMarker(mask); }
		Sparse::RowMerger &                                              GetMerger()
169
		{
Kirill Terekhov's avatar
Kirill Terekhov committed
170
#if defined(USE_OMP)
171
			return merger[omp_get_thread_num()];
Kirill Terekhov's avatar
Kirill Terekhov committed
172
#else
173
			return merger;
Kirill Terekhov's avatar
Kirill Terekhov committed
174
#endif
Kirill Terekhov's avatar
Kirill Terekhov committed
175
		}
176 177 178 179 180 181 182 183
		/// Remove global current automatizator.
		static void RemoveCurrent() {CurrentAutomatizator = NULL;}
		/// Set current global automatizator, so that variable will be optimized with row merger.
		static void MakeCurrent(Automatizator * aut) {CurrentAutomatizator = aut;}
		/// Check that there is an automatizator.
		static bool HaveCurrent() {return CurrentAutomatizator != NULL;}
		/// Retrive the automatizator.
		static Automatizator * GetCurrent() {return CurrentAutomatizator;}
Kirill Terekhov's avatar
Kirill Terekhov committed
184
	};
Kirill Terekhov's avatar
Kirill Terekhov committed
185
#endif
186
} //namespace INMOST
Kirill Terekhov's avatar
Kirill Terekhov committed
187

188 189
#endif //USE_AUTODIFF
#endif //INMOST_AUTODIFF_H_INCLUDED