Changeset 3108


Ignore:
Timestamp:
Feb 20, 2014 8:51:45 AM (6 years ago)
Author:
bradbell
Message:

sparse_jacobian.hpp: sort by color (simpler and faster).
package.sh: do not include files in the new directories.
sparse_jacobian.cpp: add missing test of sparsity pattern argument.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/bin/package.sh

    r3081 r3108  
    115115        \( -name '*.sh' \) -or \
    116116        \( -name '*.txt' \) | sed \
     117                -e '/\/new\//d' \
    117118                -e '/\.\/build\//d' \
    118119                -e '/bug\/build\//d' \
  • trunk/cppad/local/sparse_jacobian.hpp

    r3107 r3108  
    115115They specify which rows and columns of $latex F^{(1)} (x)$$ are
    116116computes and in what order.
     117Not all the non-zero entries in $latex F^{(1)} (x)$$ need be computed,
     118but all the entries specified by $icode row$$ and $icode col$$
     119must be possibly non-zero in the sparsity pattern.
    117120We use $latex K$$ to denote the value $icode%jac%.size()%$$
    118121which must also equal the size of $icode row$$ and $icode col$$.
     
    120123for $latex k = 0 , \ldots , K-1$$, it must hold that
    121124$latex row[k] < m$$ and $latex col[k] < n$$.
    122 All of the $latex (row[k], col[k])$$ pairs must correspond to an entry
    123 in the sparsity pattern $icode p$$,
    124 but not all such entries need be included
    125 in the $latex (row[k], col[k])$$ pairs.
    126125
    127126$head jac$$
     
    252251class sparse_jacobian_work {
    253252        public:
    254                 /// indices that sort the user arrays by row
    255                 /// with the extra value K at the end
    256                 CppAD::vector<size_t> sort_row;
    257                 /// indices that sort the user arrays by column
    258                 /// with the extra value K at the end
    259                 CppAD::vector<size_t> sort_col;
     253                /// indices that sort the user row and col arrays by color
     254                CppAD::vector<size_t> order;
    260255                /// results of the coloring algorithm
    261256                CppAD::vector<size_t> color;
     
    263258                void clear(void)
    264259                {
    265                         sort_row.clear();
    266                         sort_col.clear();
     260                        order.clear();
    267261                        color.clear();
    268262                }
     
    281275\tparam VectorSet
    282276is either \c sparse_pack, \c sparse_set or \c sparse_list.
     277
     278\tparam VectorSize
     279is a simple vector class with elements of type \c size_t.
    283280
    284281\param x [in]
     
    330327        size_t i, j, k, ell;
    331328
    332         CppAD::vector<size_t>& sort_col(work.sort_col);
     329        CppAD::vector<size_t>& order(work.order);
    333330        CppAD::vector<size_t>& color(work.color);
    334331
     
    359356                CPPAD_ASSERT_UNKNOWN( p_transpose.end() ==  m );
    360357
    361                 // put sorting indices in sort_col
    362                 work.sort_col.resize(K);
    363                 index_sort(col, work.sort_col);
    364 
    365358                // initialize columns that are in the returned jacobian
    366359                CppAD::vector<bool> col_used(n);
     
    372365                c2r_used.resize(n, m);
    373366                r2c_used.resize(m, n);
    374                 k = 0;
    375                 while( k < K )
    376                 {       CPPAD_ASSERT_UNKNOWN(
    377                                 row[sort_col[k]] < m && col[sort_col[k]] < n
    378                         );
    379                         CPPAD_ASSERT_UNKNOWN(
    380                                 k == 0 || col[sort_col[k-1]] <= col[sort_col[k]]
    381                         );
    382                         CPPAD_ASSERT_KNOWN(
    383                                 p_transpose.is_element(col[sort_col[k]], row[sort_col[k]]) ,
    384                                 "SparseJacobianForward: "
    385                                 "an (row, col) pair is not in sparsity pattern."
     367                for(k = 0; k < K; k++)
     368                {       CPPAD_ASSERT_KNOWN(
     369                                p_transpose.is_element(col[k], row[k]) ,
     370                                "SparseJacobianForward: a (row, col) pair "
     371                                "is not in the sparsity pattern."
    386372                        );
    387373                        col_used[ col[k] ] = true;
    388374                        c2r_used.add_element(col[k], row[k]);
    389375                        r2c_used.add_element(row[k], col[k]);
    390                         k++;
    391376                }
    392377       
     
    468453                        color[j] = ell;
    469454                }
     455
     456                // put sorting indices in color order
     457                VectorSize key(K);
     458                order.resize(K);
     459                for(k = 0; k < K; k++)
     460                        key[k] = color[ col[k] ];
     461                index_sort(key, order);
    470462        }
    471463        size_t n_color = 1;
     
    484476
    485477        // loop over colors
    486         size_t n_sweep = 0;
     478        k = 0;
    487479        for(ell = 0; ell < n_color; ell++)
    488         {       n_sweep++;
     480        {       CPPAD_ASSERT_UNKNOWN( color[ col[ order[k] ] ] == ell );
    489481
    490482                // combine all columns with this color
     
    498490
    499491                // set the corresponding components of the result
    500                 k = 0;
    501                 for(j = 0; j < n; j++) if( color[j] == ell )
    502                 {       // find first k index for this column
    503                         while( k < K && col[sort_col[k]] < j )
    504                                 k++;
    505                         // extract the row results for this column
    506                         while( k < K && col[sort_col[k]] == j )
    507                         {       jac[ sort_col[k] ] = dy[row[sort_col[k]]];
    508                                 k++;
    509                         }
     492                while( k < K && color[ col[order[k]] ] == ell )
     493                {       jac[ order[k] ] = dy[row[order[k]]];
     494                        k++;
    510495                }
    511496        }
    512         return n_sweep;
     497        return n_color;
    513498}
    514499/*!
     
    524509\tparam VectorSet
    525510is either \c sparse_pack, \c sparse_set or \c sparse_list.
     511
     512\tparam VectorSize
     513is a simple vector class with elements of type \c size_t.
    526514
    527515\param x [in]
     
    573561        size_t i, j, k, ell;
    574562
    575         CppAD::vector<size_t>& sort_row(work.sort_row);
     563        CppAD::vector<size_t>& order(work.order);
    576564        CppAD::vector<size_t>& color(work.color);
    577565
     
    601589                CPPAD_ASSERT_UNKNOWN( p.end() ==  n );
    602590
    603                 // put sorting indices in sort_row
    604                 work.sort_row.resize(K);
    605                 index_sort(row, work.sort_row);
    606 
    607591                // initialize rows that are in the returned jacobian
    608592                CppAD::vector<bool> row_used(m);
     
    614598                c2r_used.resize(n, m);
    615599                r2c_used.resize(m, n);
    616                 k = 0;
    617                 while( k < K )
    618                 {       CPPAD_ASSERT_UNKNOWN(
    619                                 row[sort_row[k]] < m && col[sort_row[k]] < n
    620                         );
    621                         CPPAD_ASSERT_UNKNOWN(
    622                                 k == 0 || row[sort_row[k-1]] <= row[sort_row[k]]
    623                         );
    624                         CPPAD_ASSERT_KNOWN(
    625                                 p.is_element(row[sort_row[k]], col[sort_row[k]]) ,
    626                                 "SparseJacobianReverse: "
    627                                 "an (row, col) pair is not in sparsity pattern."
     600                for(k = 0;  k < K; k++)
     601                {       CPPAD_ASSERT_KNOWN(
     602                                p.is_element(row[k], col[k]) ,
     603                                "SparseJacobianReverse: a (row, col) pair "
     604                                "is not in the sparsity pattern."
    628605                        );
    629606                        row_used[ row[k] ] = true;
    630607                        c2r_used.add_element(col[k], row[k]);
    631608                        r2c_used.add_element(row[k], col[k]);
    632                         k++;
    633609                }
    634610       
     
    713689                        color[i] = ell;
    714690                }
     691
     692                // put sorting indices in color order
     693                VectorSize key(K);
     694                order.resize(K);
     695                for(k = 0; k < K; k++)
     696                        key[k] = color[ row[k] ];
     697                index_sort(key, order);
    715698        }
    716699        size_t n_color = 1;
     
    729712
    730713        // loop over colors
    731         size_t n_sweep = 0;
     714        k = 0;
    732715        for(ell = 0; ell < n_color; ell++)
    733         {       n_sweep++;
     716        {       CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell );
    734717
    735718                // combine all the rows with this color
     
    743726
    744727                // set the corresponding components of the result
    745                 k = 0;
    746                 for(i = 0; i < m; i++) if( color[i] == ell )
    747                 {       // find first k index for this row
    748                         while( k < K && row[sort_row[k]] < i )
    749                                 k++;
    750                         // extract the row results for this row
    751                         while( k < K && row[sort_row[k]] == i )
    752                         {       jac[ sort_row[k] ] = dw[col[sort_row[k]]];
    753                                 k++;
    754                         }
     728                while( k < K && color[ row[order[k]] ]  == ell )
     729                {       jac[ order[k] ] = dw[col[order[k]]];
     730                        k++;
    755731                }
    756732        }
    757         return n_sweep;
     733        return n_color;
    758734}
    759735// ==========================================================================
  • trunk/test_more/sparse_jacobian.cpp

    r2870 r3108  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    319319        }
    320320        p   = f.RevSparseJac(m, s);
    321         jac = f.SparseJacobian(x);
     321        jac = f.SparseJacobian(x, p);
    322322        for(k = 0; k < 12; k++)
    323323                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     
    376376                s[i].insert(i);
    377377        p   = f.RevSparseJac(m, s);
    378         jac = f.SparseJacobian(x);
     378        jac = f.SparseJacobian(x, p);
    379379        for(k = 0; k < 12; k++)
    380380                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     
    440440        }
    441441        p   = f.ForSparseJac(n, r);
    442         jac = f.SparseJacobian(x);
     442        jac = f.SparseJacobian(x, p);
    443443        for(k = 0; k < 12; k++)
    444444                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
     
    500500                r[j].insert(j);
    501501        p   = f.ForSparseJac(n, r);
    502         jac = f.SparseJacobian(x);
     502        jac = f.SparseJacobian(x, p);
    503503        for(k = 0; k < 12; k++)
    504504                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
Note: See TracChangeset for help on using the changeset viewer.