1 // This is core/vnl/algo/tests/test_sparse_lu.cxx
2 #include <iostream>
3 #include <testlib/testlib_test.h>
4 #include <vnl/vnl_sparse_matrix.h>
5 #include <vnl/algo/vnl_sparse_lu.h>
6
test_sparse_lu()7 static void test_sparse_lu()
8 {
9 //mat0 of Kenneth S. Kunder's Sparse 1.3a release
10 vnl_sparse_matrix<double> A(4,4);
11 std::vector<int> cols0(2), cols1(3), cols2(3), cols3(2);
12 std::vector<double> vals0(2), vals1(3), vals2(3), vals3(2);
13 cols0[0]=0; cols0[1]=1;
14 vals0[0]=2.0; vals0[1]=-1.0;
15 A.set_row(0, cols0, vals0);
16 cols1[0]=0; cols1[1]=1; cols1[2]=2;
17 vals1[0]=-1.0; vals1[1]=3.0; vals1[2]=-1;
18 A.set_row(1, cols1, vals1);
19 cols2[0]=1; cols2[1]= 2; cols2[2]= 3;
20 vals2[0]=-1.0; vals2[1]=3.0; vals2[2]=-1.0;
21 A.set_row(2, cols2, vals2);
22 cols3[0]=2; cols3[1]=3;
23 vals3[0]=-1.0; vals3[1]=3.0;
24 A.set_row(3, cols3, vals3);
25 for (A.reset(); A.next();)
26 std::cout << "A[" << A.getrow() << "][" << A.getcolumn()
27 << "]= " << A.value() << '\n';
28 vnl_vector<double> b(4, 0.0), x(4);
29 b[0]=34.0;
30 vnl_sparse_lu lu(A,vnl_sparse_lu::verbose);
31 lu.solve(b, &x);
32 for (unsigned i = 0; i<4; ++i)
33 std::cout << "x[" << i << "]= " << x[i] << '\n';
34 TEST_NEAR("solution of mat0 example", x[0], 21, 1.e-03);
35 double det = lu.determinant();
36 std::cout << "determinant = " << det << '\n';
37 TEST_NEAR("determinant of mat0 example", det, 34, 1.e-03);
38 lu.solve_transpose(b,&x);
39 std::cout << "transpose solution\n";
40 for (unsigned i = 0; i<4; ++i)
41 std::cout << "x[" << i << "]= " << x[i] << '\n';
42 TEST_NEAR("transpose solution of mat0 example", x[2], 3, 1.e-03);
43 //mat5 of sparse test data
44 vnl_sparse_matrix<double> Ap(3,3);
45 Ap(0,1)=1; Ap(1,2)=1; Ap(2,0)=1;
46 vnl_vector<double> bp(3), xp(3);
47 bp[0]=2.0; bp[1]=3.0; bp[2]=1.0;
48 vnl_sparse_lu lup(Ap,vnl_sparse_lu::verbose);
49 lup.solve(bp, &xp);
50 for (unsigned i = 0; i<3; ++i)
51 std::cout << "xp[" << i << "]= " << xp[i] << '\n';
52 TEST_NEAR("solution of mat5 example", xp[2], 3, 1.e-03);
53
54 //test matrix derived from Poisson birth-death queue
55 double s = -0.01, l = 0.5, m = 0.5;
56 vnl_sparse_matrix<double> S(6,6);
57 S(0,0)=s+l; S(0,1)=-l;
58 S(1,0)=-m; S(1,1)=s+l+m; S(1,2)=-l;
59 S(2,1)=-m; S(2,2)=s+l+m;
60 S(3,3)=s+l+m; S(3,4)=-l;
61 S(4,3)=-m; S(4,4)=s+l+m; S(4,5)=-l;
62 S(5,4)=-m; S(5,5)=m+s;
63 vnl_vector<double> bbd(6),xbd(6);
64 bbd[0]=0; bbd[1]=0; bbd[2]=l; bbd[3]=m; bbd[4]=0; bbd[5]=0;
65 vnl_sparse_lu lubd(S,vnl_sparse_lu::estimate_condition_verbose);
66 lubd.solve(bbd, &xbd);
67 for (unsigned i = 0; i<6; ++i)
68 std::cout << "xbd[" << i << "]= " << xbd[i] << '\n';
69 TEST_NEAR("test solution of birth-death matrix", xbd[2], 1.06622, 1.e-04);
70 det = lubd.determinant();
71 std::cout << "birth-death determinant = " << det << '\n';
72 double cond = lubd.rcond();
73 std::cout << "birth-death condition number = " << cond << '\n';
74 TEST_NEAR("birth-death matrix condition number", cond, 0.03756, 1.e-04);
75 double upbnd = lubd.max_error_bound();
76 std::cout << "birth-death upper error bound = " << upbnd << '\n';
77 TEST_NEAR("birth-death upper error", upbnd, 5.923e-015, 1.e-016);
78 }
79
80 TESTMAIN(test_sparse_lu);
81