WSL/SLF GitLab Repository

arrays.cc 5.24 KB
Newer Older
1
2
3
4
5
6
#include <time.h>
#include <meteoio/MeteoIO.h>

using namespace std;
using namespace mio;

7
8
bool array1d(const unsigned int& n) {
	cout << "Testing Array1D<double>\n";
9
10
11
12
	bool status = true;
	srand((unsigned)time(0));
	const double range = 10.;

13
	Array1D<double> grid1(n);
14

15
16
	for(unsigned int ii=0; ii<grid1.getNx(); ii++) {
		grid1(ii) = (double)rand()/(double)RAND_MAX*range;
17
18
	}

19
20
21
22
23
24
	Array1D<double> grid2=grid1;
	grid2 -= 8.;
	grid2 /= 2.;
	grid2 += 4.;
	grid2 *= 2.;
	if(!grid2.checkEpsilonEquality(grid1, 1e-6)) {
25
26
27
28
		cout << "\terror: basic operations with constants fail!\n";
		status=false;
	}

29
30
31
32
	grid2.resize(n, 0.);
	grid2 += grid1;
	grid2 /= grid1;
	if(!IOUtils::checkEpsilonEquality(grid2.getMean(), 1., 1e-6)) {
33
34
35
		cout << "\terror: grids addition or division fails!\n";
		status=false;
	}
36
37
	grid2 *= grid1;
	if(!grid2.checkEpsilonEquality(grid1, 1e-6)) {
38
39
40
		cout << "\terror: grids multiplication fails!\n";
		status=false;
	}
41
42
	grid2 -= grid1;
	if(!IOUtils::checkEpsilonEquality(grid2.getMean(), 0., 1e-6)) {
43
44
45
46
47
48
49
		cout << "\terror: grids substraction fails!\n";
		status=false;
	}

	return status;
}

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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
bool grid2d(const unsigned int& n) {
	cout << "Testing Grid2DObject\n";
	bool status = true;
	srand((unsigned)time(0));
	const double range = 10.;
	Coords llcorner("CH1903","");
	llcorner.setXY(785425. , 191124., 1400.);

	Grid2DObject grid1(n, n, 100, llcorner);

	for(unsigned int jj=0; jj<grid1.getNy(); jj++) {
		for(unsigned int ii=0; ii<grid1.getNx(); ii++) {
			grid1.grid2D(ii,jj) = (double)rand()/(double)RAND_MAX*range;
		}
	}

	Grid2DObject grid2(n, n, 100, llcorner, grid1.grid2D);
	grid2.grid2D -= 8.;
	grid2.grid2D /= 2.;
	grid2.grid2D += 4.;
	grid2.grid2D *= 2.;
	if(!grid2.grid2D.checkEpsilonEquality(grid1.grid2D, 1e-6)) {
		cout << "\terror: basic operations with constants fail!\n";
		status=false;
	}

	grid2.set(n, n, 100, llcorner, 0.);
	grid2.grid2D += grid1.grid2D;
	grid2.grid2D /= grid1.grid2D;
	if(!IOUtils::checkEpsilonEquality(grid2.grid2D.getMean(), 1., 1e-6)) {
		cout << "\terror: grids addition or division fails!\n";
		status=false;
	}
	grid2.grid2D *= grid1.grid2D;
	if(!grid2.grid2D.checkEpsilonEquality(grid1.grid2D, 1e-6)) {
		cout << "\terror: grids multiplication fails!\n";
		status=false;
	}
	grid2.grid2D -= grid1.grid2D;
	if(!IOUtils::checkEpsilonEquality(grid2.grid2D.getMean(), 0., 1e-6)) {
		cout << "\terror: grids substraction fails!\n";
		status=false;
	}

	return status;
}

97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
bool grid3d(const unsigned int& n) {
	cout << "Testing Grid3DObject\n";
	bool status = true;
	srand((unsigned)time(0));
	const double range = 10.;
	Coords llcorner("CH1903","");
	llcorner.setXY(785425. , 191124., 1400.);

	Grid3DObject grid1(n, n, n, 100, llcorner);

	for(unsigned int kk=0; kk<grid1.getNz(); kk++) {
		for(unsigned int jj=0; jj<grid1.getNy(); jj++) {
			for(unsigned int ii=0; ii<grid1.getNx(); ii++) {
				grid1.grid3D(ii,jj,kk) = (double)rand()/(double)RAND_MAX*range;
			}
		}
	}

	Grid3DObject grid2(n, n, n, 100, llcorner, grid1.grid3D);
	grid2.grid3D -= 8.;
	grid2.grid3D /= 2.;
	grid2.grid3D += 4.;
	grid2.grid3D *= 2.;
	if(!grid2.grid3D.checkEpsilonEquality(grid1.grid3D, 1e-6)) {
		cout << "\terror: basic operations with constants fail!\n";
		status=false;
	}

	grid2.set(n, n, n, 100, llcorner, 0.);
	grid2.grid3D += grid1.grid3D;
	grid2.grid3D /= grid1.grid3D;
	if(!IOUtils::checkEpsilonEquality(grid2.grid3D.getMean(), 1., 1e-6)) {
		cout << "\terror: grids addition or division fails!\n";
		status=false;
	}
	grid2.grid3D *= grid1.grid3D;
	if(!grid2.grid3D.checkEpsilonEquality(grid1.grid3D, 1e-6)) {
		cout << "\terror: grids multiplication fails!\n";
		status=false;
	}
	grid2.grid3D -= grid1.grid3D;
	if(!IOUtils::checkEpsilonEquality(grid2.grid3D.getMean(), 0., 1e-6)) {
		cout << "\terror: grids substraction fails!\n";
		status=false;
	}

	return status;
}

146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
bool matrix(const unsigned int& n) {
	cout << "Testing Matrix\n";
	bool status=true;
	Matrix I(n,1.); //build an n*n identity matrix
	Matrix m1(n,n);
	m1.random(10.);

	const double det = m1.det();
	const Matrix m1_trans = m1.getT();
	const double det_trans = m1_trans.det();

	if(!IOUtils::checkEpsilonEquality(det, det_trans, Matrix::epsilon_mtr*fabs(det))) {
		cout << "\terror: m1.det != m1T.det\n";
		status=false;
	}

	Matrix m2 = m1-8.;
	m2 /= 2.;
	m2 += 4.;
	m2 *= 2.;
	if(m2!=m1) {
		cout << "\terror: basic operations with constants fail!\n";
		status=false;
	}

171
172
173
174
175
176
177
178
	m2 = m1+m1;
	m2 -= m1;
	if(m2 != m1) {
		cout << "\terror when adding/substracting matrix\n";
		status=false;
	}


179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
	m2 = m1.getInv();
	Matrix m3=m1*m2;
	if(m3.isIdentity()!=true) {
		cout << "\terror: m1*inv(m1) is NOT identity matrix\n";
		status=false;
	}

	m3 = m1*m1;
	if(m2*m3 != m1) {
		cout << "\terror when multiplying matrix\n";
		status=false;
	}

	Matrix m4=Matrix::solve(m1,I); //solve m1*X=I
	if(m1*m4 != I){
		cout << "\terror when solving A*X=B\n";
		status=false;
	}
197

198
199
200
201
202
203
204
205
206
207
208
209
210
	Matrix L,U;
	if(m1.LU(L,U)==false) {
		cout << "\terror: LU decomposition could NOT be computed\n";
		status=false;
	}

	return status;
}


int main() {
	const unsigned int n=50;

211
	const bool grid1d_status = array1d(n);
212
	const bool grid2d_status = grid2d(n);
213
	const bool grid3d_status = grid3d(n);
214
	const bool matrix_status = matrix(n);
215
	if(grid1d_status!=true || grid2d_status!=true || grid3d_status!=true || matrix_status!=true) throw IOException("Grid/Matrix error", AT);
216
217
	return 0;
}