1 # ifndef CPPAD_EXAMPLE_ATOMIC_THREE_MAT_MUL_HPP 2 # define CPPAD_EXAMPLE_ATOMIC_THREE_MAT_MUL_HPP 3 /* -------------------------------------------------------------------------- 4 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-19 Bradley M. Bell 5 6 CppAD is distributed under the terms of the 7 Eclipse Public License Version 2.0. 8 9 This Source Code may also be made available under the following 10 Secondary License when the conditions for such availability set forth 11 in the Eclipse Public License, Version 2.0 are satisfied: 12 GNU General Public License, Version 2.0 or later. 13 ---------------------------------------------------------------------------- */ 14 15 /* 16 $begin atomic_three_mat_mul.hpp$$ 17 $spell 18 Taylor 19 ty 20 px 21 CppAD 22 jac 23 hes 24 nr 25 nc 26 afun 27 mul 28 $$ 29 30 $section Matrix Multiply as an Atomic Operation$$ 31 32 $head See Also$$ 33 $cref atomic_two_eigen_mat_mul.hpp$$ 34 35 $head Purpose$$ 36 Use scalar $code double$$ operations in an $cref atomic_three$$ operation 37 that computes the matrix product for $code AD<double$$ operations. 38 39 $subhead parameter_x$$ 40 This example demonstrates the use of the 41 $cref/parameter_x/atomic_three/parameter_x/$$ 42 argument to the $cref atomic_three$$ virtual functions. 43 44 $subhead type_x$$ 45 This example also demonstrates the use of the 46 $cref/type_x/atomic_three/type_x/$$ 47 argument to the $cref atomic_three$$ virtual functions. 48 49 $head Matrix Dimensions$$ 50 The first three components of the argument vector $icode ax$$ 51 in the call $icode%afun%(%ax%, %ay%)%$$ 52 are parameters and contain the matrix dimensions. 53 This enables them to be different for each use of the same atomic 54 function $icode afun$$. 55 These dimensions are: 56 $table 57 $icode%ax%[0]%$$ $pre $$ 58 $cnext $icode nr_left$$ $pre $$ 59 $cnext number of rows in the left matrix and result matrix 60 $rend 61 $icode%ax%[1]%$$ $pre $$ 62 $cnext $icode n_middle$$ $pre $$ 63 $cnext columns in the left matrix and rows in right matrix 64 $rend 65 $icode%ax%[2]%$$ $pre $$ 66 $cnext $icode nc_right$$ $pre $$ 67 $cnext number of columns in the right matrix and result matrix 68 $tend 69 70 71 $head Left Matrix$$ 72 The number of elements in the left matrix is 73 $codei% 74 %n_left% = %nr_left% * %n_middle% 75 %$$ 76 The elements are in 77 $icode%ax%[3]%$$ through $icode%ax%[2+%n_left%]%$$ in row major order. 78 79 $head Right Matrix$$ 80 The number of elements in the right matrix is 81 $codei% 82 %n_right% = %n_middle% * %nc_right% 83 %$$ 84 The elements are in 85 $icode%ax%[3+%n_left%]%$$ through 86 $icode%ax%[2+%n_left%+%n_right%]%$$ in row major order. 87 88 $head Result Matrix$$ 89 The number of elements in the result matrix is 90 $codei% 91 %n_result% = %nr_left% * %nc_right% 92 %$$ 93 The elements are in 94 $icode%ay%[0]%$$ through $icode%ay%[%n_result%-1]%$$ in row major order. 95 96 $head Start Class Definition$$ 97 $srccode%cpp% */ 98 # include <cppad/cppad.hpp> 99 namespace { // Begin empty namespace 100 using CppAD::vector; 101 // 102 // matrix result = left * right 103 class atomic_mat_mul : public CppAD::atomic_three<double> { 104 /* %$$ 105 $head Constructor$$ 106 $srccode%cpp% */ 107 public: 108 // --------------------------------------------------------------------- 109 // constructor atomic_mat_mul(void)110 atomic_mat_mul(void) : CppAD::atomic_three<double>("mat_mul") 111 { } 112 private: 113 /* %$$ 114 $head Left Operand Element Index$$ 115 Index in the Taylor coefficient matrix $icode tx$$ of a left matrix element. 116 $srccode%cpp% */ left(size_t i,size_t j,size_t k,size_t nk,size_t nr_left,size_t n_middle,size_t nc_right)117 size_t left( 118 size_t i , // left matrix row index 119 size_t j , // left matrix column index 120 size_t k , // Taylor coeffocient order 121 size_t nk , // number of Taylor coefficients in tx 122 size_t nr_left , // rows in left matrix 123 size_t n_middle , // rows in left and columns in right 124 size_t nc_right ) // columns in right matrix 125 { assert( i < nr_left ); 126 assert( j < n_middle ); 127 return (3 + i * n_middle + j) * nk + k; 128 } 129 /* %$$ 130 $head Right Operand Element Index$$ 131 Index in the Taylor coefficient matrix $icode tx$$ of a right matrix element. 132 $srccode%cpp% */ right(size_t i,size_t j,size_t k,size_t nk,size_t nr_left,size_t n_middle,size_t nc_right)133 size_t right( 134 size_t i , // right matrix row index 135 size_t j , // right matrix column index 136 size_t k , // Taylor coeffocient order 137 size_t nk , // number of Taylor coefficients in tx 138 size_t nr_left , // rows in left matrix 139 size_t n_middle , // rows in left and columns in right 140 size_t nc_right ) // columns in right matrix 141 { assert( i < n_middle ); 142 assert( j < nc_right ); 143 size_t offset = 3 + nr_left * n_middle; 144 return (offset + i * nc_right + j) * nk + k; 145 } 146 /* %$$ 147 $head Result Element Index$$ 148 Index in the Taylor coefficient matrix $icode ty$$ of a result matrix element. 149 $srccode%cpp% */ result(size_t i,size_t j,size_t k,size_t nk,size_t nr_left,size_t n_middle,size_t nc_right)150 size_t result( 151 size_t i , // result matrix row index 152 size_t j , // result matrix column index 153 size_t k , // Taylor coeffocient order 154 size_t nk , // number of Taylor coefficients in ty 155 size_t nr_left , // rows in left matrix 156 size_t n_middle , // rows in left and columns in right 157 size_t nc_right ) // columns in right matrix 158 { assert( i < nr_left ); 159 assert( j < nc_right ); 160 return (i * nc_right + j) * nk + k; 161 } 162 /* %$$ 163 $head Forward Matrix Multiply$$ 164 Forward mode multiply Taylor coefficients in $icode tx$$ and sum into 165 $icode ty$$ (for one pair of left and right orders) 166 $srccode%cpp% */ forward_multiply(size_t k_left,size_t k_right,const vector<double> & tx,vector<double> & ty,size_t nr_left,size_t n_middle,size_t nc_right)167 void forward_multiply( 168 size_t k_left , // order for left coefficients 169 size_t k_right , // order for right coefficients 170 const vector<double>& tx , // domain space Taylor coefficients 171 vector<double>& ty , // range space Taylor coefficients 172 size_t nr_left , // rows in left matrix 173 size_t n_middle , // rows in left and columns in right 174 size_t nc_right ) // columns in right matrix 175 { 176 size_t nx = 3 + (nr_left + nc_right) * n_middle; 177 size_t nk = tx.size() / nx; 178 # ifndef NDEBUG 179 size_t ny = nr_left * nc_right; 180 assert( nk == ty.size() / ny ); 181 # endif 182 // 183 size_t k_result = k_left + k_right; 184 assert( k_result < nk ); 185 // 186 for(size_t i = 0; i < nr_left; i++) 187 { for(size_t j = 0; j < nc_right; j++) 188 { double sum = 0.0; 189 for(size_t ell = 0; ell < n_middle; ell++) 190 { size_t i_left = left( 191 i, ell, k_left, nk, nr_left, n_middle, nc_right 192 ); 193 size_t i_right = right( 194 ell, j, k_right, nk, nr_left, n_middle, nc_right 195 ); 196 sum += tx[i_left] * tx[i_right]; 197 } 198 size_t i_result = result( 199 i, j, k_result, nk, nr_left, n_middle, nc_right 200 ); 201 ty[i_result] += sum; 202 } 203 } 204 } 205 /* %$$ 206 $head Reverse Matrix Multiply$$ 207 Reverse mode partials of Taylor coefficients and sum into $icode px$$ 208 (for one pair of left and right orders) 209 $srccode%cpp% */ reverse_multiply(size_t k_left,size_t k_right,const vector<double> & tx,const vector<double> & ty,vector<double> & px,const vector<double> & py,size_t nr_left,size_t n_middle,size_t nc_right)210 void reverse_multiply( 211 size_t k_left , // order for left coefficients 212 size_t k_right , // order for right coefficients 213 const vector<double>& tx , // domain space Taylor coefficients 214 const vector<double>& ty , // range space Taylor coefficients 215 vector<double>& px , // partials w.r.t. tx 216 const vector<double>& py , // partials w.r.t. ty 217 size_t nr_left , // rows in left matrix 218 size_t n_middle , // rows in left and columns in right 219 size_t nc_right ) // columns in right matrix 220 { 221 size_t nx = 3 + (nr_left + nc_right) * n_middle; 222 size_t nk = tx.size() / nx; 223 # ifndef NDEBUG 224 size_t ny = nr_left * nc_right; 225 assert( nk == ty.size() / ny ); 226 # endif 227 assert( tx.size() == px.size() ); 228 assert( ty.size() == py.size() ); 229 // 230 size_t k_result = k_left + k_right; 231 assert( k_result < nk ); 232 // 233 for(size_t i = 0; i < nr_left; i++) 234 { for(size_t j = 0; j < nc_right; j++) 235 { size_t i_result = result( 236 i, j, k_result, nk, nr_left, n_middle, nc_right 237 ); 238 for(size_t ell = 0; ell < n_middle; ell++) 239 { size_t i_left = left( 240 i, ell, k_left, nk, nr_left, n_middle, nc_right 241 ); 242 size_t i_right = right( 243 ell, j, k_right, nk, nr_left, n_middle, nc_right 244 ); 245 // sum += tx[i_left] * tx[i_right]; 246 px[i_left] += tx[i_right] * py[i_result]; 247 px[i_right] += tx[i_left] * py[i_result]; 248 } 249 } 250 } 251 return; 252 } 253 /* %$$ 254 $head for_type$$ 255 Routine called by CppAD during $cref/afun(ax, ay)/atomic_three_afun/$$. 256 $srccode%cpp% */ 257 // calculate type_y for_type(const vector<double> & parameter_x,const vector<CppAD::ad_type_enum> & type_x,vector<CppAD::ad_type_enum> & type_y)258 virtual bool for_type( 259 const vector<double>& parameter_x , 260 const vector<CppAD::ad_type_enum>& type_x , 261 vector<CppAD::ad_type_enum>& type_y ) 262 { assert( parameter_x.size() == type_x.size() ); 263 bool ok = true; 264 ok &= type_x[0] == CppAD::constant_enum; 265 ok &= type_x[1] == CppAD::constant_enum; 266 ok &= type_x[2] == CppAD::constant_enum; 267 if( ! ok ) 268 return false; 269 // 270 size_t nr_left = size_t( parameter_x[0] ); 271 size_t n_middle = size_t( parameter_x[1] ); 272 size_t nc_right = size_t( parameter_x[2] ); 273 // 274 ok &= type_x.size() == 3 + (nr_left + nc_right) * n_middle; 275 ok &= type_y.size() == n_middle * nc_right; 276 if( ! ok ) 277 return false; 278 // 279 // commpute type_y 280 size_t nk = 1; // number of orders 281 size_t k = 0; // order 282 for(size_t i = 0; i < nr_left; ++i) 283 { for(size_t j = 0; j < nc_right; ++j) 284 { // compute type for result[i, j] 285 CppAD::ad_type_enum type_yij = CppAD::constant_enum; 286 for(size_t ell = 0; ell < n_middle; ++ell) 287 { // index for left(i, ell) 288 size_t i_left = left( 289 i, ell, k, nk, nr_left, n_middle, nc_right 290 ); 291 // indx for right(ell, j) 292 size_t i_right = right( 293 ell, j, k, nk, nr_left, n_middle, nc_right 294 ); 295 // multiplication on left or right by the constant zero 296 // always results in a constant 297 bool zero_left = type_x[i_left] == CppAD::constant_enum; 298 zero_left &= parameter_x[i_left] == 0.0; 299 bool zero_right = type_x[i_right] == CppAD::constant_enum; 300 zero_right &= parameter_x[i_right] == 0.0; 301 if( ! (zero_left | zero_right) ) 302 { type_yij = std::max(type_yij, type_x[i_left] ); 303 type_yij = std::max(type_yij, type_x[i_right] ); 304 } 305 } 306 size_t i_result = result( 307 i, j, k, nk, nr_left, n_middle, nc_right 308 ); 309 type_y[i_result] = type_yij; 310 } 311 } 312 return true; 313 } 314 /* %$$ 315 $head forward$$ 316 Routine called by CppAD during $cref Forward$$ mode. 317 $srccode%cpp% */ forward(const vector<double> & parameter_x,const vector<CppAD::ad_type_enum> & type_x,size_t need_y,size_t q,size_t p,const vector<double> & tx,vector<double> & ty)318 virtual bool forward( 319 const vector<double>& parameter_x , 320 const vector<CppAD::ad_type_enum>& type_x , 321 size_t need_y , 322 size_t q , 323 size_t p , 324 const vector<double>& tx , 325 vector<double>& ty ) 326 { size_t n_order = p + 1; 327 size_t nr_left = size_t( tx[ 0 * n_order + 0 ] ); 328 size_t n_middle = size_t( tx[ 1 * n_order + 0 ] ); 329 size_t nc_right = size_t( tx[ 2 * n_order + 0 ] ); 330 # ifndef NDEBUG 331 size_t nx = 3 + (nr_left + nc_right) * n_middle; 332 size_t ny = nr_left * nc_right; 333 # endif 334 assert( nx * n_order == tx.size() ); 335 assert( ny * n_order == ty.size() ); 336 size_t i, j, ell; 337 338 // initialize result as zero 339 size_t k; 340 for(i = 0; i < nr_left; i++) 341 { for(j = 0; j < nc_right; j++) 342 { for(k = q; k <= p; k++) 343 { size_t i_result = result( 344 i, j, k, n_order, nr_left, n_middle, nc_right 345 ); 346 ty[i_result] = 0.0; 347 } 348 } 349 } 350 for(k = q; k <= p; k++) 351 { // sum the produces that result in order k 352 for(ell = 0; ell <= k; ell++) 353 forward_multiply( 354 ell, k - ell, tx, ty, nr_left, n_middle, nc_right 355 ); 356 } 357 358 // all orders are implemented, so always return true 359 return true; 360 } 361 /* %$$ 362 $head reverse$$ 363 Routine called by CppAD during $cref Reverse$$ mode. 364 $srccode%cpp% */ reverse(const vector<double> & parameter_x,const vector<CppAD::ad_type_enum> & type_x,size_t p,const vector<double> & tx,const vector<double> & ty,vector<double> & px,const vector<double> & py)365 virtual bool reverse( 366 const vector<double>& parameter_x , 367 const vector<CppAD::ad_type_enum>& type_x , 368 size_t p , 369 const vector<double>& tx , 370 const vector<double>& ty , 371 vector<double>& px , 372 const vector<double>& py ) 373 { size_t n_order = p + 1; 374 size_t nr_left = size_t( tx[ 0 * n_order + 0 ] ); 375 size_t n_middle = size_t( tx[ 1 * n_order + 0 ] ); 376 size_t nc_right = size_t( tx[ 2 * n_order + 0 ] ); 377 # ifndef NDEBUG 378 size_t nx = 3 + (nr_left + nc_right) * n_middle; 379 size_t ny = nr_left * nc_right; 380 # endif 381 assert( nx * n_order == tx.size() ); 382 assert( ny * n_order == ty.size() ); 383 assert( px.size() == tx.size() ); 384 assert( py.size() == ty.size() ); 385 386 // initialize summation 387 for(size_t i = 0; i < px.size(); i++) 388 px[i] = 0.0; 389 390 // number of orders to differentiate 391 size_t k = n_order; 392 while(k--) 393 { // differentiate the produces that result in order k 394 for(size_t ell = 0; ell <= k; ell++) 395 reverse_multiply( 396 ell, k - ell, tx, ty, px, py, nr_left, n_middle, nc_right 397 ); 398 } 399 400 // all orders are implented, so always return true 401 return true; 402 } 403 /* %$$ 404 $head jac_sparsity$$ 405 $srccode%cpp% */ 406 // Jacobian sparsity routine called by CppAD jac_sparsity(const vector<double> & parameter_x,const vector<CppAD::ad_type_enum> & type_x,bool dependency,const vector<bool> & select_x,const vector<bool> & select_y,CppAD::sparse_rc<vector<size_t>> & pattern_out)407 virtual bool jac_sparsity( 408 const vector<double>& parameter_x , 409 const vector<CppAD::ad_type_enum>& type_x , 410 bool dependency , 411 const vector<bool>& select_x , 412 const vector<bool>& select_y , 413 CppAD::sparse_rc< vector<size_t> >& pattern_out ) 414 { 415 size_t n = select_x.size(); 416 size_t m = select_y.size(); 417 assert( parameter_x.size() == n ); 418 assert( type_x.size() == n ); 419 // 420 size_t nr_left = size_t( parameter_x[0] ); 421 size_t n_middle = size_t( parameter_x[1] ); 422 size_t nc_right = size_t( parameter_x[2] ); 423 size_t nk = 1; // only one order 424 size_t k = 0; // order zero 425 // 426 // count number of non-zeros in sparsity pattern 427 size_t nnz = 0; 428 for(size_t i = 0; i < nr_left; ++i) 429 { for(size_t j = 0; j < nc_right; ++j) 430 { size_t i_result = result( 431 i, j, k, nk, nr_left, n_middle, nc_right 432 ); 433 if( select_y[i_result] ) 434 { for(size_t ell = 0; ell < n_middle; ++ell) 435 { size_t i_left = left( 436 i, ell, k, nk, nr_left, n_middle, nc_right 437 ); 438 size_t i_right = right( 439 ell, j, k, nk, nr_left, n_middle, nc_right 440 ); 441 bool zero_left = 442 type_x[i_left] == CppAD::constant_enum; 443 zero_left &= parameter_x[i_left] == 0.0; 444 bool zero_right = 445 type_x[i_right] == CppAD::constant_enum; 446 zero_right &= parameter_x[i_right] == 0.0; 447 if( ! (zero_left | zero_right ) ) 448 { bool var_left = 449 type_x[i_left] == CppAD::variable_enum; 450 bool var_right = 451 type_x[i_right] == CppAD::variable_enum; 452 if( select_x[i_left] & var_left ) 453 ++nnz; 454 if( select_x[i_right] & var_right ) 455 ++nnz; 456 } 457 } 458 } 459 } 460 } 461 // 462 // fill in the sparsity pattern 463 pattern_out.resize(m, n, nnz); 464 size_t idx = 0; 465 for(size_t i = 0; i < nr_left; ++i) 466 { for(size_t j = 0; j < nc_right; ++j) 467 { size_t i_result = result( 468 i, j, k, nk, nr_left, n_middle, nc_right 469 ); 470 if( select_y[i_result] ) 471 { for(size_t ell = 0; ell < n_middle; ++ell) 472 { size_t i_left = left( 473 i, ell, k, nk, nr_left, n_middle, nc_right 474 ); 475 size_t i_right = right( 476 ell, j, k, nk, nr_left, n_middle, nc_right 477 ); 478 bool zero_left = 479 type_x[i_left] == CppAD::constant_enum; 480 zero_left &= parameter_x[i_left] == 0.0; 481 bool zero_right = 482 type_x[i_right] == CppAD::constant_enum; 483 zero_right &= parameter_x[i_right] == 0.0; 484 if( ! (zero_left | zero_right ) ) 485 { bool var_left = 486 type_x[i_left] == CppAD::variable_enum; 487 bool var_right = 488 type_x[i_right] == CppAD::variable_enum; 489 if( select_x[i_left] & var_left ) 490 pattern_out.set(idx++, i_result, i_left); 491 if( select_x[i_right] & var_right ) 492 pattern_out.set(idx++, i_result, i_right); 493 } 494 } 495 } 496 } 497 } 498 assert( idx == nnz ); 499 // 500 return true; 501 } 502 /* %$$ 503 $head hes_sparsity$$ 504 $srccode%cpp% */ 505 // Jacobian sparsity routine called by CppAD hes_sparsity(const vector<double> & parameter_x,const vector<CppAD::ad_type_enum> & type_x,const vector<bool> & select_x,const vector<bool> & select_y,CppAD::sparse_rc<vector<size_t>> & pattern_out)506 virtual bool hes_sparsity( 507 const vector<double>& parameter_x , 508 const vector<CppAD::ad_type_enum>& type_x , 509 const vector<bool>& select_x , 510 const vector<bool>& select_y , 511 CppAD::sparse_rc< vector<size_t> >& pattern_out ) 512 { 513 size_t n = select_x.size(); 514 assert( parameter_x.size() == n ); 515 assert( type_x.size() == n ); 516 // 517 size_t nr_left = size_t( parameter_x[0] ); 518 size_t n_middle = size_t( parameter_x[1] ); 519 size_t nc_right = size_t( parameter_x[2] ); 520 size_t nk = 1; // only one order 521 size_t k = 0; // order zero 522 // 523 // count number of non-zeros in sparsity pattern 524 size_t nnz = 0; 525 for(size_t i = 0; i < nr_left; ++i) 526 { for(size_t j = 0; j < nc_right; ++j) 527 { size_t i_result = result( 528 i, j, k, nk, nr_left, n_middle, nc_right 529 ); 530 if( select_y[i_result] ) 531 { for(size_t ell = 0; ell < n_middle; ++ell) 532 { // i_left depends on i, ell 533 size_t i_left = left( 534 i, ell, k, nk, nr_left, n_middle, nc_right 535 ); 536 // i_right depens on ell, j 537 size_t i_right = right( 538 ell, j, k, nk, nr_left, n_middle, nc_right 539 ); 540 bool var_left = select_x[i_left] & 541 (type_x[i_left] == CppAD::variable_enum); 542 bool var_right = select_x[i_right] & 543 (type_x[i_right] == CppAD::variable_enum); 544 if( var_left & var_right ) 545 nnz += 2; 546 } 547 } 548 } 549 } 550 // 551 // fill in the sparsity pattern 552 pattern_out.resize(n, n, nnz); 553 size_t idx = 0; 554 for(size_t i = 0; i < nr_left; ++i) 555 { for(size_t j = 0; j < nc_right; ++j) 556 { size_t i_result = result( 557 i, j, k, nk, nr_left, n_middle, nc_right 558 ); 559 if( select_y[i_result] ) 560 { for(size_t ell = 0; ell < n_middle; ++ell) 561 { size_t i_left = left( 562 i, ell, k, nk, nr_left, n_middle, nc_right 563 ); 564 size_t i_right = right( 565 ell, j, k, nk, nr_left, n_middle, nc_right 566 ); 567 bool var_left = select_x[i_left] & 568 (type_x[i_left] == CppAD::variable_enum); 569 bool var_right = select_x[i_right] & 570 (type_x[i_right] == CppAD::variable_enum); 571 if( var_left & var_right ) 572 { // Cannot possibly set the same (i_left, i_right) 573 // pair twice. 574 assert( i_left != i_right ); 575 pattern_out.set(idx++, i_left, i_right); 576 pattern_out.set(idx++, i_right, i_left); 577 } 578 } 579 } 580 } 581 } 582 assert( idx == nnz ); 583 // 584 return true; 585 } 586 /* %$$ 587 $head rev_depend$$ 588 Routine called when a function using $code mat_mul$$ is optimized. 589 $srccode%cpp% */ 590 // calculate depend_x rev_depend(const vector<double> & parameter_x,const vector<CppAD::ad_type_enum> & type_x,vector<bool> & depend_x,const vector<bool> & depend_y)591 virtual bool rev_depend( 592 const vector<double>& parameter_x , 593 const vector<CppAD::ad_type_enum>& type_x , 594 vector<bool>& depend_x , 595 const vector<bool>& depend_y ) 596 { assert( parameter_x.size() == depend_x.size() ); 597 assert( parameter_x.size() == type_x.size() ); 598 bool ok = true; 599 // 600 size_t nr_left = size_t( parameter_x[0] ); 601 size_t n_middle = size_t( parameter_x[1] ); 602 size_t nc_right = size_t( parameter_x[2] ); 603 // 604 ok &= depend_x.size() == 3 + (nr_left + nc_right) * n_middle; 605 ok &= depend_y.size() == n_middle * nc_right; 606 if( ! ok ) 607 return false; 608 // 609 // initialize depend_x 610 for(size_t ell = 0; ell < 3; ++ell) 611 depend_x[ell] = true; // always need these parameters 612 for(size_t ell = 3; ell < depend_x.size(); ++ell) 613 depend_x[ell] = false; // initialize as false 614 // 615 // commpute depend_x 616 size_t nk = 1; // number of orders 617 size_t k = 0; // order 618 for(size_t i = 0; i < nr_left; ++i) 619 { for(size_t j = 0; j < nc_right; ++j) 620 { // check depend for result[i, j] 621 size_t i_result = result( 622 i, j, k, nk, nr_left, n_middle, nc_right 623 ); 624 if( depend_y[i_result] ) 625 { for(size_t ell = 0; ell < n_middle; ++ell) 626 { // index for left(i, ell) 627 size_t i_left = left( 628 i, ell, k, nk, nr_left, n_middle, nc_right 629 ); 630 // indx for right(ell, j) 631 size_t i_right = right( 632 ell, j, k, nk, nr_left, n_middle, nc_right 633 ); 634 bool zero_left = 635 type_x[i_left] == CppAD::constant_enum; 636 zero_left &= parameter_x[i_left] == 0.0; 637 bool zero_right = 638 type_x[i_right] == CppAD::constant_enum; 639 zero_right &= parameter_x[i_right] == 0.0; 640 if( ! zero_right ) 641 depend_x[i_left] = true; 642 if( ! zero_left ) 643 depend_x[i_right] = true; 644 } 645 } 646 } 647 } 648 return true; 649 } 650 /* %$$ 651 $head End Class Definition$$ 652 $srccode%cpp% */ 653 }; // End of mat_mul class 654 } // End empty namespace 655 /* %$$ 656 $comment end nospell$$ 657 $end 658 */ 659 660 661 # endif 662