335 const Teuchos::ParameterList& params)
339 using Teuchos::ParameterList;
340 using Teuchos::Array;
343 typedef typename Cijk_type::i_iterator i_iterator;
344 typedef typename Cijk_type::ik_iterator ik_iterator;
345 typedef typename Cijk_type::ikj_iterator ikj_iterator;
347 const size_type i_tile_size = params.get<OrdinalType>(
"Tile Size");
352 i_iterator i_end = Cijk.
i_end();
353 for (i_iterator i_it=
i_begin; i_it!=i_end; ++i_it) {
354 OrdinalType i = index(i_it);
356 ik_iterator k_end = Cijk.
k_end(i_it);
357 for (ik_iterator k_it =
k_begin; k_it != k_end; ++k_it) {
358 OrdinalType k = index(k_it);
360 ikj_iterator j_end = Cijk.
j_end(k_it);
361 for (ikj_iterator j_it =
j_begin; j_it != j_end; ++j_it) {
362 OrdinalType
j = index(j_it);
364 ValueType c = Stokhos::value(j_it);
365 Cijk_sym.add_term(i,
j, k, c);
370 Cijk_sym.fillComplete();
375 size_type num_i_parts = (basis_size + i_tile_size-1) / i_tile_size;
378 Array<ITile> i_tiles(num_i_parts);
379 for (
size_type i=0; i<num_i_parts; ++i) {
380 i_tiles[i].lower = i*its;
381 i_tiles[i].upper = std::min(basis_size, (i+1)*its);
382 i_tiles[i].parts.resize(1);
383 i_tiles[i].parts[0].lower = basis_size;
384 i_tiles[i].parts[0].upper = 0;
389 for (i_iterator i_it=Cijk_sym.i_begin(); i_it!=Cijk_sym.i_end(); ++i_it) {
390 OrdinalType i = index(i_it);
394 while (idx < num_i_parts && i >= i_tiles[idx].lower) ++idx;
396 TEUCHOS_ASSERT(idx >= 0 && idx < num_i_parts);
398 ik_iterator
k_begin = Cijk_sym.k_begin(i_it);
399 ik_iterator k_end = Cijk_sym.k_end(i_it);
400 for (ik_iterator k_it =
k_begin; k_it != k_end; ++k_it) {
401 OrdinalType
j = index(k_it);
403 if (
j < i_tiles[idx].parts[0].lower)
404 i_tiles[idx].parts[0].lower =
j;
405 if (
j > i_tiles[idx].parts[0].upper)
406 i_tiles[idx].parts[0].upper =
j;
409 for (
size_type idx=0; idx<num_i_parts; ++idx) {
410 size_type lower = i_tiles[idx].parts[0].lower;
411 size_type upper = i_tiles[idx].parts[0].upper;
413 size_type num_j_parts = (range + j_tile_size-1) / j_tile_size;
417 Array<JTile> j_tiles(num_j_parts);
419 j_tiles[
j].lower = lower +
j*jts;
420 j_tiles[
j].upper = std::min(upper+1, lower + (
j+1)*jts);
421 j_tiles[
j].parts.resize(1);
422 j_tiles[
j].parts[0].lower = basis_size;
423 j_tiles[
j].parts[0].upper = 0;
425 i_tiles[idx].parts.swap(j_tiles);
429 for (i_iterator i_it=Cijk_sym.i_begin(); i_it!=Cijk_sym.i_end(); ++i_it) {
430 OrdinalType i = index(i_it);
434 while (idx < num_i_parts && i >= i_tiles[idx].lower) ++idx;
436 TEUCHOS_ASSERT(idx >= 0 && idx < num_i_parts);
438 ik_iterator
k_begin = Cijk_sym.k_begin(i_it);
439 ik_iterator k_end = Cijk_sym.k_end(i_it);
440 for (ik_iterator k_it =
k_begin; k_it != k_end; ++k_it) {
441 OrdinalType
j = index(k_it);
444 size_type num_j_parts = i_tiles[idx].parts.size();
446 while (jdx < num_j_parts && j >= i_tiles[idx].parts[jdx].lower) ++jdx;
448 TEUCHOS_ASSERT(jdx >= 0 && jdx < num_j_parts);
450 ikj_iterator
j_begin = Cijk_sym.j_begin(k_it);
451 ikj_iterator j_end = Cijk_sym.j_end(k_it);
452 for (ikj_iterator j_it =
j_begin; j_it != j_end; ++j_it) {
453 OrdinalType k = index(j_it);
454 ValueType cijk = Stokhos::value(j_it);
458 i_tiles[idx].parts[jdx].parts[0].parts.push_back(
coord);
459 if (k < i_tiles[idx].parts[jdx].parts[0].lower)
460 i_tiles[idx].parts[jdx].parts[0].lower = k;
461 if (k > i_tiles[idx].parts[jdx].parts[0].upper)
462 i_tiles[idx].parts[jdx].parts[0].upper = k;
470 for (
size_type idx=0; idx<num_i_parts; ++idx) {
471 size_type num_j_parts = i_tiles[idx].parts.size();
472 for (
size_type jdx=0; jdx<num_j_parts; ++jdx) {
473 size_type lower = i_tiles[idx].parts[jdx].parts[0].lower;
474 size_type upper = i_tiles[idx].parts[jdx].parts[0].upper;
476 size_type num_k_parts = (range + j_tile_size-1) / j_tile_size;
480 Array<KTile> k_tiles(num_k_parts);
481 for (
size_type k=0; k<num_k_parts; ++k) {
482 k_tiles[k].lower = lower + k*kts;
483 k_tiles[k].upper = std::min(upper+1, lower + (k+1)*kts);
485 size_type num_k = i_tiles[idx].parts[jdx].parts[0].parts.size();
487 size_type i = i_tiles[idx].parts[jdx].parts[0].parts[l].i;
488 size_type j = i_tiles[idx].parts[jdx].parts[0].parts[l].j;
489 size_type k = i_tiles[idx].parts[jdx].parts[0].parts[l].k;
490 value_type cijk = i_tiles[idx].parts[jdx].parts[0].parts[l].cijk;
494 while (kdx < num_k_parts && k >= k_tiles[kdx].lower) ++kdx;
496 TEUCHOS_ASSERT(kdx >= 0 && kdx < num_k_parts);
500 k_tiles[kdx].parts.push_back(
coord);
502 if (
j != k) ++num_coord;
506 Array<KTile> k_tiles2;
507 for (
size_type k=0; k<num_k_parts; ++k) {
508 if (k_tiles[k].parts.size() > 0)
509 k_tiles2.push_back(k_tiles[k]);
511 i_tiles[idx].parts[jdx].parts.swap(k_tiles2);
518 size_type max_num_j_parts = 0, max_num_k_parts = 0;
519 Array< Array< Array< Array<size_type> > > > coord_work(num_i_parts);
520 for (
size_type idx=0; idx<num_i_parts; ++idx) {
521 size_type num_j_parts = i_tiles[idx].parts.size();
522 max_num_j_parts = std::max(max_num_j_parts, num_j_parts);
523 coord_work[idx].resize(num_j_parts);
524 for (
size_type jdx=0; jdx<num_j_parts; ++jdx) {
525 size_type num_k_parts = i_tiles[idx].parts[jdx].parts.size();
526 max_num_k_parts = std::max(max_num_k_parts, num_k_parts);
527 coord_work[idx][jdx].resize(num_k_parts);
528 for (
size_type kdx=0; kdx<num_k_parts; ++kdx) {
529 size_type num_rows = i_tiles[idx].upper - i_tiles[idx].lower + 1;
530 total_num_rows += num_rows;
531 max_num_rows = std::max(max_num_rows, num_rows);
532 coord_work[idx][jdx][kdx].resize(num_rows, 0);
534 size_type nc = i_tiles[idx].parts[jdx].parts[kdx].parts.size();
536 size_type i = i_tiles[idx].parts[jdx].parts[kdx].parts[c].i;
538 ++(coord_work[idx][jdx][kdx][i-
i_begin]);
546 for (
size_type idx=0; idx<num_i_parts; ++idx) {
547 size_type num_j_parts = i_tiles[idx].parts.size();
548 for (
size_type jdx=0; jdx<num_j_parts; ++jdx) {
549 size_type num_k_parts = i_tiles[idx].parts[jdx].parts.size();
550 for (
size_type kdx=0; kdx<num_k_parts; ++kdx) {
551 size_type sz = coord_work[idx][jdx][kdx].size();
553 const size_t rem = coord_work[idx][jdx][kdx][i] %
tensor_align;
556 coord_work[idx][jdx][kdx][i] += pad;
580 max_num_j_parts, max_num_k_parts,
583 max_num_j_parts, max_num_k_parts,
591 typename value_array_type::HostMirror host_value =
593 typename coord_array_type::HostMirror host_coord =
595 typename coord2_array_type::HostMirror host_coord2 =
597 typename i_begin_type::HostMirror host_i_begin =
599 typename i_size_type::HostMirror host_i_size =
601 typename num_j_type::HostMirror host_num_j =
603 typename j_begin_type::HostMirror host_j_begin =
605 typename j_size_type::HostMirror host_j_size =
607 typename num_k_type::HostMirror host_num_k =
609 typename k_begin_type::HostMirror host_k_begin =
611 typename k_size_type::HostMirror host_k_size =
613 typename row_map_type::HostMirror host_row_map =
615 typename num_entry_type::HostMirror host_num_entry =
620 for (
size_type idx=0; idx<num_i_parts; ++idx) {
621 size_type num_j_parts = i_tiles[idx].parts.size();
622 for (
size_type jdx=0; jdx<num_j_parts; ++jdx) {
623 size_type num_k_parts = i_tiles[idx].parts[jdx].parts.size();
624 for (
size_type kdx=0; kdx<num_k_parts; ++kdx) {
625 size_type nc = coord_work[idx][jdx][kdx].size();
626 host_row_map(idx,jdx,kdx,0) = sum;
628 sum += coord_work[idx][jdx][kdx][t];
629 host_row_map(idx,jdx,kdx,t+1) = sum;
630 host_num_entry(idx,jdx,kdx,t) = 0;
637 for (
size_type idx=0; idx<num_i_parts; ++idx) {
638 size_type num_j_parts = i_tiles[idx].parts.size();
639 for (
size_type jdx=0; jdx<num_j_parts; ++jdx) {
640 size_type num_k_parts = i_tiles[idx].parts[jdx].parts.size();
641 for (
size_type kdx=0; kdx<num_k_parts; ++kdx) {
642 size_type nc = coord_work[idx][jdx][kdx].size();
644 coord_work[idx][jdx][kdx][t] = host_row_map(idx,jdx,kdx,t);
651 for (
size_type idx=0; idx<num_i_parts; ++idx) {
652 host_i_begin(idx) = i_tiles[idx].lower;
653 host_i_size(idx) = i_tiles[idx].upper - i_tiles[idx].lower;
654 TEUCHOS_ASSERT(host_i_size(idx) <= i_tile_size);
655 size_type num_j_parts = i_tiles[idx].parts.size();
656 host_num_j(idx) = num_j_parts;
657 for (
size_type jdx=0; jdx<num_j_parts; ++jdx) {
658 host_j_begin(idx,jdx) = i_tiles[idx].parts[jdx].lower;
659 host_j_size(idx,jdx) = i_tiles[idx].parts[jdx].upper -
660 i_tiles[idx].parts[jdx].lower;
662 size_type num_k_parts = i_tiles[idx].parts[jdx].parts.size();
663 host_num_k(idx,jdx) = num_k_parts;
664 for (
size_type kdx=0; kdx<num_k_parts; ++kdx) {
665 host_k_begin(idx,jdx,kdx) = i_tiles[idx].parts[jdx].parts[kdx].lower;
666 host_k_size(idx,jdx,kdx) = i_tiles[idx].parts[jdx].parts[kdx].upper -
667 i_tiles[idx].parts[jdx].parts[kdx].lower;
670 size_type nc = i_tiles[idx].parts[jdx].parts[kdx].parts.size();
672 Coord s = i_tiles[idx].parts[jdx].parts[kdx].parts[t];
678 const size_type row = i - host_i_begin(idx);
679 const size_type n = coord_work[idx][jdx][kdx][row];
680 ++coord_work[idx][jdx][kdx][row];
682 host_value(n) = (
j != k) ? c : 0.5*c;
683 host_coord2(n,0) =
j - host_j_begin(idx,jdx);
684 host_coord2(n,1) = k - host_k_begin(idx,jdx,kdx);
685 host_coord(n) = (host_coord2(n,1) << 16) | host_coord2(n,0);
687 ++host_num_entry(idx,jdx,kdx,row);
710 for (
size_type idx=0; idx<num_i_parts; ++idx) {
711 size_type num_j_parts = i_tiles[idx].parts.size();
712 for (
size_type jdx=0; jdx<num_j_parts; ++jdx) {
713 size_type num_k_parts = i_tiles[idx].parts[jdx].parts.size();
714 for (
size_type kdx=0; kdx<num_k_parts; ++kdx) {
715 for (
size_type i = 0; i < host_i_size(idx); ++i) {
716 tensor.
m_flops += 5*host_num_entry(idx,jdx,kdx,i) + 1;