source: trunk/cppad/core/sparse_hes.hpp @ 3941

Last change on this file since 3941 was 3941, checked in by bradbell, 2 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: b7056a15a1028d7be587b70a3ecc44b1f42dc05e
end hash code: c8c4cc081accff3628e7e66370ec01e4c99afe8d

commit c8c4cc081accff3628e7e66370ec01e4c99afe8d
Author: Brad Bell <bradbell@…>
Date: Thu Jun 1 23:16:39 2017 -0600

Changes automatically generated by the autotools.

commit f4392bc3eee8f6d0ccd45a0bb3be51181e211680
Author: Brad Bell <bradbell@…>
Date: Thu Jun 1 23:11:56 2017 -0600

  1. Add colpack_jac.cpp example (rename colpack_jac.cpp->colpack_jacobian.cpp).
  2. Add colpack_hescpp example (rename colpack_hes.cpp->colpack_hessian.cpp).


test_one.sh.in: adapt to using test_boolofvoid for testing.
sparse_hes.hpp: fix bug in cppad.symmetric case.

commit 086b8a8709b0c9cb01ce2cf8bc7910e903105ff7
Author: Brad Bell <bradbell@…>
Date: Thu Jun 1 08:54:59 2017 -0600

  1. Fix bug in use of colpack (see kludge in comments).
  2. Fix colpack.symmetric (not general) and add colpack.general.
  3. Deprecate colpack.star.
  4. More autotools from install to deprecated.
  5. Advance to cppad-20170601.

commit 23f26c060648f5c6fc62a1598c659aeccc5ca46f
Author: Brad Bell <bradbell@…>
Date: Tue May 30 08:14:04 2017 -0700

Advance to cppad-20170530.

commit 97f8c08509865d1bfb7ec2e5cd557ddc979f8412
Author: Brad Bell <bradbell@…>
Date: Tue May 30 07:38:47 2017 -0700

debug_rel branch:
There is a problem with speed sparse_hessian debug that goes back to master.
Supresss debug in cppad speed tests until it is fixed.

commit 39ea0d7d9c041784ccd26ce80d19a7ab02752818
Author: Brad Bell <bradbell@…>
Date: Mon May 29 22:34:22 2017 -0700

debug_rel branch:
run_cmake.sh: fix debug_none case.
CMakeLists.txt: use cppad_debug_which to determine debug or release.
CMakeLists.txt: let set_compile_flags determkine build type.

commit 191553e54dca407207789cf0d7c6c27fe6188775
Author: Brad Bell <bradbell@…>
Date: Mon May 29 19:53:08 2017 -0700

debug_rel branch:
Use set_compile_flags in introduction.

commit fba276a84e58d9a0d0944168d5706b7446beb32c
Author: Brad Bell <bradbell@…>
Date: Mon May 29 19:46:30 2017 -0700

debug_rel branch:
Use set_compile_flags in eample/multi_thread subdirectories.

commit 66c8cdb266fa3af29b211b8c870a3aed7a13b021
Author: Brad Bell <bradbell@…>
Date: Mon May 29 18:56:48 2017 -0700

debug_rel branch:
Use set_compile_flags in speed directory.

commit c431b15ee7714d3106234bc527ba2f9a836750e1
Author: Brad Bell <bradbell@…>
Date: Mon May 29 18:36:51 2017 -0700

debug_rel branch:
Convert cppad_ipopt to use set_compile_flags and cppad_debug_which.


CMakeLists.txt: alwasy compile for release to reduce testing time.

commit 2c95b0019f1b665fb14b9f00b049e8b5fb11f89d
Author: Brad Bell <bradbell@…>
Date: Mon May 29 16:55:07 2017 -0700

debug_rel branch:
Add cppad_debug_which to the cmake command line.

commit fd8d1498cf6dc092deca41f764cbb2a001a4cc88
Author: Brad Bell <bradbell@…>
Date: Mon May 29 08:14:23 2017 -0700

debug_rel branch:
Change random_debug_release -> set_compile_flags.

commit 159f5a5aa09012213a52f4ed1c9f0607129a5fe7
Author: Brad Bell <bradbell@…>
Date: Mon May 29 06:50:43 2017 -0700

debug_rel branch:
Update the autotools automatically generated build files.


batch_edit.sh: Start comments about a plan for editing all the source files.
get_sacado.sh: advance to trilions-12.10.11.
makefile.am: advance to trilinos-12.10.1

commit 302153317cd296ec6f927c3202cf96bf38594bbb
Author: Brad Bell <bradbell@…>
Date: Mon May 29 05:20:00 2017 -0700

debug_rel branch:
Add error message if sacado configuration file does not exist.

commit 3f01a631ae808c3a1359e53e1cd55e9a0ea88711
Author: Brad Bell <bradbell@…>
Date: Mon May 29 04:24:00 2017 -0700

debug_rel branch:
CMakeLists.txt: automate naming of libraries Sacado needs.
checkpoint.cpp: fix warnings.

commit dd240928c0c8b6972a8197c985ccc01f08b8886b
Author: Brad Bell <bradbell@…>
Date: Sun May 28 08:25:20 2017 -0700

debug_rel branch:
sparse_sub_hes.cpp: add missing cases found by clang compiler.

commit 30a0c35f1ac50ec425be9a2b7b026284026eccd7
Author: Brad Bell <bradbell@…>
Date: Sun May 28 07:57:36 2017 -0700

debug_rel branch:
eigen_cholesky.hpp: fix compiler warning.
harmonic_time.cpp: remove include that is not used.
forward_active.cpp: fix compiler warning.

commit 4876d14e49dc235865b1574fb38a55cf5ea7a422
Author: Brad Bell <bradbell@…>
Date: Sun May 28 06:19:48 2017 -0700

debug_rel branch:
random_debug_release.cmake: fix comment, remove message replaced by random_choice_0123 in output.
multiple_solution.cpp: fix warnings with clang compiler.
eigen_cholesky.hpp: fix warnings with clang compiler.
compare_change.cpp: fix CPPAD_DEBUG_AND_RELEASE case.

commit 2c51a18f35188d04d2f94069382439580e23f4ac
Author: Brad Bell <bradbell@…>
Date: Sat May 27 21:04:37 2017 -0700

debug_rel branch:
Advance version to cppad-20170527.

commit 4500887b362537774b05e954ad2a95b65a7b8ba0
Author: Brad Bell <bradbell@…>
Date: Sat May 27 09:04:56 2017 -0700

debug_rel branch:
Ramdomly select debug or release flags in example directory.


CMakeLists.txt: always debug for multi_threed examples.

commit 140b5269a0b1a30643894e5a7a8c9a5eb1310301
Author: Brad Bell <bradbell@…>
Date: Sat May 27 08:10:51 2017 -0700

debug_rel branch:
Changing how we set all debug and release flags.

commit e6fb2639db1288fb75de4030b5906df1e41756f9
Author: Brad Bell <bradbell@…>
Date: Sat May 27 07:30:24 2017 -0700

debug_rel branch:
Replace use of cppad_extra_debug by CPPAD_DEBUG_AND_RELEASE.

commit fbbfd0f6e94862174a8a7a17308489ffddb28084
Author: Brad Bell <bradbell@…>
Date: Sat May 27 05:55:58 2017 -0700

debug_rel branch:
Improve random selection of which files are build for release or debug.


forward.cpp: use new -DCPPAD_DEBUG_AND_RELEASE flag.

commit 284be366fb5e2f685a0c71ea6a0e3f74584bf187
Author: Brad Bell <bradbell@…>
Date: Thu May 25 07:39:32 2017 -0700

debug_rel branch:
Add test that failed before change to player.


player.hpp: Fix so it has the same size in debug and release more.
checkpoint.cpp: fix warning when compiling for release.
run_cmake.sh: prepare to use random number to switch debug and release set.
CMakeLists.txt: switch to only test debug (for now).

commit f32375b77e3825628fee6cb160f691a32c48b796
Author: Brad Bell <bradbell@…>
Date: Wed May 24 12:04:27 2017 -0700

debug_rel branch:
forward.cpp: fix a warning during release build.

commit 5fcc7eb78ae8de9f1dbc6c4f0c76fe38e8aeba95
Author: Brad Bell <bradbell@…>
Date: Wed May 24 10:11:12 2017 -0700

debug_rel branch:
CMakeLists.txt: make easy to mix debug and release builds.
eigen_mat_inv.hpp: fix release version warning.

commit 696266f3d62079f5e3bfb1a0f60a7e4f8134e068
Author: Brad Bell <bradbell@…>
Date: Wed May 24 05:43:29 2017 -0700

push_git2svn.py: user ./build in place of ./build/work.
testvector.hpp: improve comments about replacing CPPAD_TESTVECTOR.

File size: 15.4 KB
Line 
1# ifndef CPPAD_CORE_SPARSE_HES_HPP
2# define CPPAD_CORE_SPARSE_HES_HPP
3/* --------------------------------------------------------------------------
4CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
5
6CppAD is distributed under multiple licenses. This distribution is under
7the terms of the
8                    Eclipse Public License Version 1.0.
9
10A copy of this license is included in the COPYING file of this distribution.
11Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
12-------------------------------------------------------------------------- */
13
14/*
15$begin sparse_hes$$
16$spell
17        const
18        Taylor
19        rc
20        rcv
21        nr
22        nc
23        hes
24        std
25        cppad
26        colpack
27        cmake
28        Jacobian
29$$
30
31$section Computing Sparse Hessians$$
32
33$head Syntax$$
34$icode%n_sweep% = %f%.sparse_hes(
35        %x%, %w%, %subset%, %pattern%, %coloring%, %work%
36)%$$
37
38$head Purpose$$
39We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
40function corresponding to $icode f$$.
41Here $icode n$$ is the $cref/domain/seq_property/Domain/$$ size,
42and $icode m$$ is the $cref/range/seq_property/Range/$$ size, or $icode f$$.
43The syntax above takes advantage of sparsity when computing the Hessian
44$latex \[
45        H(x) = \dpow{2}{x} \sum_{i=0}^{m-1} w_i F_i (x)
46\] $$
47In the sparse case, this should be faster and take less memory than
48$cref Hessian$$.
49The matrix element $latex H_{i,j} (x)$$ is the second partial of
50$latex w^\R{T} F (x)$$ with respect to $latex x_i$$ and $latex x_j$$.
51
52$head SizeVector$$
53The type $icode SizeVector$$ is a $cref SimpleVector$$ class with
54$cref/elements of type/SimpleVector/Elements of Specified Type/$$
55$code size_t$$.
56
57$head BaseVector$$
58The type $icode BaseVector$$ is a $cref SimpleVector$$ class with
59$cref/elements of type/SimpleVector/Elements of Specified Type/$$
60$code size_t$$.
61
62$head f$$
63This object has prototype
64$codei%
65        ADFun<%Base%> %f%
66%$$
67Note that the Taylor coefficients stored in $icode f$$ are affected
68by this operation; see
69$cref/uses forward/sparse_hes/Uses Forward/$$ below.
70
71$head x$$
72This argument has prototype
73$codei%
74        const %BaseVector%& %x%
75%$$
76and its size is $icode n$$.
77It specifies the point at which to evaluate the Hessian
78$latex H(x)$$.
79
80$head w$$
81This argument has prototype
82$codei%
83        const %BaseVector%& %w%
84%$$
85and its size is $icode m$$.
86It specifies the weight for each of the components of $latex F(x)$$;
87i.e. $latex w_i$$ is the weight for $latex F_i (x)$$.
88
89$head subset$$
90This argument has prototype
91$codei%
92        sparse_rcv<%SizeVector%, %BaseVector%>& %subset%
93%$$
94Its row size and column size is $icode n$$; i.e.,
95$icode%subset%.nr() == %n%$$ and $icode%subset%.nc() == %n%$$.
96It specifies which elements of the Hessian are computed.
97$list number$$
98The input value of its value vector
99$icode%subset%.val()%$$ does not matter.
100Upon return it contains the value of the corresponding elements
101of the Hessian.
102$lnext
103All of the row, column pairs in $icode subset$$ must also appear in
104$icode pattern$$; i.e., they must be possibly non-zero.
105$lnext
106The Hessian is symmetric, so one has a choice as to which off diagonal
107elements to put in $icode subset$$.
108It will probably be more efficient if one makes this choice so that
109the there are more entries in each non-zero column of $icode subset$$;
110see $cref/n_sweep/sparse_hes/n_sweep/$$ below.
111$lend
112
113$head pattern$$
114This argument has prototype
115$codei%
116        const sparse_rc<%SizeVector%>& %pattern%
117%$$
118Its row size and column size is $icode n$$; i.e.,
119$icode%pattern%.nr() == %n%$$ and $icode%pattern%.nc() == %n%$$.
120It is a sparsity pattern for the Hessian $latex H(x)$$.
121If the $th i$$ row ($th j$$ column) does not appear in $icode subset$$,
122the $th i$$ row ($th j$$ column) of $icode pattern$$ does not matter
123and need not be computed.
124This argument is not used (and need not satisfy any conditions),
125when $cref/work/sparse_hes/work/$$ is non-empty.
126
127$subhead subset$$
128If the $th i$$ row and $th i$$ column do not appear in $icode subset$$,
129the $th i$$ row and column of $icode pattern$$ do not matter.
130In this case the $th i-th$$ row and column may have no entries in
131$icode pattern$$ even though they are possibly non-zero in $latex H(x)$$.
132(This can be used to reduce the amount of computation required to find
133$icode pattern$$.)
134
135$head coloring$$
136The coloring algorithm determines which rows and columns
137can be computed during the same sweep.
138This field has prototype
139$codei%
140        const std::string& %coloring%
141%$$
142This value only matters when work is empty; i.e.,
143after the $icode work$$ constructor or $icode%work%.clear()%$$.
144
145$subhead cppad.symmetric$$
146This coloring takes advantage of the fact that the Hessian matrix
147is symmetric when find a coloring that requires fewer
148$cref/sweeps/sparse_hes/n_sweep/$$.
149
150$subhead cppad.general$$
151This is the same as the sparse Jacobian
152$cref/cppad/sparse_jac/coloring/cppad/$$ method
153which does not take advantage of symmetry.
154
155$subhead colpack.symmetric$$
156If $cref colpack_prefix$$ was specified on the
157$cref/cmake command/cmake/CMake Command/$$ line,
158you can set $icode coloring$$ to $code colpack.symmetric$$.
159This also takes advantage of the fact that the Hessian matrix is symmetric.
160
161$subhead colpack.general$$
162If $cref colpack_prefix$$ was specified on the
163$cref/cmake command/cmake/CMake Command/$$ line,
164you can set $icode coloring$$ to $code colpack.general$$.
165This is the same as the sparse Jacobian
166$cref/colpack/sparse_jac/coloring/colpack/$$ method
167which does not take advantage of symmetry.
168
169$subhead colpack.star Deprecated 2017-06-01$$
170The $code colpack.star$$ method is deprecated.
171It is the same as the $code colpack.symmetric$$ method
172which should be used instead.
173
174
175$head work$$
176This argument has prototype
177$codei%
178        sparse_hes_work& %work%
179%$$
180We refer to its initial value,
181and its value after $icode%work%.clear()%$$, as empty.
182If it is empty, information is stored in $icode work$$.
183This can be used to reduce computation when
184a future call is for the same object $icode f$$,
185and the same subset of the Hessian.
186If either of these values change, use $icode%work%.clear()%$$ to
187empty this structure.
188
189$head n_sweep$$
190The return value $icode n_sweep$$ has prototype
191$codei%
192        size_t %n_sweep%
193%$$
194It is the number of first order forward sweeps
195used to compute the requested Hessian values.
196Each first forward sweep is followed by a second order reverse sweep
197so it is also the number of reverse sweeps.
198It is also the number of colors determined by the coloring method
199mentioned above.
200This is proportional to the total computational work,
201not counting the zero order forward sweep,
202or combining multiple columns and rows into a single sweep.
203
204$head Uses Forward$$
205After each call to $cref Forward$$,
206the object $icode f$$ contains the corresponding
207$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
208After a call to $code sparse_hes$$
209the zero order coefficients correspond to
210$codei%
211        %f%.Forward(0, %x%)
212%$$
213All the other forward mode coefficients are unspecified.
214
215$head Example$$
216$children%
217        example/sparse/sparse_hes.cpp
218%$$
219The files $cref sparse_hes.cpp$$
220is an example and test of $code sparse_hes$$.
221It returns $code true$$, if it succeeds, and $code false$$ otherwise.
222
223$head Subset Hessian$$
224The routine
225$cref sparse_sub_hes.cpp$$
226is an example and test that compute a subset of a sparse Hessian.
227It returns $code true$$, for success, and $code false$$ otherwise.
228
229$end
230*/
231# include <cppad/core/cppad_assert.hpp>
232# include <cppad/local/sparse_internal.hpp>
233# include <cppad/local/color_general.hpp>
234# include <cppad/local/color_symmetric.hpp>
235
236/*!
237\file sparse_hes.hpp
238Sparse Hessian calculation routines.
239*/
240namespace CppAD { // BEGIN_CPPAD_NAMESPACE
241
242/*!
243Class used to hold information used by Sparse Hessian routine in this file,
244so it does not need to be recomputed every time.
245*/
246class sparse_hes_work {
247        public:
248                /// row and column indicies for return values
249                /// (some may be reflected by symmetric coloring algorithms)
250                CppAD::vector<size_t> row;
251                CppAD::vector<size_t> col;
252                /// indices that sort the row and col arrays by color
253                CppAD::vector<size_t> order;
254                /// results of the coloring algorithm
255                CppAD::vector<size_t> color;
256
257                /// constructor
258                sparse_hes_work(void)
259                { }
260                /// inform CppAD that this information needs to be recomputed
261                void clear(void)
262                {
263                        row.clear();
264                        col.clear();
265                        order.clear();
266                        color.clear();
267                }
268};
269// ----------------------------------------------------------------------------
270/*!
271Calculate sparse Hessians using forward mode
272
273\tparam Base
274the base type for the recording that is stored in the ADFun object.
275
276\tparam SizeVector
277a simple vector class with elements of type size_t.
278
279\tparam BaseVector
280a simple vector class with elements of type Base.
281
282\param x
283a vector of length n, the number of independent variables in f
284(this ADFun object).
285
286\param w
287a vector of length m, the number of dependent variables in f
288(this ADFun object).
289
290\param subset
291specifices the subset of the sparsity pattern where the Hessian is evaluated.
292subset.nr() == n,
293subset.nc() == n.
294
295\param pattern
296is a sparsity pattern for the Hessian of w^T * f;
297pattern.nr() == n,
298pattern.nc() == n,
299where m is number of dependent variables in f.
300
301\param coloring
302determines which coloring algorithm is used.
303This must be cppad.symmetric, cppad.general, colpack.symmetic,
304or colpack.star.
305
306\param work
307this structure must be empty, or contain the information stored
308by a previous call to sparse_hes.
309The previous call must be for the same ADFun object f
310and the same subset.
311
312\return
313This is the number of first order forward
314(and second order reverse) sweeps used to compute thhe Hessian.
315*/
316template <class Base>
317template <class SizeVector, class BaseVector>
318size_t ADFun<Base>::sparse_hes(
319        const BaseVector&                    x        ,
320        const BaseVector&                    w        ,
321        sparse_rcv<SizeVector , BaseVector>& subset   ,
322        const sparse_rc<SizeVector>&         pattern  ,
323        const std::string&                   coloring ,
324        sparse_hes_work&                     work     )
325{       size_t n = Domain();
326        //
327        CPPAD_ASSERT_KNOWN(
328                subset.nr() == n,
329                "sparse_hes: subset.nr() not equal domain dimension for f"
330        );
331        CPPAD_ASSERT_KNOWN(
332                subset.nc() == n,
333                "sparse_hes: subset.nc() not equal domain dimension for f"
334        );
335        CPPAD_ASSERT_KNOWN(
336                size_t( x.size() ) == n,
337                "sparse_hes: x.size() not equal domain dimension for f"
338        );
339        CPPAD_ASSERT_KNOWN(
340                size_t( w.size() ) == Range(),
341                "sparse_hes: w.size() not equal range dimension for f"
342        );
343        //
344        // work information
345        vector<size_t>& row(work.row);
346        vector<size_t>& col(work.col);
347        vector<size_t>& color(work.color);
348        vector<size_t>& order(work.order);
349        //
350        // subset information
351        const SizeVector& subset_row( subset.row() );
352        const SizeVector& subset_col( subset.col() );
353        //
354        // point at which we are evaluationg the Hessian
355        Forward(0, x);
356        //
357        // number of elements in the subset
358        size_t K = subset.nnz();
359        //
360        // check for case were there is nothing to do
361        // (except for call to Forward(0, x)
362        if( K == 0 )
363                return 0;
364        //
365# ifndef NDEBUG
366        if( color.size() != 0 )
367        {       CPPAD_ASSERT_KNOWN(
368                        color.size() == n,
369                        "sparse_hes: work is non-empty and conditions have changed"
370                );
371                CPPAD_ASSERT_KNOWN(
372                        row.size() == K,
373                        "sparse_hes: work is non-empty and conditions have changed"
374                );
375                CPPAD_ASSERT_KNOWN(
376                        col.size() == K,
377                        "sparse_hes: work is non-empty and conditions have changed"
378                );
379                //
380                for(size_t k = 0; k < K; k++)
381                {       bool ok = row[k] == subset_row[k] && col[k] == subset_col[k];
382                        ok     |= row[k] == subset_col[k] && col[k] == subset_row[k];
383                        CPPAD_ASSERT_KNOWN(
384                                ok,
385                                "sparse_hes: work is non-empty and conditions have changed"
386                        );
387                }
388        }
389# endif
390        //
391        // check for case where input work is empty
392        if( color.size() == 0 )
393        {       // compute work color and order vectors
394                CPPAD_ASSERT_KNOWN(
395                        pattern.nr() == n,
396                        "sparse_hes: pattern.nr() not equal domain dimension for f"
397                );
398                CPPAD_ASSERT_KNOWN(
399                        pattern.nc() == n,
400                        "sparse_hes: pattern.nc() not equal domain dimension for f"
401                );
402                //
403                // initialize work row, col to be same as subset row, col
404                row.resize(K);
405                col.resize(K);
406                // cannot assign vectors becasue may be of different types
407                // (SizeVector and CppAD::vector<size_t>)
408                for(size_t k = 0; k < K; k++)
409                {       row[k] = subset_row[k];
410                        col[k] = subset_col[k];
411                }
412                //
413                // convert pattern to an internal version of its transpose
414                vector<size_t> internal_index(n);
415                for(size_t j = 0; j < n; j++)
416                        internal_index[j] = j;
417                bool transpose   = true;
418                bool zero_empty  = false;
419                bool input_empty = true;
420                local::sparse_list internal_pattern;
421                internal_pattern.resize(n, n);
422                local::set_internal_sparsity(zero_empty, input_empty,
423                        transpose, internal_index, internal_pattern, pattern
424                );
425                //
426                // execute coloring algorithm
427                // (we are using transpose becasue coloring groups rows, not columns)
428                color.resize(n);
429                if( coloring == "cppad.general" )
430                        local::color_general_cppad(internal_pattern, col, row, color);
431                else if( coloring == "cppad.symmetric" )
432                        local::color_symmetric_cppad(internal_pattern, col, row, color);
433                else if( coloring == "colpack.general" )
434                {
435# if CPPAD_HAS_COLPACK
436                        local::color_general_colpack(internal_pattern, col, row, color);
437# else
438                        CPPAD_ASSERT_KNOWN(
439                                false,
440                                "sparse_hes: coloring = colpack.star "
441                                "and colpack_prefix not in cmake command line."
442                        );
443# endif
444                }
445                else if(
446                        coloring == "colpack.symmetric" ||
447                        coloring == "colpack.star"
448                )
449                {
450# if CPPAD_HAS_COLPACK
451                        local::color_symmetric_colpack(internal_pattern, col, row, color);
452# else
453                        CPPAD_ASSERT_KNOWN(
454                                false,
455                                "sparse_hes: coloring = colpack.symmetic or colpack.star "
456                                "and colpack_prefix not in cmake command line."
457                        );
458# endif
459                }
460                else CPPAD_ASSERT_KNOWN(
461                        false,
462                        "sparse_hes: coloring is not valid."
463                );
464                //
465                // put sorting indices in color order
466                SizeVector key(K);
467                order.resize(K);
468                for(size_t k = 0; k < K; k++)
469                        key[k] = color[ col[k] ];
470                index_sort(key, order);
471        }
472        // Base versions of zero and one
473        Base one(1.0);
474        Base zero(0.0);
475        //
476        size_t n_color = 1;
477        for(size_t j = 0; j < n; j++) if( color[j] < n )
478                n_color = std::max(n_color, color[j] + 1);
479        //
480        // initialize the return Hessian values as zero
481        for(size_t k = 0; k < K; k++)
482                subset.set(k, zero);
483        //
484        // direction vector for calls to first order forward
485        BaseVector dx(n);
486        //
487        // return values for calls to second order reverse
488        BaseVector ddw(2 * n);
489        //
490        // loop over colors
491        size_t k = 0;
492        for(size_t ell = 0; ell < n_color; ell++)
493        if( k  == K )
494        {       // kludge because colpack returns colors that are not used
495                // (it does not know about the subset corresponding to row, col)
496                CPPAD_ASSERT_UNKNOWN(
497                        coloring == "colpack.general" ||
498                        coloring == "colpack.symmetric" ||
499                        coloring == "colpack.star"
500                );
501        }
502        else if( color[ col[ order[k] ] ] != ell )
503        {       // kludge because colpack returns colors that are not used
504                // (it does not know about the subset corresponding to row, col)
505                CPPAD_ASSERT_UNKNOWN(
506                        coloring == "colpack.general" ||
507                        coloring == "colpack.symmetic" ||
508                        coloring == "colpack.star"
509                );
510        }
511        else
512        {       CPPAD_ASSERT_UNKNOWN( color[ col[ order[k] ] ] == ell );
513                //
514                // combine all columns with this color
515                for(size_t j = 0; j < n; j++)
516                {       dx[j] = zero;
517                        if( color[j] == ell )
518                                dx[j] = one;
519                }
520                // call forward mode for all these rows at once
521                Forward(1, dx);
522                //
523                // evaluate derivative of w^T * F'(x) * dx
524                ddw = Reverse(2, w);
525                //
526                // set the corresponding components of the result
527                while( k < K && color[ col[order[k]] ] == ell )
528                {       size_t index = row[ order[k] ] * 2 + 1;
529                        subset.set(order[k], ddw[index] );
530                        k++;
531                }
532        }
533        // check that all the required entries have been set
534        CPPAD_ASSERT_UNKNOWN( k == K );
535        return n_color;
536}
537
538} // END_CPPAD_NAMESPACE
539
540# endif
Note: See TracBrowser for help on using the repository browser.