checks.cpp 1.86 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
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
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
79
80
#include "checks.h"

using namespace INMOST;


//shortcuts
typedef Storage::bulk       bulk;
typedef Storage::enumerator enumerator;
typedef Storage::real       real;

bool check_flux_properties(int i, const variable & flux, const std::string & name, bool print)
{
	bool test = true;
	const Sparse::Row & r = flux.GetRow();
	for(int k = 0; k < (int)r.Size(); ++k)
	{
		if( r.GetIndex(k) == i ) //diagonal element
		{
			if( r.GetValue(k) < 0.0 )
			{
				if( print ) std::cout << "flux " << name << " gives negative diagonal contribution to " << i << "-th variable" << std::endl;
				test = false;
			}
		}
		else //off-diagonal element
		{
			if( r.GetValue(k) > 0.0 )
			{
				if( print ) std::cout << "flux " << name << " gives positive off-diagonal contribution to " << i << "-th variable" << std::endl;
				test = false;
			}
		}
	}
	if( !test && print ) 
	{
		r.Print();
		//scanf("%*c");
	}
	return test;
}


bool check_matrix_properties(const Sparse::Matrix & A, bool print)
{
	bool is_m_matrix = true;
	enumerator mbeg,mend;
	A.GetInterval(mbeg,mend);
	for(enumerator i = mbeg; i < mend; ++i)
	{
		const Sparse::Row & r = A[i];
		real rowsum = 0.0;
		for(enumerator k = 0; k < r.Size(); ++k)
		{
			rowsum += r.GetValue(k);
			if( r.GetIndex(k) != i ) 
			{
				if( r.GetValue(k) > 0.0 )
				{
					if( print ) std::cout << "Off-diagonal element (" << i << "," << r.GetIndex(k) << ") is positive: " << r.GetValue(k) << std::endl;
					is_m_matrix = false;
				}
			}
			else 
			{
				if( r.GetValue(k) < 0.0 )
				{
					if( print ) std::cout << "Diagonal element (" << i << "," << r.GetIndex(k) << ") is negative: " << r.GetValue(k) << std::endl;
					is_m_matrix = false;
				}
			}
		}
		if( rowsum < -1.0e-9 )
		{
			std::cout << "Row-sum is " << rowsum << std::endl;
			is_m_matrix = false;
		}
	}
	if( !is_m_matrix ) std::cout << "Not an M-matrix" << std::endl;
	return is_m_matrix;
}