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

Last change on this file since 3320 was 3320, checked in by bradbell, 6 years ago
  1. g++ 4.8.2 has shadow warnings by default, but eigen and fadbad do not

these warnings, so supress then in these cases.

  1. Move check that arguments come before result into on place,

CPPAD_ASSERT_ARG_BEFORE_RESULT (only one argument case so far).

main.cpp: fix shadowing of index variable.
CMakeLists.txt: adapt to change in teuchos library name.
sparse_jacobian.cpp: fix a shadowed variable.
check_svn_id.sh: ignore svn_commit.sh.
gpl_license.sh: ignore svn_commit.sh.

  • Property svn:keywords set to Id
File size: 6.0 KB
Line 
1/* $Id: sparse_jacobian.cpp 3320 2014-09-11 23:06:21Z 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                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.