source: trunk/speed/adolc/sparse_jacobian.cpp @ 3138

Last change on this file since 3138 was 3138, checked in by bradbell, 6 years ago

Change link_sparse_jacobian to return n_sweep and report it for all sizes.

  • Property svn:keywords set to Id
File size: 6.0 KB
Line 
1/* $Id: sparse_jacobian.cpp 3138 2014-03-02 18:46:11Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12# include <cstring>
13# include <cppad/vector.hpp>
14
15/*
16$begin adolc_sparse_jacobian.cpp$$
17$spell
18        const
19        sparsedrivers.cpp
20        colpack
21        boolsparsity
22        adouble
23        int int_n
24        cppad.hpp
25        onetape
26        typedef
27        alloc
28        jac
29        nnz
30        cind
31        bool
32        CppAD
33        adolc
34        sparse_jacobian
35$$
36
37$section adolc Speed: Sparse Jacobian$$
38
39$index link_sparse_jacobian, adolc$$
40$index adolc, link_sparse_jacobian$$
41$index speed, adolc$$
42$index adolc, speed$$
43$index sparse, speed adolc$$
44$index jacobian, speed adolc$$
45
46$head Specifications$$
47See $cref link_sparse_jacobian$$.
48
49$head Implementation$$
50
51$codep */
52# include <adolc/adolc.h>
53# include <adolc/adolc_sparse.h>
54# include <cppad/vector.hpp>
55# include <cppad/speed/uniform_01.hpp>
56# include <cppad/speed/sparse_jac_fun.hpp>
57
58// list of possible options
59extern bool global_memory, global_onetape, global_atomic, global_optimize;
60extern bool global_colpack, global_boolsparsity;
61
62bool link_sparse_jacobian(
63        size_t                           size     , 
64        size_t                           repeat   , 
65        size_t                           m        ,
66        const CppAD::vector<size_t>&     row      ,
67        const CppAD::vector<size_t>&     col      ,
68              CppAD::vector<double>&     x_return ,
69              CppAD::vector<double>&     jacobian ,
70              size_t&                    n_sweep  )
71{
72        if( global_atomic || (! global_colpack) )
73                return false; 
74        if( global_memory || global_optimize )
75                return false; 
76        // -----------------------------------------------------
77        // setup
78        typedef unsigned int*    SizeVector;
79        typedef double*          DblVector;
80        typedef adouble          ADScalar;
81        typedef ADScalar*        ADVector;
82
83        size_t i, j;                // temporary indices
84        size_t n = size;            // number of independent variables
85        size_t order = 0;          // derivative order corresponding to function
86
87        // set up for thread_alloc memory allocator (fast and checks for leaks)
88        using CppAD::thread_alloc; // the allocator
89        size_t capacity;           // capacity of an allocation
90
91        // tape identifier
92        int tag  = 0;
93        // AD domain space vector
94        ADVector a_x = thread_alloc::create_array<ADScalar>(n, capacity);
95        // AD range space vector
96        ADVector a_y = thread_alloc::create_array<ADScalar>(m, capacity);
97        // argument value in double
98        DblVector x = thread_alloc::create_array<double>(n, capacity);
99        // function value in double
100        DblVector y = thread_alloc::create_array<double>(m, capacity);
101
102       
103        // options that control sparse_jac
104        int        options[4];
105        extern bool global_boolsparsity;
106        if( global_boolsparsity )
107                options[0] = 1;  // sparsity by propagation of bit pattern
108        else
109                options[0] = 0;  // sparsity pattern by index domains
110        options[1] = 0; // (0 = safe mode, 1 = tight mode)
111        options[2] = 0; // see changing to -1 and back to 0 below
112        options[3] = 0; // (0 = column compression, 1 = row compression)
113
114        // structure that holds some of the work done by sparse_jac
115        int        nnz;                   // number of non-zero values
116        SizeVector rind   = CPPAD_NULL;   // row indices
117        SizeVector cind   = CPPAD_NULL;   // column indices
118        DblVector  values = CPPAD_NULL;   // Jacobian values
119
120        // initialize all entries as zero
121        for(i = 0; i < m; i++)
122        {       for(j = 0; j < n; j++)
123                        jacobian[ i * n + j ] = 0.;
124        }
125
126        // choose a value for x
127        CppAD::uniform_01(n, x);
128
129        // declare independent variables
130        int keep = 0; // keep forward mode results
131        trace_on(tag, keep);
132        for(j = 0; j < n; j++)
133                a_x[j] <<= x[j];
134
135        // AD computation of f (x)
136        CppAD::sparse_jac_fun<ADScalar>(m, n, a_x, row, col, order, a_y);
137
138        // create function object f : x -> y
139        for(i = 0; i < m; i++)
140                a_y[i] >>= y[i];
141        trace_off();
142
143        // Retrieve n_sweep using undocumented feature of sparsedrivers.cpp
144        int same_pattern = 0;
145        options[2]       = -1;
146        n_sweep = sparse_jac(tag, int(m), int(n), 
147                same_pattern, x, &nnz, &rind, &cind, &values, options
148        );
149        options[2]       = 0;
150        // ----------------------------------------------------------------------
151        if( ! global_onetape ) while(repeat--)
152        {       // choose a value for x
153                CppAD::uniform_01(n, x);
154
155                // declare independent variables
156                trace_on(tag, keep);
157                for(j = 0; j < n; j++)
158                        a_x[j] <<= x[j];
159
160                // AD computation of f (x)
161                CppAD::sparse_jac_fun<ADScalar>(m, n, a_x, row, col, order, a_y);
162
163                // create function object f : x -> y
164                for(i = 0; i < m; i++)
165                        a_y[i] >>= y[i];
166                trace_off();
167
168                // is this a repeat call with the same sparsity pattern
169                int same_pattern = 0;
170
171                // calculate the jacobian at this x
172                rind   = CPPAD_NULL;
173                cind   = CPPAD_NULL;
174                values = CPPAD_NULL;
175                sparse_jac(tag, int(m), int(n), 
176                        same_pattern, x, &nnz, &rind, &cind, &values, options
177                );
178                int int_n = int(n);
179                for(int k = 0; k < nnz; k++)
180                        jacobian[ rind[k] * int_n + cind[k] ] = values[k];
181
182                // free raw memory allocated by sparse_jac
183                free(rind);
184                free(cind);
185                free(values);
186        }
187        else
188        {       while(repeat--)
189                {       // choose a value for x
190                        CppAD::uniform_01(n, x);
191
192                        // calculate the jacobian at this x
193                        sparse_jac(tag, int(m), int(n), 
194                                same_pattern, x, &nnz, &rind, &cind, &values, options
195                        );
196                        same_pattern = 1;
197                }
198                int int_n = int(n);
199                for(int k = 0; k < nnz; k++)
200                        jacobian[ rind[k] * int_n + cind[k] ] = values[k];
201
202                // free raw memory allocated by sparse_jac
203                free(rind);
204                free(cind);
205                free(values);
206        }
207        // --------------------------------------------------------------------
208        // return argument
209        for(j = 0; j < n; j++)
210                x_return[j] = x[j];
211
212        // tear down
213        thread_alloc::delete_array(a_x);
214        thread_alloc::delete_array(a_y);
215        thread_alloc::delete_array(x);
216        thread_alloc::delete_array(y);
217        return true;
218}
219/* $$
220$end
221*/
Note: See TracBrowser for help on using the repository browser.