source: trunk/cppad/local/rev_sparse_jac.hpp @ 3746

Last change on this file since 3746 was 3746, checked in by bradbell, 4 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 57d3215cae5e9da7b4c92e89f038d70870ef7327
end hash code: 9aebc1ca2c0949dec7c2d156517db26e60f28159

commit 9aebc1ca2c0949dec7c2d156517db26e60f28159
Author: Brad Bell <bradbell@…>
Date: Sun Nov 8 20:15:38 2015 -0800

Remove invisible white space.

commit a92ac50e9f4c8d0007ea5a245b3e23145dfcebfe
Author: Brad Bell <bradbell@…>
Date: Sun Nov 8 20:15:31 2015 -0800

Use vectorBool with partial sparsity patterns per pass to reduce memory requirements.


solve_callback.hpp: remove invisible white space.
rev_sparse_jac.hpp: fix bug (argument transposed).
bool_sparsity.cpp: remove invisible white space.

commit c09744b13ba2c70d6ffa857206d45560154d800a
Author: Brad Bell <bradbell@…>
Date: Sun Nov 8 03:22:57 2015 -0800

check_for_nan.hpp: fix minor type in documentation.

  • Property svn:keywords set to Id
File size: 19.8 KB
Line 
1/* $Id: rev_sparse_jac.hpp 3746 2015-11-09 04:50:27Z bradbell $ */
2# ifndef CPPAD_REV_SPARSE_JAC_INCLUDED
3# define CPPAD_REV_SPARSE_JAC_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin RevSparseJac$$
18$spell
19        optimizer
20        nz
21        CondExpRel
22        std
23        VecAD
24        var
25        Jacobian
26        Jac
27        const
28        Bool
29        Dep
30        proportional
31$$
32
33$section Jacobian Sparsity Pattern: Reverse Mode$$
34
35$index RevSparseJac$$
36$index reverse, sparse Jacobian$$
37$index sparse, reverse Jacobian$$
38$index pattern, reverse Jacobian$$
39
40$head Syntax$$
41$icode%s% = %f%.RevSparseJac(%q%, %r%)
42%$$
43$icode%s% = %f%.RevSparseJac(%q%, %r%, %transpose%, %dependency%)%$$
44
45$head Purpose$$
46We use $latex F : B^n \rightarrow B^m$$ to denote the
47$cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$.
48For a fixed matrix $latex R \in B^{q \times m}$$,
49the Jacobian of $latex R * F( x )$$
50with respect to $latex x$$ is
51$latex \[
52        S(x) = R * F^{(1)} ( x )
53\] $$
54Given a
55$cref/sparsity pattern/glossary/Sparsity Pattern/$$
56for $latex R$$,
57$code RevSparseJac$$ returns a sparsity pattern for the $latex S(x)$$.
58
59$head f$$
60The object $icode f$$ has prototype
61$codei%
62        ADFun<%Base%> %f%
63%$$
64
65$head x$$
66the sparsity pattern is valid for all values of the independent
67variables in $latex x \in B^n$$
68(even if it has $cref CondExp$$ or $cref VecAD$$ operations).
69
70$head q$$
71The argument $icode q$$ has prototype
72$codei%
73        size_t %q%
74%$$
75It specifies the number of rows in
76$latex R \in B^{q \times m}$$ and the
77Jacobian $latex S(x) \in B^{q \times n}$$.
78
79$head transpose$$
80The argument $icode transpose$$ has prototype
81$codei%
82        bool %transpose%
83%$$
84The default value $code false$$ is used when $icode transpose$$ is not present.
85
86$head dependency$$
87The argument $icode dependency$$ has prototype
88$codei%
89        bool %dependency%
90%$$
91If $icode dependency$$ is true,
92the $cref/dependency pattern/dependency.cpp/Dependency Pattern/$$
93(instead of sparsity pattern) is computed.
94
95$head r$$
96The argument $icode s$$ has prototype
97$codei%
98        const %VectorSet%& %r%
99%$$
100see $cref/VectorSet/RevSparseJac/VectorSet/$$ below.
101
102$subhead transpose false$$
103If $icode r$$ has elements of type $code bool$$,
104its size is $latex q * m$$.
105If it has elements of type $code std::set<size_t>$$,
106its size is $icode q$$ and all its set elements are between
107zero and $latex m - 1$$.
108It specifies a
109$cref/sparsity pattern/glossary/Sparsity Pattern/$$
110for the matrix $latex R \in B^{q \times m}$$.
111
112$subhead transpose true$$
113If $icode r$$ has elements of type $code bool$$,
114its size is $latex m * q$$.
115If it has elements of type $code std::set<size_t>$$,
116its size is $icode m$$ and all its set elements are between
117zero and $latex q - 1$$.
118It specifies a
119$cref/sparsity pattern/glossary/Sparsity Pattern/$$
120for the matrix $latex R^\R{T} \in B^{m \times q}$$.
121
122$head s$$
123The return value $icode s$$ has prototype
124$codei%
125        %VectorSet% %s%
126%$$
127see $cref/VectorSet/RevSparseJac/VectorSet/$$ below.
128
129$subhead transpose false$$
130If it has elements of type $code bool$$,
131its size is $latex q * n$$.
132If it has elements of type $code std::set<size_t>$$,
133its size is $icode q$$ and all its set elements are between
134zero and $latex n - 1$$.
135It specifies a
136$cref/sparsity pattern/glossary/Sparsity Pattern/$$
137for the matrix $latex S(x) \in {q \times n}$$.
138
139$subhead transpose true$$
140If it has elements of type $code bool$$,
141its size is $latex n * q$$.
142If it has elements of type $code std::set<size_t>$$,
143its size is $icode n$$ and all its set elements are between
144zero and $latex q - 1$$.
145It specifies a
146$cref/sparsity pattern/glossary/Sparsity Pattern/$$
147for the matrix $latex S(x)^\R{T} \in {n \times q}$$.
148
149$head VectorSet$$
150The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
151$cref/elements of type/SimpleVector/Elements of Specified Type/$$
152$code bool$$ or $code std::set<size_t>$$;
153see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
154of the difference.
155
156$head Entire Sparsity Pattern$$
157Suppose that $latex q = m$$ and
158$latex R$$ is the $latex m \times m$$ identity matrix.
159In this case,
160the corresponding value for $icode s$$ is a
161sparsity pattern for the Jacobian $latex S(x) = F^{(1)} ( x )$$.
162
163$head Example$$
164$children%
165        example/rev_sparse_jac.cpp
166%$$
167The file
168$cref rev_sparse_jac.cpp$$
169contains an example and test of this operation.
170It returns true if it succeeds and false otherwise.
171
172$end
173-----------------------------------------------------------------------------
174*/
175
176# include <cppad/local/std_set.hpp>
177
178namespace CppAD { // BEGIN_CPPAD_NAMESPACE
179/*!
180\file rev_sparse_jac.hpp
181Reverse mode Jacobian sparsity patterns.
182*/
183
184// =========================================================================
185// RevSparseJacBool
186/*!
187Calculate Jacobian vector of bools sparsity patterns using reverse mode.
188
189The C++ source code corresponding to this operation is
190\verbatim
191        s = f.RevSparseJac(q, r)
192\endverbatim
193
194\tparam Base
195is the base type for this recording.
196
197\tparam VectorSet
198is a simple vector class with elements of type \c bool.
199
200\param transpose
201are the sparsity patterns transposed.
202
203\param dependency
204Are the derivatives with respect to left and right of the expression below
205considered to be non-zero:
206\code
207        CondExpRel(left, right, if_true, if_false)
208\endcode
209This is used by the optimizer to obtain the correct dependency relations.
210
211\param q
212is the number of rows in the matrix \f$ R \f$.
213
214\param r
215is a sparsity pattern for the matrix \f$ R \f$.
216
217\param s
218the input value of \a s must be a vector with size <tt>p*n</tt>
219where \c n is the number of independent variables
220corresponding to the operation sequence stored in \a play.
221The input value of the components of \c s does not matter.
222On output, \a s is the sparsity pattern for the matrix
223\f[
224        S(x) = R * F^{(1)} (x)
225\f]
226where \f$ F \f$ is the function corresponding to the operation sequence
227and \a x is any argument value.
228
229\param total_num_var
230is the total number of variable in this recording.
231
232\param dep_taddr
233maps dependendent variable index
234to the corresponding variable in the tape.
235
236\param ind_taddr
237maps independent variable index
238to the corresponding variable in the tape.
239
240\param play
241is the recording that defines the function we are computing the sparsity
242pattern for.
243*/
244
245template <class Base, class VectorSet>
246void RevSparseJacBool(
247        bool                   transpose        ,
248        bool                   dependency       ,
249        size_t                 q                ,
250        const VectorSet&       r                ,
251        VectorSet&             s                ,
252        size_t                 total_num_var    ,
253        CppAD::vector<size_t>& dep_taddr        ,
254        CppAD::vector<size_t>& ind_taddr        ,
255        CppAD::player<Base>&   play             )
256{
257        // temporary indices
258        size_t i, j;
259
260        // check VectorSet is Simple Vector class with bool elements
261        CheckSimpleVector<bool, VectorSet>();
262
263        // range and domain dimensions for F
264        size_t m = dep_taddr.size();
265        size_t n = ind_taddr.size();
266
267        CPPAD_ASSERT_KNOWN(
268                q > 0,
269                "RevSparseJac: q is not greater than zero"
270        );
271        CPPAD_ASSERT_KNOWN(
272                size_t(r.size()) == q * m,
273                "RevSparseJac: size of r is not equal to\n"
274                "q times range dimension for ADFun object."
275        );
276
277        // vector of sets that will hold the results
278        sparse_pack    var_sparsity;
279        var_sparsity.resize(total_num_var, q);
280
281        // The sparsity pattern corresponding to the dependent variables
282        for(i = 0; i < m; i++)
283        {       CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
284                if( transpose )
285                {       for(j = 0; j < q; j++) if( r[ i * q + j ] )
286                                var_sparsity.add_element( dep_taddr[i], j );
287                }
288                else
289                {       for(j = 0; j < q; j++) if( r[ j * m + i ] )
290                                var_sparsity.add_element( dep_taddr[i], j );
291                }
292        }
293
294        // evaluate the sparsity patterns
295        RevJacSweep(
296                dependency,
297                n,
298                total_num_var,
299                &play,
300                var_sparsity
301        );
302
303        // return values corresponding to dependent variables
304        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q * n );
305        for(j = 0; j < n; j++)
306        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) );
307
308                // ind_taddr[j] is operator taddr for j-th independent variable
309                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp );
310
311                // extract the result from var_sparsity
312                if( transpose )
313                {       for(i = 0; i < q; i++)
314                                s[ j * q + i ] = false;
315                }
316                else
317                {       for(i = 0; i < q; i++)
318                                s[ i * n + j ] = false;
319                }
320                CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
321                var_sparsity.begin(j+1);
322                i = var_sparsity.next_element();
323                while( i < q )
324                {       if( transpose )
325                                s[ j * q + i ] = true;
326                        else    s[ i * n + j ] = true;
327                        i  = var_sparsity.next_element();
328                }
329        }
330}
331// =========================================================================
332// RevSparseJacSet
333/*!
334Calculate Jacobian vector of sets sparsity patterns using reverse mode.
335
336The C++ source code corresponding to this operation is
337\verbatim
338        s = f.RevSparseJac(q, r)
339\endverbatim
340
341\tparam Base
342see \c RevSparseJacBool.
343
344\tparam VectorSet
345is a simple vector class with elements of type \c std::set<size_t>.
346
347\param transpose
348see \c RevSparseJacBool.
349
350\param dependency
351see \c RevSparseJacBool.
352
353\param q
354see \c RevSparseJacBool.
355
356\param r
357see \c RevSparseJacBool.
358
359\param s
360see \c RevSparseJacBool.
361
362\param total_num_var
363see \c RevSparseJacBool.
364
365\param dep_taddr
366see \c RevSparseJacBool.
367
368\param ind_taddr
369see \c RevSparseJacBool.
370
371\param play
372see \c RevSparseJacBool.
373*/
374template <class Base, class VectorSet>
375void RevSparseJacSet(
376        bool                   transpose        ,
377        bool                   dependency       ,
378        size_t                 q                ,
379        const VectorSet&       r                ,
380        VectorSet&             s                ,
381        size_t                 total_num_var    ,
382        CppAD::vector<size_t>& dep_taddr        ,
383        CppAD::vector<size_t>& ind_taddr        ,
384        CppAD::player<Base>&   play             )
385{
386        // temporary indices
387        size_t i, j;
388        std::set<size_t>::const_iterator itr;
389
390        // check VectorSet is Simple Vector class with sets for elements
391        CheckSimpleVector<std::set<size_t>, VectorSet>(
392                one_element_std_set<size_t>(), two_element_std_set<size_t>()
393        );
394
395        // domain dimensions for F
396        size_t n = ind_taddr.size();
397        size_t m = dep_taddr.size();
398
399        CPPAD_ASSERT_KNOWN(
400                q > 0,
401                "RevSparseJac: q is not greater than zero"
402        );
403        CPPAD_ASSERT_KNOWN(
404                size_t(r.size()) == q || transpose,
405                "RevSparseJac: size of r is not equal to q and transpose is false."
406        );
407        CPPAD_ASSERT_KNOWN(
408                size_t(r.size()) == m || ! transpose,
409                "RevSparseJac: size of r is not equal to m and transpose is true."
410        );
411
412        // vector of lists that will hold the results
413        CPPAD_INTERNAL_SPARSE_SET    var_sparsity;
414        var_sparsity.resize(total_num_var, q);
415
416        // The sparsity pattern corresponding to the dependent variables
417        if( transpose )
418        {       for(i = 0; i < m; i++)
419                {       itr = r[i].begin();
420                        while(itr != r[i].end())
421                        {       j = *itr++;
422                                CPPAD_ASSERT_KNOWN(
423                                j < q,
424                                "RevSparseJac: transpose is true and element of the set\n"
425                                "r[i] has value greater than or equal q."
426                                );
427                                CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var );
428                                var_sparsity.add_element( dep_taddr[i], j );
429                        }
430                }
431        }
432        else
433        {       for(i = 0; i < q; i++)
434                {       itr = r[i].begin();
435                        while(itr != r[i].end())
436                        {       j = *itr++;
437                                CPPAD_ASSERT_KNOWN(
438                                j < m,
439                                "RevSparseJac: transpose is false and element of the set\n"
440                                "r[i] has value greater than or equal range dimension."
441                                );
442                                CPPAD_ASSERT_UNKNOWN( dep_taddr[j] < total_num_var );
443                                var_sparsity.add_element( dep_taddr[j], i );
444                        }
445                }
446        }
447        // evaluate the sparsity patterns
448        RevJacSweep(
449                dependency,
450                n,
451                total_num_var,
452                &play,
453                var_sparsity
454        );
455
456        // return values corresponding to dependent variables
457        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == q || transpose );
458        CPPAD_ASSERT_UNKNOWN( size_t(s.size()) == n || ! transpose );
459        for(j = 0; j < n; j++)
460        {       CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) );
461
462                // ind_taddr[j] is operator taddr for j-th independent variable
463                CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp );
464
465                CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
466                var_sparsity.begin(j+1);
467                i = var_sparsity.next_element();
468                while( i < q )
469                {       if( transpose )
470                                s[j].insert(i);
471                        else    s[i].insert(j);
472                        i = var_sparsity.next_element();
473                }
474        }
475}
476// =========================================================================
477// RevSparseJacCase
478
479/*!
480Private helper function for \c RevSparseJac(q, r, transpose).
481
482All of the description in the public member function
483\c RevSparseJac(q, r, transpose) apply.
484
485\param set_type
486is a \c bool value.
487This argument is used to dispatch to the proper source code
488depending on the value of \c VectorSet::value_type.
489
490\param transpose
491See \c RevSparseJac(q, r, transpose, dependency)
492
493\param dependency
494See \c RevSparseJac(q, r, transpose, dependency)
495
496\param q
497See \c RevSparseJac(q, r, transpose, dependency)
498
499\param r
500See \c RevSparseJac(q, r, transpose, dependency)
501
502\param s
503is the return value for the corresponding call to
504RevSparseJac(q, r, transpose).
505*/
506
507template <class Base>
508template <class VectorSet>
509void ADFun<Base>::RevSparseJacCase(
510        bool                set_type          ,
511        bool                transpose         ,
512        bool                dependency        ,
513        size_t              q                 ,
514        const VectorSet&    r                 ,
515        VectorSet&          s                 )
516{       size_t n = Domain();
517
518        // dimension of the result vector
519        s.resize( q * n );
520
521        // store results in s
522        RevSparseJacBool(
523                transpose      ,
524                dependency     ,
525                q              ,
526                r              ,
527                s              ,
528                num_var_tape_  ,
529                dep_taddr_     ,
530                ind_taddr_     ,
531                play_
532        );
533}
534
535/*!
536Private helper function for \c RevSparseJac(q, r, transpose).
537
538All of the description in the public member function
539\c RevSparseJac(q, r, transpose) apply.
540
541\param set_type
542is a \c std::set<size_t> object.
543This argument is used to dispatch to the proper source code
544depending on the value of \c VectorSet::value_type.
545
546\param transpose
547See \c RevSparseJac(q, r, transpose, dependency)
548
549\param dependency
550See \c RevSparseJac(q, r, transpose, dependency)
551
552\param q
553See \c RevSparseJac(q, r, transpose, dependency)
554
555\param r
556See \c RevSparseJac(q, r, transpose, dependency)
557
558\param s
559is the return value for the corresponding call to RevSparseJac(q, r, transpose)
560*/
561
562template <class Base>
563template <class VectorSet>
564void ADFun<Base>::RevSparseJacCase(
565        const std::set<size_t>&      set_type          ,
566        bool                         transpose         ,
567        bool                         dependency        ,
568        size_t                       q                 ,
569        const VectorSet&             r                 ,
570        VectorSet&                   s                 )
571{       // dimension of the result vector
572        if( transpose )
573                s.resize( Domain() );
574        else    s.resize( q );
575
576        // store results in r
577        RevSparseJacSet(
578                transpose      ,
579                dependency     ,
580                q              ,
581                r              ,
582                s              ,
583                num_var_tape_  ,
584                dep_taddr_     ,
585                ind_taddr_     ,
586                play_
587        );
588}
589// =========================================================================
590// RevSparseJac
591/*!
592User API for Jacobian sparsity patterns using reverse mode.
593
594The C++ source code corresponding to this operation is
595\verbatim
596        s = f.RevSparseJac(q, r, transpose, dependency)
597\endverbatim
598
599\tparam Base
600is the base type for this recording.
601
602\tparam VectorSet
603is a simple vector with elements of type \c bool.
604or \c std::set<size_t>.
605
606\param q
607is the number of rows in the matrix \f$ R \f$.
608
609\param r
610is a sparsity pattern for the matrix \f$ R \f$.
611
612\param transpose
613are the sparsity patterns for \f$ R \f$ and \f$ S(x) \f$ transposed.
614
615\param dependency
616Are the derivatives with respect to left and right of the expression below
617considered to be non-zero:
618\code
619        CondExpRel(left, right, if_true, if_false)
620\endcode
621This is used by the optimizer to obtain the correct dependency relations.
622
623
624\return
625If \c transpose is false (true), the return value is a sparsity pattern
626for \f$ S(x) \f$ (\f$ S(x)^T \f$) where
627\f[
628        S(x) = R * F^{(1)} (x)
629\f]
630and \f$ F \f$ is the function corresponding to the operation sequence
631and \a x is any argument value.
632If \c VectorSet::value_type is \c bool,
633the return value has size \f$ q * n \f$ ( \f$ n * q \f$).
634If \c VectorSet::value_type is \c std::set<size_t>,
635the return value has size \f$ q \f$ ( \f$ n \f$)
636and with all its elements between zero and \f$ n - 1 \f$ (\f$ q - 1 \f$).
637*/
638template <class Base>
639template <class VectorSet>
640VectorSet ADFun<Base>::RevSparseJac(
641        size_t              q          ,
642        const VectorSet&    r          ,
643        bool                transpose  ,
644        bool                dependency )
645{
646        VectorSet s;
647        typedef typename VectorSet::value_type Set_type;
648
649        RevSparseJacCase(
650                Set_type()    ,
651                transpose     ,
652                dependency    ,
653                q             ,
654                r             ,
655                s
656        );
657        return s;
658}
659// ===========================================================================
660// RevSparseJacCheckpoint
661/*!
662Reverse mode Jacobian sparsity calculation used by checkpoint functions.
663
664\tparam Base
665is the base type for this recording.
666
667\param transpose
668is true (false) s is equal to \f$ S(x) \f$ (\f$ S(x)^T \f$)
669where
670\f[
671        S(x) = R * F^{(1)} (x)
672\f]
673where \f$ F \f$ is the function corresponding to the operation sequence
674and \f$ x \f$ is any argument value.
675
676\param q
677is the number of rows in the matrix \f$ R \f$.
678
679\param r
680is a sparsity pattern for the matrix \f$ R \f$.
681
682\param transpose
683are the sparsity patterns for \f$ R \f$ and \f$ S(x) \f$ transposed.
684
685\param dependency
686Are the derivatives with respect to left and right of the expression below
687considered to be non-zero:
688\code
689        CondExpRel(left, right, if_true, if_false)
690\endcode
691This is used by the optimizer to obtain the correct dependency relations.
692
693\param s
694The input size and elements of s do not matter.
695On output, s is the sparsity pattern for the matrix \f$ S(x) \f$
696or \f$ S(x)^T \f$ depending on transpose.
697
698*/
699template <class Base>
700void ADFun<Base>::RevSparseJacCheckpoint(
701        size_t                        q          ,
702        CPPAD_INTERNAL_SPARSE_SET&    r          ,
703        bool                          transpose  ,
704        bool                          dependency ,
705        CPPAD_INTERNAL_SPARSE_SET&    s          )
706{       size_t n = Domain();
707        size_t m = Range();
708
709# ifndef NDEBUG
710        if( transpose )
711        {       CPPAD_ASSERT_UNKNOWN( r.n_set() == m );
712                CPPAD_ASSERT_UNKNOWN( r.end()   == q );
713        }
714        else
715        {       CPPAD_ASSERT_UNKNOWN( r.n_set() == q );
716                CPPAD_ASSERT_UNKNOWN( r.end()   == m );
717        }
718        for(size_t i = 0; i < m; i++)
719                CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
720# endif
721
722        // holds reverse Jacobian sparsity pattern for all variables
723        CPPAD_INTERNAL_SPARSE_SET var_sparsity;
724        var_sparsity.resize(num_var_tape_, q);
725
726        // set sparsity pattern for dependent variables
727        if( transpose )
728        {       for(size_t i = 0; i < m; i++)
729                {       r.begin(i);
730                        size_t j = r.next_element();
731                        while( j < q )
732                        {       var_sparsity.add_element( dep_taddr_[i], j );
733                                j = r.next_element();
734                        }
735                }
736        }
737        else
738        {       for(size_t j = 0; j < q; j++)
739                {       r.begin(j);
740                        size_t i = r.next_element();
741                        while( i < m )
742                        {       var_sparsity.add_element( dep_taddr_[i], j );
743                                i = r.next_element();
744                        }
745                }
746        }
747
748        // evaluate the sparsity pattern for all variables
749        RevJacSweep(
750                dependency,
751                n,
752                num_var_tape_,
753                &play_,
754                var_sparsity
755        );
756
757        // dimension the return value
758        if( transpose )
759                s.resize(n, m);
760        else
761                s.resize(m, n);
762
763        // return values corresponding to independent variables
764        for(size_t j = 0; j < n; j++)
765        {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == (j+1) );
766
767                // ind_taddr_[j] is operator taddr for j-th independent variable
768                CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );
769
770                // extract the result from var_sparsity
771                CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
772                var_sparsity.begin(j+1);
773                size_t i = var_sparsity.next_element();
774                while( i < q )
775                {       if( transpose )
776                                s.add_element(j, i);
777                        else
778                                s.add_element(i, j);
779                        i  = var_sparsity.next_element();
780                }
781        }
782
783}
784
785} // END_CPPAD_NAMESPACE
786# endif
Note: See TracBrowser for help on using the repository browser.