1 #include <stan/math/prim.hpp>
2 #include <gtest/gtest.h>
3
generate_large_L_tri_mat()4 stan::math::matrix_d generate_large_L_tri_mat() {
5 stan::math::matrix_d x;
6 double vals[10000];
7
8 vals[0] = 0.1;
9 for (int i = 1; i < 10000; ++i)
10 vals[i] = vals[i - 1] + 0.1123456;
11
12 x = Eigen::Map<Eigen::Matrix<double, 100, 100> >(vals);
13 x *= 1e10;
14
15 return x;
16 }
17
test_multiply_lower_tri_self_transpose(const stan::math::matrix_d & x)18 void test_multiply_lower_tri_self_transpose(const stan::math::matrix_d& x) {
19 using stan::math::multiply_lower_tri_self_transpose;
20 stan::math::matrix_d y = multiply_lower_tri_self_transpose(x);
21 stan::math::matrix_d xp = x;
22 for (int m = 0; m < xp.rows(); ++m)
23 for (int n = m + 1; n < xp.cols(); ++n)
24 xp(m, n) = 0;
25
26 stan::math::matrix_d xxt = xp * xp.transpose();
27 EXPECT_EQ(y.rows(), xxt.rows());
28 EXPECT_EQ(y.cols(), xxt.cols());
29 for (int m = 0; m < y.rows(); ++m)
30 for (int n = 0; n < y.cols(); ++n)
31 EXPECT_FLOAT_EQ(xxt(m, n), y(m, n));
32 }
33
TEST(MathMatrixPrimMat,multiply_lower_tri_self_transpose)34 TEST(MathMatrixPrimMat, multiply_lower_tri_self_transpose) {
35 using stan::math::check_symmetric;
36 using stan::math::multiply_lower_tri_self_transpose;
37 static const char* function
38 = "stan::math::multiply_lower_tri_self_transpose(%1%)";
39 stan::math::matrix_d x;
40 test_multiply_lower_tri_self_transpose(x);
41
42 x = stan::math::matrix_d(1, 1);
43 x << 3.0;
44 test_multiply_lower_tri_self_transpose(x);
45
46 x = stan::math::matrix_d(2, 2);
47 x << 1.0, 0.0, 2.0, 3.0;
48 test_multiply_lower_tri_self_transpose(x);
49
50 x = stan::math::matrix_d(3, 3);
51 x << 1.0, 0.0, 0.0, 2.0, 3.0, 0.0, 4.0, 5.0, 6.0;
52 test_multiply_lower_tri_self_transpose(x);
53
54 x = stan::math::matrix_d(3, 3);
55 x << 1.0, 0.0, 100000.0, 2.0, 3.0, 0.0, 4.0, 5.0, 6.0;
56 test_multiply_lower_tri_self_transpose(x);
57
58 x = stan::math::matrix_d(3, 2);
59 x << 1.0, 0.0, 2.0, 3.0, 4.0, 5.0;
60 test_multiply_lower_tri_self_transpose(x);
61
62 x = stan::math::matrix_d(2, 3);
63 x << 1.0, 0.0, 0.0, 2.0, 3.0, 0.0;
64 test_multiply_lower_tri_self_transpose(x);
65
66 x = generate_large_L_tri_mat();
67 EXPECT_NO_THROW(check_symmetric(function, "Symmetric matrix",
68 multiply_lower_tri_self_transpose(x)));
69 }
70