Line data Source code
1 : #include <gmock/gmock.h>
2 : #include <gtest/gtest.h>
3 : #include <cassert>
4 : #include <cfloat>
5 :
6 : extern "C"
7 : {
8 : #include "dist.h"
9 : #include "qrfactorization.h"
10 : }
11 :
12 : #define DOUBLE_NEAR(x) DoubleNear((x), 4.0 * DBL_EPSILON)
13 : #define DOUBLE_NEAR_MULTIPLIER(x, multiplier) \
14 : DoubleNear((x), (multiplier)*DBL_EPSILON)
15 : #define map(i, j) ((j)-1) * lda + ((i)-1)
16 : #define a(i, j) a[ map(i, j) ]
17 : #define q(i, j) q[ map(i, j) ]
18 : #define r(i, j) r[ map(i, j) ]
19 : #define map1(i) (i) - 1
20 : #define x(i) x[ map1(i) ]
21 : #define y(i) y[ map1(i) ]
22 : #define ipiv(i) ipiv[ map1(i) ]
23 : #define wrk(i) wrk[ map1(i) ]
24 : #define tau(i) tau[ map1(i) ]
25 :
26 : using namespace std;
27 : using namespace testing;
28 :
29 : class TestThatHouseholderG : public Test
30 : {
31 : public:
32 : int n = 10, incx = 1;
33 : double alpha = 2.0;
34 : double tau = 2.0;
35 : double *xvector = NULL;
36 : double epmach = DBL_EPSILON;
37 : double safmin = DBL_MIN;
38 : const double cutlo = 4.44089e-16;
39 : const double cuthi = 1.30438e19;
40 :
41 10 : void SetUp() override
42 : {
43 10 : xvector = (double *)calloc((size_t)n, (size_t)sizeof(double));
44 : // if (xvector == NULL)
45 : //{
46 : // printf("Failed to allocate xvector memory!\n");
47 : // assert(xvector);
48 : //}
49 10 : }
50 :
51 10 : void TearDown() override { free(xvector); }
52 : };
53 :
54 2 : TEST_F(TestThatHouseholderG, CanAccessEPMACHandSAFMIN)
55 : {
56 1 : EXPECT_THAT(0.2220446049250313E-15, DOUBLE_NEAR(epmach));
57 1 : EXPECT_THAT(0.2225073858507201E-307, DOUBLE_NEAR(safmin));
58 1 : }
59 :
60 2 : TEST_F(TestThatHouseholderG, CanAccessTheHouseholderReflectorMethod)
61 : {
62 1 : ASSERT_THAT(housg(n, &alpha, xvector, incx, &tau, safmin), Eq(0));
63 : }
64 :
65 2 : TEST_F(TestThatHouseholderG, IsDefinedForINCXEqualsZeroButFlagsNegativeTwo)
66 : {
67 1 : int myIncx = 0;
68 1 : EXPECT_THAT(housg(n, &alpha, xvector, myIncx, &tau, safmin), Eq(-2));
69 1 : EXPECT_THAT(tau, DoubleEq(0.0));
70 1 : }
71 :
72 2 : TEST_F(TestThatHouseholderG, IsDefinedForNegativeDimensionsButFlagsNegativeOne)
73 : {
74 1 : int dim = -10;
75 1 : EXPECT_THAT(housg(dim, &alpha, xvector, incx, &tau, safmin), Eq(-1));
76 1 : EXPECT_THAT(tau, DoubleEq(0.0));
77 1 : }
78 :
79 2 : TEST_F(TestThatHouseholderG, IsIdentityForOneDimensionalApplications)
80 : {
81 1 : int dim = 1;
82 1 : EXPECT_THAT(housg(dim, &alpha, xvector, incx, &tau, safmin), Eq(0));
83 1 : EXPECT_THAT(tau, DoubleEq(0.0));
84 1 : }
85 :
86 2 : TEST_F(TestThatHouseholderG, IsIdentityForXEqualsZero)
87 : {
88 1 : EXPECT_THAT(housg(n, &alpha, xvector, incx, &tau, safmin), Eq(0));
89 1 : EXPECT_THAT(tau, DoubleEq(0.0));
90 1 : }
91 :
92 2 : TEST_F(TestThatHouseholderG, ProductOfHTransposeHIsIdentityAndVerifyNorm)
93 : {
94 1 : xvector[ 0 ] = 1.0;
95 1 : xvector[ 1 ] = 1.0;
96 1 : xvector[ 2 ] = 1.0;
97 1 : xvector[ 3 ] = 1.0;
98 1 : xvector[ 4 ] = 1.0;
99 1 : xvector[ 5 ] = 1.0;
100 1 : xvector[ 6 ] = 1.0;
101 1 : xvector[ 7 ] = 1.0;
102 1 : xvector[ 8 ] = 1.0;
103 1 : EXPECT_THAT(ddist2(n - 1, xvector, incx, xvector, incx, 1), DOUBLE_NEAR(3.0));
104 :
105 1 : EXPECT_THAT(housg(n, &alpha, xvector, incx, &tau, safmin), Eq(0));
106 1 : double vnorm = ddist2(n - 1, xvector, incx, xvector, incx, 1);
107 1 : EXPECT_THAT(-2.0 * tau + tau * tau * (1.0 + vnorm * vnorm), DOUBLE_NEAR(0.0));
108 1 : }
109 :
110 2 : TEST_F(TestThatHouseholderG,
111 : ProductOfHTransposeHIsIdentityAndVerifyNormNoUnderflow)
112 : {
113 1 : xvector[ 0 ] = cutlo / 9;
114 1 : xvector[ 1 ] = cutlo / 8;
115 1 : xvector[ 2 ] = cutlo / 7;
116 1 : xvector[ 3 ] = cutlo / 6;
117 1 : xvector[ 4 ] = cutlo / 5;
118 1 : xvector[ 5 ] = cutlo / 4;
119 1 : xvector[ 6 ] = cutlo / 3;
120 1 : xvector[ 7 ] = cutlo / 2;
121 1 : xvector[ 8 ] = cutlo;
122 1 : EXPECT_THAT(ddist2(n - 1, xvector, incx, xvector, incx, 1),
123 : DOUBLE_NEAR(5.510583948830441e-16));
124 :
125 1 : EXPECT_THAT(housg(n, &alpha, xvector, incx, &tau, safmin), Eq(0));
126 1 : double vnorm = ddist2(n - 1, xvector, incx, xvector, incx, 1);
127 1 : EXPECT_THAT(-2.0 * tau + tau * tau * (1.0 + vnorm * vnorm), DOUBLE_NEAR(0.0));
128 1 : }
129 :
130 2 : TEST_F(TestThatHouseholderG,
131 : ProductOfHTransposeHIsIdentityAndVerifyNormNoOverflow)
132 : {
133 1 : xvector[ 0 ] = cuthi * 2;
134 1 : xvector[ 1 ] = cuthi * 3;
135 1 : xvector[ 2 ] = cuthi * 4;
136 1 : xvector[ 3 ] = cuthi * 5;
137 1 : xvector[ 4 ] = cuthi * 6;
138 1 : xvector[ 5 ] = cuthi * 7;
139 1 : xvector[ 6 ] = cuthi * 8;
140 1 : xvector[ 7 ] = cuthi * 9;
141 1 : xvector[ 8 ] = cuthi * 10;
142 1 : EXPECT_THAT(ddist2(n - 1, xvector, incx, xvector, incx, 1),
143 : DOUBLE_NEAR(2.556052344553217e+20 + 32768));
144 :
145 1 : EXPECT_THAT(housg(n, &alpha, xvector, incx, &tau, safmin), Eq(0));
146 1 : double vnorm = ddist2(n - 1, xvector, incx, xvector, incx, 1);
147 1 : EXPECT_THAT(-2.0 * tau + tau * tau * (1.0 + vnorm * vnorm), DOUBLE_NEAR(0.0));
148 1 : }
149 :
150 2 : TEST_F(TestThatHouseholderG,
151 : ProductOfHTransposeHIsIdentityAndVerifyNormNoUnderflowNoOverflow)
152 : {
153 1 : xvector[ 0 ] = 1.0;
154 1 : xvector[ 1 ] = 1.0;
155 1 : xvector[ 2 ] = cutlo;
156 1 : xvector[ 3 ] = 1.0;
157 1 : xvector[ 4 ] = 1.0;
158 1 : xvector[ 5 ] = cuthi;
159 1 : xvector[ 6 ] = 1.0;
160 1 : xvector[ 7 ] = cutlo;
161 1 : xvector[ 8 ] = 1.0;
162 1 : EXPECT_THAT(ddist2(n - 1, xvector, incx, xvector, incx, 1),
163 : DOUBLE_NEAR(1.304380000000000e+19));
164 :
165 1 : EXPECT_THAT(housg(n, &alpha, xvector, incx, &tau, safmin), Eq(0));
166 1 : double vnorm = ddist2(n - 1, xvector, incx, xvector, incx, 1);
167 1 : EXPECT_THAT(-2.0 * tau + tau * tau * (1.0 + vnorm * vnorm), DOUBLE_NEAR(0.0));
168 1 : }
169 :
170 : class TestThatHouseholderL : public Test
171 : {
172 : public:
173 : int m = 4, lda = 4, n = 3, ier = 0, incx = 1;
174 : double tau = 2.0;
175 : double *a = NULL;
176 : double *r = NULL;
177 : double *xvector = NULL;
178 : double epmach = DBL_EPSILON;
179 : double safmin = DBL_MIN;
180 :
181 10 : void SetUp() override
182 : {
183 10 : xvector = (double *)calloc((size_t)m, (size_t)sizeof(double));
184 : // if (xvector == NULL)
185 : //{
186 : // printf("Failed to allocate xvector memory!\n");
187 : // assert(xvector);
188 : //}
189 :
190 10 : a = (double *)calloc((size_t)(lda * n), (size_t)sizeof(double));
191 : // if (a == NULL)
192 : //{
193 : // printf("Failed to allocate matrix A memory!\n");
194 : // assert(a);
195 : //}
196 :
197 10 : a(1, 1) = 1.0, a(1, 2) = 2.0, a(1, 3) = 3.0;
198 10 : a(2, 1) = 1.0, a(2, 2) = 5.0, a(2, 3) = 6.0;
199 10 : a(3, 1) = 1.0, a(3, 2) = 8.0, a(3, 3) = 9.0;
200 10 : a(4, 1) = 1.0, a(4, 2) = 11.0, a(4, 3) = 12.0;
201 :
202 10 : r = (double *)calloc((size_t)(lda * n), (size_t)sizeof(double));
203 : // if (r == NULL)
204 : //{
205 : // printf("Failed to allocate space for result r!\n");
206 : // assert(r);
207 : //}
208 10 : }
209 :
210 10 : void TearDown() override
211 : {
212 10 : free(r);
213 10 : free(a);
214 10 : free(xvector);
215 10 : }
216 : };
217 :
218 2 : TEST_F(TestThatHouseholderL, CanAccessTheHouseholderLeftMultiplyMethod)
219 : {
220 1 : ASSERT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(0));
221 : }
222 :
223 2 : TEST_F(TestThatHouseholderL, IsNotDefinedForNegativeRowDimensionM)
224 : {
225 1 : m = -1;
226 1 : EXPECT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(1));
227 1 : EXPECT_THAT(ier, Eq(1));
228 1 : }
229 :
230 2 : TEST_F(TestThatHouseholderL, IsNotDefinedForNegativeColumnDimensionN)
231 : {
232 1 : n = -1;
233 1 : EXPECT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(1));
234 1 : EXPECT_THAT(ier, Eq(1));
235 1 : }
236 :
237 2 : TEST_F(TestThatHouseholderL, IsNotDefinedForINCXEqualsZero)
238 : {
239 1 : incx = 0;
240 1 : EXPECT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(1));
241 1 : EXPECT_THAT(ier, Eq(1));
242 1 : }
243 :
244 2 : TEST_F(TestThatHouseholderL,
245 : IsNotDefinedForLeadingDimensionLDALessThanNumberOfRowsM)
246 : {
247 1 : lda = 1;
248 1 : EXPECT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(1));
249 1 : EXPECT_THAT(ier, Eq(1));
250 1 : }
251 :
252 2 : TEST_F(TestThatHouseholderL, IsDefinedForRowDimensionMEqualsZero)
253 : {
254 1 : m = 0;
255 1 : EXPECT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(0));
256 1 : EXPECT_THAT(ier, Eq(0));
257 1 : }
258 :
259 2 : TEST_F(TestThatHouseholderL, IsDefinedForColumnDimensionNEqualsZero)
260 : {
261 1 : n = 0;
262 1 : EXPECT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(0));
263 1 : EXPECT_THAT(ier, Eq(0));
264 1 : }
265 :
266 2 : TEST_F(TestThatHouseholderL, IsDefinedForTAUEqualsZero)
267 : {
268 1 : tau = 0.0;
269 1 : EXPECT_THAT(housl(m, n, xvector, incx, tau, a, lda, &ier), Eq(0));
270 1 : EXPECT_THAT(ier, Eq(0));
271 1 : }
272 :
273 2 : TEST_F(TestThatHouseholderL, ImplementsHouseholderMultiplicationOnLeft)
274 : {
275 1 : xvector[ 0 ] = a(1, 1);
276 1 : xvector[ 1 ] = a(2, 1);
277 1 : xvector[ 2 ] = a(3, 1);
278 1 : xvector[ 3 ] = a(4, 1);
279 :
280 1 : double *alpha = xvector;
281 1 : double *v = xvector + 1;
282 :
283 : // On input of housg, alpha is the first component of xvector
284 1 : EXPECT_THAT(*alpha, DoubleEq(xvector[ 0 ]));
285 : // On input of housg, v is components 2:m of xvector
286 4 : for (int row = 0; row < m - 1; row++)
287 3 : EXPECT_THAT(*(v + row), DoubleEq(xvector[ row + 1 ]));
288 :
289 1 : EXPECT_THAT(housg(m, alpha, v, incx, &tau, safmin), Eq(0));
290 :
291 : // On output of housg, beta is stored in alpha
292 1 : double beta = *alpha;
293 : // On output of housg, beta is stored in the first component of xvector
294 1 : EXPECT_THAT(beta, DoubleEq(xvector[ 0 ]));
295 : // On output of housg, xvector is components 2:m of v
296 4 : for (int row = 0; row < m - 1; row++)
297 3 : EXPECT_THAT(*(v + row), DoubleEq(xvector[ row + 1 ]));
298 :
299 1 : int i = 1;
300 1 : EXPECT_THAT(
301 : housl(m - i + 1, n - i, xvector, incx, tau, a + map(i, i + 1), lda, &ier),
302 : Eq(0));
303 :
304 : // First column remains a(:,1) here which does not hold (beta v)
305 1 : r(1, 1) = 1.0, r(1, 2) = -13.0, r(1, 3) = -15.0;
306 1 : r(2, 1) = 1.0, r(2, 2) = 0.0, r(2, 3) = 0.0;
307 1 : r(3, 1) = 1.0, r(3, 2) = 3.0, r(3, 3) = 3.0;
308 1 : r(4, 1) = 1.0, r(4, 2) = 6.0, r(4, 3) = 6.0;
309 :
310 4 : for (int col = 1; col <= n; col++)
311 15 : for (int row = 1; row <= m; row++)
312 12 : EXPECT_THAT(a(row, col), DOUBLE_NEAR_MULTIPLIER(r(row, col), 100));
313 1 : }
314 :
315 2 : TEST_F(TestThatHouseholderL,
316 : ImplementsHouseholderMultiplicationOnLeftWithoutAdditionalStorage)
317 : {
318 1 : double *alpha = a + map(1, 1);
319 1 : double *v = a + map(2, 1);
320 :
321 : // On input of housg, alpha is the first component of xvector
322 1 : EXPECT_THAT(*alpha, DoubleEq(a(1, 1)));
323 : // On input of housg, v is components 2:m of xvector
324 4 : for (int row = 0; row < m - 1; row++)
325 3 : EXPECT_THAT(*(v + row), DoubleEq(a(row + 2, 1)));
326 :
327 1 : EXPECT_THAT(housg(m, alpha, v, incx, &tau, safmin), Eq(0));
328 :
329 : // On output of housg, beta is stored in alpha
330 1 : double beta = *alpha;
331 : // On output of housg, beta is stored in the first component of xvector
332 1 : EXPECT_THAT(beta, DoubleEq(a(1, 1)));
333 : // On output of housg, xvector is components 2:m of v
334 4 : for (int row = 0; row < m - 1; row++)
335 3 : EXPECT_THAT(*(v + row), DoubleEq(a(row + 2, 1)));
336 :
337 1 : int i = 1;
338 1 : EXPECT_THAT(housl(m - i + 1, n - i, a + map(i, i), incx, tau,
339 : a + map(i, i + 1), lda, &ier),
340 : Eq(0));
341 :
342 : // First column remains a(:,1) holds (beta v)
343 1 : r(1, 1) = beta, r(1, 2) = -13.0, r(1, 3) = -15.0;
344 1 : r(2, 1) = *(v + 0), r(2, 2) = 0.0, r(2, 3) = 0.0;
345 1 : r(3, 1) = *(v + 1), r(3, 2) = 3.0, r(3, 3) = 3.0;
346 1 : r(4, 1) = *(v + 2), r(4, 2) = 6.0, r(4, 3) = 6.0;
347 :
348 4 : for (int col = 1; col <= n; col++)
349 15 : for (int row = 1; row <= m; row++)
350 12 : EXPECT_THAT(a(row, col), DOUBLE_NEAR_MULTIPLIER(r(row, col), 100));
351 1 : }
352 :
353 : class TestThatQRF : public Test
354 : {
355 : public:
356 : int m = 4, lda = 4, n = 3, min_m_n = 0, ier = 0;
357 : double *a = NULL;
358 : int *ipiv = NULL;
359 : double *tau = NULL;
360 : double *wrk = NULL;
361 : double *r = NULL;
362 : double epmach = DBL_EPSILON;
363 : double safmin = DBL_MIN;
364 :
365 7 : void SetUp() override
366 : {
367 7 : a = (double *)calloc((size_t)(lda * n), (size_t)sizeof(double));
368 : // if (a == NULL)
369 : //{
370 : // printf("Failed to allocate matrix A memory!\n");
371 : // assert(a);
372 : //}
373 :
374 7 : a(1, 1) = 1.0, a(1, 2) = 2.0, a(1, 3) = 3.0;
375 7 : a(2, 1) = 1.0, a(2, 2) = 5.0, a(2, 3) = 6.0;
376 7 : a(3, 1) = 1.0, a(3, 2) = 8.0, a(3, 3) = 9.0;
377 7 : a(4, 1) = 1.0, a(4, 2) = 11.0, a(4, 3) = 12.0;
378 :
379 7 : min_m_n = (m < n) ? m : n;
380 :
381 7 : tau = (double *)calloc((size_t)min_m_n, (size_t)sizeof(double));
382 : // if (tau == NULL)
383 : //{
384 : // printf("Failed to allocate tau memory!\n");
385 : // assert(tau);
386 : //}
387 :
388 7 : wrk = (double *)calloc((size_t)n, (size_t)sizeof(double));
389 : // if (wrk == NULL)
390 : //{
391 : // printf("Failed to allocate wrk memory!\n");
392 : // assert(wrk);
393 : //}
394 :
395 7 : ipiv = (int *)calloc((size_t)n, (size_t)sizeof(int));
396 : // if (ipiv == NULL)
397 : //{
398 : // printf("Failed to allocate ipiv memory!\n");
399 : // assert(ipiv);
400 : //}
401 :
402 7 : r = (double *)calloc((size_t)(lda * n), (size_t)sizeof(double));
403 : // if (r == NULL)
404 : //{
405 : // printf("Failed to allocate space for result r!\n");
406 : // assert(r);
407 : //}
408 7 : }
409 :
410 7 : void TearDown() override
411 : {
412 7 : free(r);
413 7 : free(ipiv);
414 7 : free(wrk);
415 7 : free(tau);
416 7 : free(a);
417 7 : }
418 : };
419 :
420 2 : TEST_F(TestThatQRF, CanAccessTheQRFactorizationMethod)
421 : {
422 1 : ASSERT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(0));
423 : }
424 :
425 2 : TEST_F(TestThatQRF, IsNotDefinedForNegativeRowDimensionM)
426 : {
427 1 : m = -1;
428 1 : EXPECT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(1));
429 1 : EXPECT_THAT(ier, Eq(1));
430 1 : }
431 :
432 2 : TEST_F(TestThatQRF, IsNotDefinedForNegativeColumnDimensionN)
433 : {
434 1 : n = -1;
435 1 : EXPECT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(1));
436 1 : EXPECT_THAT(ier, Eq(1));
437 1 : }
438 :
439 2 : TEST_F(TestThatQRF, IsNotDefinedForLeadingDimensionLDALessThanNumberOfRowsM)
440 : {
441 1 : lda = 1;
442 1 : EXPECT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(1));
443 1 : EXPECT_THAT(ier, Eq(1));
444 1 : }
445 :
446 2 : TEST_F(TestThatQRF, IsDefinedForRowDimensionMEqualsZero)
447 : {
448 1 : m = 0;
449 1 : EXPECT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(0));
450 1 : EXPECT_THAT(ier, Eq(0));
451 1 : }
452 :
453 2 : TEST_F(TestThatQRF, IsDefinedForColumnDimensionNEqualsZero)
454 : {
455 1 : n = 0;
456 1 : EXPECT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(0));
457 1 : EXPECT_THAT(ier, Eq(0));
458 1 : }
459 :
460 2 : TEST_F(TestThatQRF, WillQRFactorA)
461 : {
462 1 : EXPECT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(0));
463 1 : EXPECT_THAT(ier, Eq(0));
464 :
465 : // >> M = [[1 2 3];[1 5 6];[1 8 9];[1 11 12]]
466 : // M =
467 : // 1 2 3
468 : // 1 5 6
469 : // 1 8 9
470 : // 1 11 12
471 : // >> [Q,R,E] = qr(M)
472 : // Q =
473 : // Columns 1 through 3
474 : // -1.825741858350554e-01 -8.164965809277261e-01 5.477210117775342e-01
475 : // -3.651483716701108e-01 -4.082482904638630e-01 -7.312645785255080e-01
476 : // -5.477225575051662e-01 -8.814441016525298e-17 -1.806338782815866e-01
477 : // -7.302967433402215e-01 4.082482904638629e-01 3.641774450295606e-01
478 : // Column 4
479 : // 1.301252240837732e-03
480 : // 4.065121353587260e-01
481 : // -8.169280274399654e-01
482 : // 4.091146398404016e-01
483 : // R =
484 : // -1.643167672515498e+01 -1.825741858350554e+00 -1.460593486680443e+01
485 : // 0 -8.164965809277261e-01 8.164965809277228e-01
486 : // 0 0 2.174446367394293e-15
487 : // 0 0 0
488 : // E =
489 : // 0 1 0
490 : // 0 0 1
491 : // 1 0 0
492 : // >> transpose(Q)*Q
493 : // ans =
494 : // Columns 1 through 3
495 : // 1.000000000000000e+00 1.110223024625157e-16 -7.307522642552300e-17
496 : // 1.110223024625157e-16 1.000000000000000e+00 -1.110223024625157e-16
497 : // -7.307522642552300e-17 -1.110223024625157e-16 1.000000000000000e+00
498 : // 0 8.283304597789254e-17 8.326672684688674e-17
499 : // Column 4
500 : // 0
501 : // 8.283304597789254e-17
502 : // 8.326672684688674e-17
503 : // 9.999999999999999e-01
504 : // >> Q*R-M*E
505 : // ans =
506 : // -4.440892098500626e-16 0 4.440892098500626e-15
507 : // 0 0 8.881784197001252e-16
508 : // 0 0 1.776356839400250e-15
509 : // 0 0 1.776356839400250e-15
510 : // Check R
511 1 : r(1, 1) = -1.643167672515498e+01;
512 1 : r(1, 2) = -1.825741858350554e+00;
513 1 : r(1, 3) = -1.460593486680443e+01;
514 1 : r(2, 1) = 0;
515 1 : r(2, 2) = -8.164965809277261e-01;
516 1 : r(2, 3) = 8.164965809277228e-01;
517 1 : r(3, 1) = 0;
518 1 : r(3, 2) = 0;
519 1 : r(3, 3) = 2.174446367394293e-15;
520 1 : r(4, 1) = 0;
521 1 : r(4, 2) = 0;
522 1 : r(4, 3) = 0;
523 4 : for (int col = 1; col <= n; col++)
524 9 : for (int row = 1; row <= col; row++)
525 6 : EXPECT_THAT(a(row, col), DOUBLE_NEAR_MULTIPLIER(r(row, col), 100));
526 : // Check E
527 1 : EXPECT_THAT(ipiv[ 0 ], Eq(3));
528 1 : EXPECT_THAT(ipiv[ 1 ], Eq(1));
529 1 : EXPECT_THAT(ipiv[ 2 ], Eq(2));
530 1 : }
531 :
532 : class TestThatQRS : public Test
533 : {
534 : public:
535 : int m = 4, lda = 4, n = 2, min_m_n = 0, ier = 0;
536 : double *a = NULL;
537 : int *ipiv = NULL;
538 : double *tau = NULL;
539 : double *wrk = NULL;
540 : double *r = NULL;
541 : double *x = NULL;
542 : double *y = NULL;
543 : double epmach = DBL_EPSILON;
544 : double safmin = DBL_MIN;
545 :
546 8 : void SetUp() override
547 : {
548 8 : a = (double *)calloc((size_t)(lda * n), (size_t)sizeof(double));
549 : // if (a == NULL)
550 : //{
551 : // printf("Failed to allocate matrix A memory!\n");
552 : // assert(a);
553 : //}
554 :
555 8 : a(1, 1) = 1.0, a(1, 2) = 2.0;
556 8 : a(2, 1) = 1.0, a(2, 2) = 5.0;
557 8 : a(3, 1) = 1.0, a(3, 2) = 8.0;
558 8 : a(4, 1) = 1.0, a(4, 2) = 11.0;
559 :
560 8 : min_m_n = (m < n) ? m : n;
561 :
562 8 : tau = (double *)calloc((size_t)min_m_n, (size_t)sizeof(double));
563 : // if (tau == NULL)
564 : //{
565 : // printf("Failed to allocate tau memory!\n");
566 : // assert(tau);
567 : //}
568 :
569 8 : wrk = (double *)calloc((size_t)n, (size_t)sizeof(double));
570 : // if (wrk == NULL)
571 : //{
572 : // printf("Failed to allocate wrk memory!\n");
573 : // assert(wrk);
574 : //}
575 :
576 8 : ipiv = (int *)calloc((size_t)n, (size_t)sizeof(int));
577 : // if (ipiv == NULL)
578 : //{
579 : // printf("Failed to allocate ipiv memory!\n");
580 : // assert(ipiv);
581 : //}
582 :
583 8 : r = (double *)calloc((size_t)(lda * n), (size_t)sizeof(double));
584 : // if (r == NULL)
585 : //{
586 : // printf("Failed to allocate space for result r!\n");
587 : // assert(r);
588 : //}
589 :
590 8 : x = (double *)calloc((size_t)n, (size_t)sizeof(double));
591 : // if (x == NULL)
592 : //{
593 : // printf("Failed to allocate space for result x!\n");
594 : // assert(x);
595 : //}
596 :
597 8 : y = (double *)calloc((size_t)m, (size_t)sizeof(double));
598 : // if (y == NULL)
599 : //{
600 : // printf("Failed to allocate space for result y!\n");
601 : // assert(y);
602 : //}
603 8 : }
604 :
605 8 : void TearDown() override
606 : {
607 8 : free(y);
608 8 : free(x);
609 8 : free(r);
610 8 : free(ipiv);
611 8 : free(wrk);
612 8 : free(tau);
613 8 : free(a);
614 8 : }
615 : };
616 :
617 2 : TEST_F(TestThatQRS, CanAccessTheLeastSquaredSolutionMethod)
618 : {
619 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(0));
620 1 : EXPECT_THAT(ier, Eq(0));
621 1 : }
622 :
623 2 : TEST_F(TestThatQRS, IsNotDefinedForNegativeRowDimensionM)
624 : {
625 1 : m = -1;
626 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(1));
627 1 : EXPECT_THAT(ier, Eq(1));
628 1 : }
629 :
630 2 : TEST_F(TestThatQRS, IsNotDefinedForNegativeColumnDimensionN)
631 : {
632 1 : n = -1;
633 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(1));
634 1 : EXPECT_THAT(ier, Eq(1));
635 1 : }
636 :
637 2 : TEST_F(TestThatQRS, IsNotDefinedForRowDimensionMLessThanNumberColumnDimensionN)
638 : {
639 1 : m = 1;
640 1 : n = 2;
641 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(1));
642 1 : EXPECT_THAT(ier, Eq(1));
643 1 : }
644 :
645 2 : TEST_F(TestThatQRS, IsNotDefinedForLeadingDimensionLDALessThanNumberOfRowsM)
646 : {
647 1 : lda = 1;
648 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(1));
649 1 : EXPECT_THAT(ier, Eq(1));
650 1 : }
651 :
652 2 : TEST_F(TestThatQRS, IsDefinedForRowDimensionMEqualsZero)
653 : {
654 1 : m = 0;
655 1 : n = 0;
656 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(0));
657 1 : EXPECT_THAT(ier, Eq(0));
658 1 : }
659 :
660 2 : TEST_F(TestThatQRS, IsDefinedForColumnDimensionNEqualsZero)
661 : {
662 1 : n = 0;
663 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(0));
664 1 : EXPECT_THAT(ier, Eq(0));
665 1 : }
666 :
667 2 : TEST_F(TestThatQRS, WillQRFactorANDComputeLinearLeastSquarsSolution)
668 : {
669 1 : EXPECT_THAT(qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier), Eq(0));
670 1 : EXPECT_THAT(ier, Eq(0));
671 1 : y(1) = 3.0;
672 1 : y(2) = 5.0;
673 1 : y(3) = 2.0;
674 1 : y(4) = 1.0;
675 1 : EXPECT_THAT(qrs(m, n, a, lda, tau, y, x, &ier), Eq(0));
676 1 : EXPECT_THAT(ier, Eq(0));
677 1 : EXPECT_THAT(x(ipiv(1)), DOUBLE_NEAR_MULTIPLIER(4.7, 100.0));
678 1 : EXPECT_THAT(x(ipiv(2)), DOUBLE_NEAR_MULTIPLIER(-0.3, 100.0));
679 1 : }
680 :
681 : class TestThatQR : public Test
682 : {
683 : public:
684 : int m{3}, lda{3}, n{3}, min_m_n{0}, ier{0}, ldq{3};
685 : vector<int> ipiv{0, 0, 0};
686 : vector<double> a{2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0};
687 : vector<double> q{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
688 : vector<double> tau{0.0, 0.0, 0.0};
689 : vector<double> wrk{0.0, 0.0, 0.0};
690 : vector<double> x{0.0, 0.0, 0.0};
691 : vector<double> y{-1.0, -1.0, -1.0};
692 : double epmach{DBL_EPSILON};
693 : double safmin{DBL_MIN};
694 :
695 : vector<double> a_check{
696 : 2.4494897427831783, -0.57979589711327117, 0.28989794855663559,
697 : -1.6329931618554516, -1.5275252316519461, -0.39985928207919952,
698 : -1.6329931618554518, 1.0910894511799623, 1.0690449676496974};
699 : vector<double> q_check{
700 : -0.40824829046386291, 0.81649658092772592, -0.40824829046386296,
701 : -0.87287156094396945, -0.21821789023599258, 0.43643578047198495,
702 : 0.26726124191242445, 0.53452248382484879, 0.80178372573727308};
703 : vector<int> ipiv_check{2, 1, 3};
704 : vector<double> x_check{-1.5, -2, -1.5};
705 :
706 9 : void SetUp() override { min_m_n = (m < n) ? m : n; }
707 :
708 9 : void TearDown() override {}
709 : };
710 :
711 4 : inline void EXPECT_DARRAY_EQ(const int n, double *expected, double *actual)
712 : {
713 40 : for (int i = 0; i < n; i++)
714 : {
715 36 : EXPECT_THAT(*(actual + i), DOUBLE_NEAR_MULTIPLIER(*(expected + i), 100.0));
716 : }
717 4 : }
718 :
719 1 : inline void EXPECT_IARRAY_EQ(const int n, int *expected, int *actual)
720 : {
721 4 : for (int i = 0; i < n; i++)
722 : {
723 3 : EXPECT_THAT(*(actual + i), Eq(*(expected + i)));
724 : }
725 1 : }
726 :
727 2 : TEST_F(TestThatQR, WillFactorANDComputeSolutionForSquareSystem)
728 : {
729 1 : EXPECT_THAT(
730 : qrf(m, n, &a[ 0 ], lda, &ipiv[ 0 ], &tau[ 0 ], &wrk[ 0 ], safmin, &ier),
731 : Eq(0));
732 1 : EXPECT_THAT(ier, Eq(0));
733 :
734 1 : EXPECT_DARRAY_EQ(m * n, &a[ 0 ], &a_check[ 0 ]);
735 1 : EXPECT_IARRAY_EQ(m, &ipiv[ 0 ], &ipiv_check[ 0 ]);
736 :
737 1 : EXPECT_THAT(qorg(m, n, 1, m, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
738 : Eq(0));
739 1 : EXPECT_THAT(ier, Eq(0));
740 :
741 1 : EXPECT_DARRAY_EQ(m * n, &a[ 0 ], &a_check[ 0 ]);
742 1 : EXPECT_DARRAY_EQ(m * n, &q[ 0 ], &q_check[ 0 ]);
743 :
744 1 : EXPECT_THAT(qrs(m, n, &a[ 0 ], lda, &tau[ 0 ], &y[ 0 ], &x[ 0 ], &ier),
745 : Eq(0));
746 1 : EXPECT_THAT(ier, Eq(0));
747 :
748 1 : EXPECT_DARRAY_EQ(m * n, &a[ 0 ], &a_check[ 0 ]);
749 4 : for (int i = 0; i < n; i++)
750 : {
751 3 : EXPECT_THAT(x[ ipiv[ i ] - 1 ],
752 : DOUBLE_NEAR_MULTIPLIER(x_check[ i ], 100.0));
753 : }
754 1 : }
755 :
756 2 : TEST_F(TestThatQR, WillProvideOrthogonalMartixQUnlessMLessThanZero)
757 : {
758 1 : int ier{0};
759 1 : int m{-1};
760 1 : int n{2};
761 1 : int k{1};
762 1 : int ell{m};
763 1 : int lda{1};
764 1 : int ldq{2};
765 : // All matricies not accessed in this test
766 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
767 : Eq(1));
768 1 : }
769 :
770 2 : TEST_F(TestThatQR, WillProvideOrthogonalMartixQUnlessNLessThanZero)
771 : {
772 1 : int ier{0};
773 1 : int m{2};
774 1 : int n{-1};
775 1 : int k{1};
776 1 : int ell{m};
777 1 : int lda{1};
778 1 : int ldq{2};
779 : // All matricies not accessed in this test
780 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
781 : Eq(1));
782 1 : }
783 :
784 2 : TEST_F(TestThatQR, WillProvideOrthogonalMartixQUnlessMLessThanN)
785 : {
786 1 : int ier{0};
787 1 : int m{2};
788 1 : int n{3};
789 1 : int k{1};
790 1 : int ell{m};
791 1 : int lda{m};
792 1 : int ldq{m};
793 : // All matricies not accessed in this test
794 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
795 : Eq(1));
796 1 : }
797 :
798 2 : TEST_F(TestThatQR, WillProvideOrthogonalMartixQUnlessKGreaterThanL)
799 : {
800 1 : int ier{0};
801 1 : int m{3};
802 1 : int n{2};
803 1 : int k{m + 1};
804 1 : int ell{m};
805 1 : int lda{m};
806 1 : int ldq{m};
807 : // All matricies not accessed in this test
808 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
809 : Eq(1));
810 1 : }
811 :
812 2 : TEST_F(TestThatQR, WillProvideOrthogonalMartixQUnlessLDALessThanM)
813 : {
814 1 : int ier{0};
815 1 : int m{3};
816 1 : int n{2};
817 1 : int k{m};
818 1 : int ell{m};
819 1 : int lda{m - 1};
820 1 : int ldq{m};
821 : // All matricies not accessed in this test
822 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
823 : Eq(1));
824 1 : }
825 :
826 2 : TEST_F(TestThatQR, WillProvideOrthogonalMartixQUnlessLDQLessThanOne)
827 : {
828 1 : int ier{0};
829 1 : int m{3};
830 1 : int n{2};
831 1 : int k{m};
832 1 : int ell{m};
833 1 : int lda{m};
834 1 : int ldq{0};
835 : // All matricies not accessed in this test
836 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
837 : Eq(1));
838 1 : }
839 :
840 2 : TEST_F(TestThatQR, WillProvideOrthogonalMartixQWithQuickReturnNEqualsZero)
841 : {
842 1 : int ier{0};
843 1 : int m{3};
844 1 : int n{0};
845 1 : int k{m};
846 1 : int ell{m};
847 1 : int lda{m};
848 1 : int ldq{m};
849 : // All matricies not accessed in this test
850 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
851 : Eq(0));
852 1 : }
853 :
854 2 : TEST_F(TestThatQR,
855 : WillProvideOrthogonalMartixQWithQuickReturnNEqualsZeroANDMEqualsZero)
856 : {
857 1 : int ier{0};
858 1 : int m{0};
859 1 : int n{0};
860 1 : int k{m};
861 1 : int ell{m};
862 1 : int lda{1};
863 1 : int ldq{1};
864 : // All matricies not accessed in this test
865 1 : EXPECT_THAT(qorg(m, n, k, ell, &a[ 0 ], lda, &tau[ 0 ], &q[ 0 ], ldq, &ier),
866 : Eq(0));
867 1 : }
|
