@@ -650,81 +650,75 @@ SSVectorBase<R>& SSVectorBase<R>::assign2productShort(const SVSetBase<S>& A,
650650{
651651 assert (x.isSetup ());
652652
653+ clear ();
654+
653655 if (x.size () == 0 ) // x can be setup but have size 0 => this := zero vector
654- {
655- clear ();
656656 return *this ;
657- }
658657
659658 // compute x[0] * A[0]
660659 int curidx = x.idx [0 ];
661660 const T x0 = x.val [curidx];
662661 const SVectorBase<S>& A0 = A[curidx];
663- int nonzero_idx = 0 ;
664662 int xsize = x.size ();
665663 int Aisize;
666664
667- num = A0.size ();
668-
669- if (isZero (x0, this ->tolerances ()->epsilon ()) || num == 0 )
665+ // If x[0] == 0, do nothing.
666+ if (isNotZero (x0, this ->tolerances ()->epsilon ()))
670667 {
671- // A[0] == 0 or x[0] == 0 => this := zero vector
672- clear ();
673- }
674- else
675- {
676- for (int j = 0 ; j < num; ++j)
668+ Aisize = A0.size ();
669+
670+ for (int j = 0 ; j < Aisize; ++j)
677671 {
678672 const Nonzero<S>& elt = A0.element (j);
679673 const R product = x0 * elt.val ;
680674
681- // store the value in any case
682- idx[nonzero_idx] = elt.idx ;
683- VectorBase<R>::val[elt.idx ] = product;
684-
685675 // count only non-zero values; not 'isNotZero(product, epsilon)'
686676 if (product != 0 )
687- ++nonzero_idx;
677+ {
678+ assert (num < len);
679+ idx[num++] = elt.idx ;
680+ VectorBase<R>::val[elt.idx ] = product;
681+ }
688682 }
689683 }
690684
691685 // Compute the other x[i] * A[i] and add them to the existing vector.
692686 for (int i = 1 ; i < xsize; ++i)
693687 {
694688 curidx = x.idx [i];
695- const T xi = x.val [curidx];
689+ const T xi = x.val [curidx];
696690 const SVectorBase<S>& Ai = A[curidx];
697691
698- // If A[i] == 0 or x[i] == 0, do nothing.
699- Aisize = Ai.size ();
700-
701- if (isNotZero (xi, this ->tolerances ()->epsilon ()) || Aisize == 0 )
692+ // If x[i] == 0, do nothing.
693+ if (isNotZero (xi, this ->tolerances ()->epsilon ()))
702694 {
695+ Aisize = Ai.size ();
696+
703697 // Compute x[i] * A[i] and add it to the existing vector.
704698 for (int j = 0 ; j < Aisize; ++j)
705699 {
706700 const Nonzero<S>& elt = Ai.element (j);
707- idx[nonzero_idx] = elt.idx ;
708- R oldval = VectorBase<R>::val[elt.idx ];
709-
710- // An old value of exactly 0 means the position is still unused.
711- // It will be used now (either by a new nonzero or by a SOPLEX_VECTOR_MARKER),
712- // so increase the counter. If oldval != 0, we just
713- // change an existing NZ-element, so don't increase the counter.
714- if (oldval == 0 )
715- ++nonzero_idx;
716-
717- // Add the current product x[i] * A[i][j]; if oldval was
718- // SOPLEX_VECTOR_MARKER before, it does not hurt because SOPLEX_VECTOR_MARKER is really small.
719- oldval += xi * elt.val ;
720-
721- // If the new value is exactly 0, mark the index as used
722- // by setting a value which is nearly 0; otherwise, store
723- // the value. Values below epsilon will be removed later.
724- if (oldval == 0 )
701+ const R oldval = VectorBase<R>::val[elt.idx ];
702+ const R newval = oldval + xi * elt.val ;
703+
704+ // If the value becomes exactly 0, mark the index as used by setting a value which is
705+ // nearly 0. Values below epsilon will be removed later.
706+ if (oldval != 0 && newval == 0 )
725707 VectorBase<R>::val[elt.idx ] = SOPLEX_VECTOR_MARKER;
708+ // Otherwise, store the value. If oldval was SOPLEX_VECTOR_MARKER before, it does not
709+ // hurt because SOPLEX_VECTOR_MARKER is really small.
726710 else
727- VectorBase<R>::val[elt.idx ] = oldval;
711+ {
712+ VectorBase<R>::val[elt.idx ] = newval;
713+
714+ // If the value changes from exactly 0, the position is still unused, so increase
715+ // the counter. Otherwise, we just change the value of an existing element.
716+ if (oldval == 0 && newval != 0 )
717+ {
718+ assert (num < len);
719+ idx[num++] = elt.idx ;
720+ }
721+ }
728722 }
729723 }
730724 }
@@ -733,21 +727,18 @@ SSVectorBase<R>& SSVectorBase<R>::assign2productShort(const SVSetBase<S>& A,
733727 // zeroing all values which are nearly 0, and setting #num# appropriately.
734728 int nz_counter = 0 ;
735729
736- for (int i = 0 ; i < nonzero_idx ; ++i)
730+ for (int i = 0 ; i < num ; ++i)
737731 {
738732 curidx = idx[i];
739733
740734 if (isZero (VectorBase<R>::val[curidx], this ->tolerances ()->epsilon ()))
741735 VectorBase<R>::val[curidx] = 0 ;
742736 else
743- {
744- idx[nz_counter] = curidx;
745- ++nz_counter;
746- }
747-
748- num = nz_counter;
737+ idx[nz_counter++] = curidx;
749738 }
750739
740+ num = nz_counter;
741+
751742 assert (isConsistent ());
752743
753744 return *this ;
0 commit comments