1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 Bradley M. Bell
3 
4 CppAD is distributed under the terms of the
5              Eclipse Public License Version 2.0.
6 
7 This Source Code may also be made available under the following
8 Secondary License when the conditions for such availability set forth
9 in the Eclipse Public License, Version 2.0 are satisfied:
10       GNU General Public License, Version 2.0 or later.
11 ---------------------------------------------------------------------------- */
12 
13 /*
14 $begin sparse2eigen.cpp$$
15 $spell
16     Eigen
17     Cpp
18 $$
19 
20 $section Converting CppAD Sparse Matrix to Eigen Format: Example and Test$$
21 
22 
23 $srcthisfile%0%// BEGIN C++%// END C++%1%$$
24 
25 $end
26 */
27 // BEGIN C++
28 # include <cppad/utility/sparse2eigen.hpp>
sparse2eigen(void)29 bool sparse2eigen(void)
30 {   bool ok = true;
31     //
32     typedef CppAD::vector<size_t>                 s_vector;
33     typedef CppAD::vector<double>                 d_vector;
34     typedef CppAD::sparse_rc<s_vector>            sparse_rc;
35     typedef CppAD::sparse_rcv<s_vector, d_vector> sparse_rcv;
36     //
37     // sparsity pattern
38     size_t nr  = 3;
39     size_t nc  = 4;
40     size_t nnz = 2 * nr;
41     sparse_rc pattern(nr, nc, nnz);
42     for(size_t i = 0; i < nr; i++)
43     {   for(size_t j = i; j <= i + 1; j++)
44         {   size_t k = i + j;
45             pattern.set(k, i, j);
46         }
47     }
48     //
49     // sparse matrix
50     sparse_rcv source(pattern);
51     for(size_t i = 0; i < nr; i++)
52     {   for(size_t j = i; j <= i + 1; j++)
53         {   size_t k = i + j;
54             double v = double(k) / 2.0;
55             source.set(k, v);
56         }
57     }
58     //
59     //  example using row major order
60     {   Eigen::SparseMatrix<double, Eigen::RowMajor> destination;
61         CppAD::sparse2eigen(source, destination);
62         //
63         typedef
64         Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator iterator;
65         //
66         // check result
67         const double* d_value = destination.valuePtr();
68         size_t k = 0;
69         for(int i = 0; i < destination.outerSize(); ++i)
70         {
71             for(iterator itr(destination, i); itr; ++itr)
72             {   // check row
73                 ok &= itr.row() == i;
74                 //
75                 // check column
76                 if( k % 2 == 0 )
77                     ok &= itr.col() == i;
78                 else
79                     ok &= itr.col() == (i+1);
80                 //
81                 // check value
82                 ok  &= itr.value() == double( itr.row() + itr.col() ) / 2.0;
83                 ok  &= itr.value() == d_value[k];
84                 //
85                 ++k;
86             }
87         }
88         ok &= k == nnz;
89     }
90     //
91     //  example using column major order
92     {   Eigen::SparseMatrix<double, Eigen::ColMajor> destination;
93         CppAD::sparse2eigen(source, destination);
94         //
95         typedef
96         Eigen::SparseMatrix<double, Eigen::ColMajor>::InnerIterator iterator;
97         //
98         // check result
99         const double* d_value = destination.valuePtr();
100         size_t k = 0;
101         for(int j = 0; j < destination.outerSize(); ++j)
102         {   for(iterator itr(destination, j); itr; ++itr)
103             {   // check column
104                 ok &= itr.col() == j;
105                 //
106                 // check row
107                 if( j == 0 )
108                 {   assert( k == 0 );
109                     ok &= itr.row() == 0;
110                 }
111                 else if( k % 2 == 1 )
112                     ok &= itr.row() == j - 1;
113                 else
114                     ok &= itr.row() == j;
115                 //
116                 // check value
117                 ok  &= itr.value() == double( itr.row() + itr.col() ) / 2.0;
118                 ok  &= itr.value() == d_value[k];
119                 //
120                 ++k;
121             }
122         }
123         ok &= k == nnz;
124     }
125     //
126     return ok;
127 }
128 // END C++
129