1 /* ========================================================================== */
2 /* === Source/Mongoose_QPBoundary.cpp ======================================= */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * Mongoose Graph Partitioning Library Copyright (C) 2017-2018,
7 * Scott P. Kolodziej, Nuri S. Yeralan, Timothy A. Davis, William W. Hager
8 * Mongoose is licensed under Version 3 of the GNU General Public License.
9 * Mongoose is also available under other licenses; contact authors for details.
10 * -------------------------------------------------------------------------- */
11
12 /*
13 Move all components of x to boundary of the feasible region
14
15 0 <= x <= 1, a'x = b, lo <= b <= hi
16
17 while decreasing the cost function. The algorithm has the following parts
18
19 1. For each i in the free set, see if x_i can be feasibly pushed to either
20 boundary while decreasing the cost.
21
22 2. For each i in the bound set, see if x_i can be feasibly flipped to
23 opposite boundary while decreasing the cost.
24
25 3. For each i in the free list with a_{ij} = 0 and with j free,
26 move either x_i or x_j to the boundary while decreasing
27 the cost. The adjustments has the form x_i = s/a_i and x_j = -s/a_j
28 where s is a scalar factor. These adjustments must decrease cost.
29
30 4. For the remaining i in the free list, take pair x_i and x_j and
31 apply adjustments of the same form as in #2 above to push at least one
32 component to boundary. The quadratic terms can only decrease the
33 cost function. We choose the sign of s such that g_i x_i + g_j x_j <= 0.
34 Hence, this adjustment cannot increase the cost.
35 */
36
37 /* ========================================================================== */
38
39 #include "Mongoose_QPBoundary.hpp"
40 #include "Mongoose_Debug.hpp"
41 #include "Mongoose_Internal.hpp"
42
43 #define EMPTY (-1)
44
45 namespace Mongoose
46 {
47
QPBoundary(EdgeCutProblem * graph,const EdgeCut_Options * options,QPDelta * QP)48 void QPBoundary(EdgeCutProblem *graph, const EdgeCut_Options *options, QPDelta *QP)
49 {
50 (void)options; // Unused variable
51 /* ---------------------------------------------------------------------- */
52 /* Step 0. read in the needed arrays */
53 /* ---------------------------------------------------------------------- */
54
55 /* input and output */
56
57 //--- FreeSet
58 Int nFreeSet = QP->nFreeSet;
59 Int *FreeSet_list = QP->FreeSet_list; /* list for free indices */
60 Int *FreeSet_status = QP->FreeSet_status;
61 /* FreeSet_status [i] = +1, -1, or 0
62 if x_i = 1, 0, or 0 < x_i < 1 */
63 //---
64
65 PR(("Mongoose_QPBoundary nFreeSet %ld\n", nFreeSet));
66
67 if (nFreeSet == 0)
68 {
69 // quick return if FreeSet is empty
70 return;
71 }
72
73 double *x = QP->x; /* current estimate of solution */
74 double *grad = QP->gradient; /* gradient at current x */
75 Int ib = QP->ib; /* ib = +1, -1, or 0 ,
76 if b = hi, lo, or lo < b < hi, respectively. Note there are cases
77 where roundoff occurs, and ib can be zero even though b == lo or
78 b == hi. The value of be can even be < lo or > hi, but only by a tiny
79 amount of roundoff error. This is OK. */
80
81 double b = QP->b; /* current value for a'x */
82
83 /* problem specification for the graph G */
84 Int n = graph->n; /* problem dimension */
85 double *Ex = graph->x; /* numerical values for edge weights */
86 Int *Ei = graph->i; /* adjacent vertices for each vertex */
87 Int *Ep = graph->p; /* points into Ex or Ei */
88 double *a = graph->w; /* a'x = b, lo <= b <= hi */
89
90 double lo = QP->lo;
91 double hi = QP->hi;
92
93 /* work array */
94 double *D = QP->D; /* diagonal of quadratic */
95
96 PR(("\n----- QPBoundary start: [\n"));
97 DEBUG(QPcheckCom(graph, options, QP, 1, QP->nFreeSet, QP->b)); // check b
98
99 /* ---------------------------------------------------------------------- */
100 /* Step 1. if lo < b < hi, then for each free k, */
101 /* see if x_k can be pushed to 0 or 1 */
102 /* ---------------------------------------------------------------------- */
103
104 DEBUG(FreeSet_dump("QPBoundary start", n, FreeSet_list, nFreeSet,
105 FreeSet_status, 0, x));
106
107 PR(("Boundary 1 start: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo,
108 b, hi, b - lo, hi - b));
109
110 Int kfree2 = 0;
111 for (Int kfree = 0; kfree < nFreeSet; kfree++)
112 {
113 // Once b becomes bounded, the remainder of the FreeSet is unchanged,
114 // and no further changes are made to x. However, this loop must still
115 // continue, so as to compact the FreeSet from deletions made by earlier
116 // iterations.
117
118 // get the next k from the FreeSet
119 Int k = FreeSet_list[kfree];
120
121 PR(("Step 1: k %ld x[k] %g ib %ld b %g\n", k, x[k], ib, b));
122
123 // only modify x[k] if ib == 0 (which means lo < b < hi)
124 if (ib == 0)
125 {
126 double delta_xk;
127 double ak = (a) ? a[k] : 1;
128 if (grad[k] > 0.0)
129 {
130 // decrease x [k]
131 delta_xk = (b - lo) / ak; // note that delta_xk > 0
132 if (delta_xk < x[k])
133 {
134 // x [k] decreases by delta_xk but does not hit zero
135 // b hits the lower bound, lo
136 ib = -1;
137 b = lo;
138 x[k] -= delta_xk;
139 //--- keep k in the FreeSet
140 FreeSet_list[kfree2++] = k;
141 }
142 else
143 {
144 // x [k] hits lower bound of zero
145 // b does not hit lo; still between lower and upper bound
146 delta_xk = x[k];
147 x[k] = 0.;
148 FreeSet_status[k] = -1;
149 b -= delta_xk * ak;
150 //--- remove k from the FreeSet by not incrementing kfree2
151 }
152 }
153 else
154 {
155 // increase x [k]
156 delta_xk = (b - hi) / ak; // note that delta_xk < 0
157 if (delta_xk < x[k] - 1.)
158 {
159 // x [k] hits upper bound of one
160 // b does not reach hi; still between lower and upper bound
161 delta_xk = x[k] - 1.;
162 x[k] = 1.;
163 FreeSet_status[k] = +1;
164 b -= delta_xk * ak;
165 //--- remove k from the FreeSet by not incrementing kfree2
166 }
167 else
168 {
169 // x [k] increases by -delta_xk but does not hit one
170 // b hits the upper bound, hi.
171 ib = +1;
172 b = hi;
173 x[k] -= delta_xk;
174 //--- keep k in the FreeSet
175 FreeSet_list[kfree2++] = k;
176 }
177 }
178 // x [k] has dropped by delta_xk, so update the gradient
179 for (Int p = Ep[k]; p < Ep[k + 1]; p++)
180 {
181 grad[Ei[p]] += delta_xk * ((Ex) ? Ex[p] : 1);
182 }
183 grad[k] += delta_xk * D[k];
184 }
185 else
186 {
187 // b is at lo or hi and thus x [k] is not changed.
188 // Once this happens, the remainder of this loop does this next
189 // step only, and no further changes are made to x and the FreeSet.
190 //--- keep k in the FreeSet
191 FreeSet_list[kfree2++] = k;
192 }
193 }
194
195 // update the size of the FreeSet, after pruning
196 nFreeSet = kfree2;
197
198 /* ---------------------------------------------------------------------- */
199 /* Step 2. Examine flips of x_k from 0 to 1 or from 1 to 0 */
200 /* ---------------------------------------------------------------------- */
201
202 PR(("Boundary step 2:\n"));
203
204 for (Int k = 0; k < n; k++)
205 {
206 Int FreeSet_status_k = FreeSet_status[k];
207 if (FreeSet_status_k == 0)
208 {
209 // k is in FreeSet so it cannot be simply flipped 0->1 or 1->0
210 continue;
211 }
212
213 // k not in FreeSet, so no changes here to FreeSet
214
215 double ak = (a) ? a[k] : 1;
216 if (FreeSet_status_k > 0) /* try changing x_k from 1 to 0 */
217 {
218 if (b - ak >= lo)
219 {
220 if (0.5 * D[k] + grad[k] >= 0) /* flip lowers cost */
221 {
222 b -= ak;
223 ib = (b <= lo ? -1 : 0);
224 x[k] = 0.0;
225 FreeSet_status[k] = -1;
226 }
227 }
228 }
229 else /* try changing x_k from 0 to 1 */
230 {
231 if (b + ak <= hi)
232 {
233 if (grad[k] - 0.5 * D[k] <= 0) /* flip lowers cost */
234 {
235 b += ak;
236 ib = (b >= hi ? 1 : 0);
237 x[k] = 1.0;
238 FreeSet_status[k] = +1;
239 }
240 }
241 }
242
243 if (FreeSet_status_k != FreeSet_status[k])
244 {
245 if (FreeSet_status_k == 1) /* x [k] was 1, now it is 0 */
246 {
247 for (Int p = Ep[k]; p < Ep[k + 1]; p++)
248 {
249 grad[Ei[p]] += (Ex) ? Ex[p] : 1;
250 }
251 grad[k] += D[k];
252 }
253 else /* x [k] was 0, now it is 1 */
254 {
255 for (Int p = Ep[k]; p < Ep[k + 1]; p++)
256 {
257 grad[Ei[p]] -= (Ex) ? Ex[p] : 1;
258 }
259 grad[k] -= D[k];
260 }
261 }
262 // DEBUG (QPcheckCom (graph, options, QP, 1, nFreeSet, b)) ; //
263 // check b
264 }
265
266 /* ---------------------------------------------------------------------- */
267 // quick return if FreeSet is now empty
268 /* ---------------------------------------------------------------------- */
269
270 if (nFreeSet == 0)
271 {
272 PR(("Boundary quick: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo,
273 b, hi, b - lo, hi - b));
274 QP->nFreeSet = nFreeSet;
275 QP->b = b;
276 QP->ib = ib;
277 PR(("------- QPBoundary end ]\n"));
278 return;
279 }
280
281 /* ---------------------------------------------------------------------- */
282 /* Step 3. Search for a_{ij} = 0 in the free index set */
283 /* ---------------------------------------------------------------------- */
284
285 // look for where both i and j are in the FreeSet,
286 // but i and j are not adjacent in the graph G.
287
288 DEBUG(FreeSet_dump("step 3", n, FreeSet_list, nFreeSet, FreeSet_status, 0,
289 x));
290
291 // for each j in FreeSet, except for the last one
292 for (Int jfree = 0; jfree < nFreeSet - 1; jfree++)
293 {
294
295 // get j from the FreeSet
296 Int j = FreeSet_list[jfree];
297 if (j == EMPTY)
298 {
299 // j has already been deleted, skip it
300 continue;
301 }
302
303 /* -------------------------------------------------------------- */
304 /* find i and j both free and where a_{ij} = 0 */
305 /* -------------------------------------------------------------- */
306
307 // mark all vertices i adjacent to j in the FreeSet
308 for (Int p = Ep[j]; p < Ep[j + 1]; p++)
309 {
310 Int i = Ei[p];
311 ASSERT(i != j); // graph has no self edges
312 graph->mark(i);
313 }
314 graph->mark(j);
315
316 // for each i that follows after j in the FreeSet
317 for (Int ifree = jfree + 1; ifree < nFreeSet; ifree++)
318 {
319
320 // get i from the FreeSet
321 Int i = FreeSet_list[ifree];
322 if (i == EMPTY)
323 {
324 // i has already been deleted it; skip it
325 continue;
326 }
327
328 if (!graph->isMarked(i))
329 {
330 // vertex i is not adjacent to j in the graph G
331 double aj = (a) ? a[j] : 1;
332 double ai = (a) ? a[i] : 1;
333 double xi = x[i];
334 double xj = x[j];
335
336 /* cost change if x_j increases dx_j = s/a_j, dx_i = s/a_i */
337 double s;
338 Int bind1, bind2;
339 if (aj * (1. - xj) < ai * xi) // x_j hits upper bound
340 {
341 s = aj * (1. - xj);
342 bind1 = 1;
343 }
344 else /* x_i hits lower bound */
345 {
346 s = ai * xi;
347 bind1 = 0;
348 }
349 double dxj = s / aj;
350 double dxi = -s / ai;
351 double c1 = (grad[j] - .5 * D[j] * dxj) * dxj
352 + (grad[i] - .5 * D[i] * dxi) * dxi;
353
354 /* cost change if x_j decreases dx_j = s/a_j, dx_i = s/a_i */
355 if (aj * xj < ai * (1. - xi)) // x_j hits lower bound
356 {
357 s = -aj * xj;
358 bind2 = -1;
359 }
360 else /* x_i hits upper bound */
361 {
362 s = -ai * (1. - xi);
363 bind2 = 0;
364 }
365 dxj = s / aj;
366 dxi = -s / ai;
367 double c2 = (grad[j] - 0.5 * D[j] * dxj) * dxj
368 + (grad[i] - 0.5 * D[i] * dxi) * dxi;
369
370 Int new_FreeSet_status;
371 if (c1 < c2) /* increase x_j */
372 {
373 if (bind1 == 1)
374 {
375 // j is bound (not i) and x_j becomes 1
376 dxj = 1. - xj;
377 dxi = -aj * dxj / ai;
378 x[j] = 1.;
379 x[i] += dxi;
380 new_FreeSet_status = +1; /* j is bound at 1 */
381 }
382 else // bind1 is zero
383 {
384 // i is bound (not j) and x_i becomes 0
385 dxi = -xi;
386 dxj = -ai * dxi / aj;
387 x[i] = 0.;
388 x[j] += dxj;
389 new_FreeSet_status = -1; /* i is bound at 0 */
390 }
391 }
392 else
393 {
394 if (bind2 == -1)
395 {
396 // j is bound (not i) and x_j becomes 0
397 bind1 = 1;
398 x[j] = 0.;
399 x[i] += dxi;
400 new_FreeSet_status = -1; /* j is bound at 0 */
401 }
402 else /* x_i = 1 */
403 {
404 // i is bound (not j) and x_i becomes 1
405 bind1 = 0;
406 x[i] = 1;
407 x[j] += dxj;
408 new_FreeSet_status = +1; /* i is bound at 1 */
409 }
410 }
411
412 for (Int p = Ep[j]; p < Ep[j + 1]; p++)
413 {
414 grad[Ei[p]] -= ((Ex) ? Ex[p] : 1) * dxj;
415 }
416 for (Int p = Ep[i]; p < Ep[i + 1]; p++)
417 {
418 grad[Ei[p]] -= ((Ex) ? Ex[p] : 1) * dxi;
419 }
420 grad[j] -= D[j] * dxj;
421 grad[i] -= D[i] * dxi;
422
423 // Remove either i or j from the FreeSet. Note that it
424 // is possible for both x[i] and x[j] to reach their bounds
425 // at the same time. Only one is removed from the FreeSet;
426 // the other will be removed later.
427
428 if (bind1)
429 {
430 // remove j from the FreeSet by setting its place to EMPTY
431 PR(("(b1):remove j = %ld from the FreeSet\n", j));
432 ASSERT(j == FreeSet_list[jfree]);
433 ASSERT(FreeSet_status[j] == 0);
434 FreeSet_list[jfree] = EMPTY;
435 FreeSet_status[j] = new_FreeSet_status;
436 ASSERT(FreeSet_status[j] != 0);
437 //---
438 // no longer consider j, so skip all of remainder of i loop
439 break;
440 }
441 else
442 {
443 // remove i from the FreeSet by setting its place to EMPTY
444 PR(("(b2):remove i = %ld from the FreeSet\n", i));
445 ASSERT(i == FreeSet_list[ifree]);
446 ASSERT(FreeSet_status[i] == 0);
447 FreeSet_list[ifree] = EMPTY;
448 FreeSet_status[i] = new_FreeSet_status;
449 ASSERT(FreeSet_status[i] != 0);
450 //---
451 // keep j, and consider it with the next i
452 continue;
453 }
454 }
455 }
456
457 // clear the marks from all the vertices
458 graph->clearMarkArray();
459 }
460
461 // remove deleted vertices from the FreeSet
462 kfree2 = 0;
463 for (Int kfree = 0; kfree < nFreeSet; kfree++)
464 {
465 Int k = FreeSet_list[kfree];
466 if (k != EMPTY)
467 {
468 // keep k in the FreeSet
469 FreeSet_list[kfree2++] = k;
470 ASSERT(0 <= k && k < n);
471 ASSERT(FreeSet_status[k] == 0);
472 }
473 }
474 nFreeSet = kfree2;
475
476 DEBUG(FreeSet_dump("step 3 done", n, FreeSet_list, nFreeSet, FreeSet_status,
477 1, x));
478
479 DEBUG(QPcheckCom(graph, options, QP, 1, nFreeSet, b)); // check b
480
481 #ifndef NDEBUG
482 // the vertices in the FreeSet now form a single clique. Check this.
483 // this test is for debug mode only
484 ASSERT(nFreeSet >= 1); // we can have 1 or more vertices still in FreeSet
485 for (Int kfree = 0; kfree < nFreeSet; kfree++)
486 {
487 // j must be adjacent to all other vertices in the FreeSet
488 Int j = FreeSet_list[kfree];
489 Int nfree_neighbors = 0;
490 for (Int p = Ep[j]; p < Ep[j + 1]; p++)
491 {
492 Int i = Ei[p];
493 ASSERT(i != j);
494 if (FreeSet_status[i] == 0)
495 nfree_neighbors++;
496 }
497 ASSERT(nfree_neighbors == nFreeSet - 1);
498 }
499 #endif
500
501 /* ---------------------------------------------------------------------- */
502 /* Step 4. dxj = s/aj, dxi = -s/ai, choose s with g_j dxj + g_i dxi <= 0 */
503 /* ---------------------------------------------------------------------- */
504
505 DEBUG(FreeSet_dump("step 4 starts", n, FreeSet_list, nFreeSet,
506 FreeSet_status, 0, x));
507
508 // consider pairs of vertices in the FreeSet, until only one is left
509 while (nFreeSet > 1)
510 {
511 /* free variables: 0 < x_j < 1 */
512 /* choose s so that first derivative terms decrease */
513
514 // i and j are the last two vertex in the FreeSet_list, as in:
515 // FreeSet_list = [ .... i j ]
516 // at the end of this iteration, one will be deleted, thus becoming
517 // FreeSet_list = [ .... j ]
518 // or
519 // FreeSet_list = [ .... i ]
520 Int j = FreeSet_list[nFreeSet - 1];
521 ASSERT(FreeSet_status[j] == 0);
522 Int i = FreeSet_list[nFreeSet - 2];
523 ASSERT(FreeSet_status[i] == 0);
524
525 double ai = (a) ? a[i] : 1;
526 double aj = (a) ? a[j] : 1;
527 double xi = x[i];
528 double xj = x[j];
529
530 Int new_FreeSet_status;
531 Int bind1;
532 double dxj, dxi, s = grad[j] / aj - grad[i] / ai;
533
534 if (s < 0.) /* increase x_j */
535 {
536 if (aj * (1. - xj) < ai * xi) /* x_j hits upper bound */
537 {
538 dxj = 1. - xj;
539 dxi = -aj * dxj / ai;
540 x[j] = 1.;
541 x[i] += dxi;
542 new_FreeSet_status = +1;
543 bind1 = 1; /* x_j is bound at 1 */
544 }
545 else /* x_i hits lower bound */
546 {
547 dxi = -xi;
548 dxj = -ai * dxi / aj;
549 x[i] = 0.;
550 x[j] += dxj;
551 new_FreeSet_status = -1;
552 bind1 = 0; /* x_i is bound at 0 */
553 }
554 }
555 else /* decrease x_j */
556 {
557 if (aj * xj < ai * (1. - xi)) /* x_j hits lower bound */
558 {
559 dxj = -xj;
560 dxi = -aj * dxj / ai;
561 x[j] = 0;
562 x[i] += dxi;
563 new_FreeSet_status = -1;
564 bind1 = 1; /* x_j is bound */
565 }
566 else /* x_i hits upper bound */
567 {
568 dxi = 1 - xi;
569 dxj = -ai * dxi / aj;
570 x[i] = 1;
571 x[j] += dxj;
572 new_FreeSet_status = +1;
573 bind1 = 0; /* x_i is bound */
574 }
575 }
576
577 for (Int k = Ep[j]; k < Ep[j + 1]; k++)
578 {
579 grad[Ei[k]] -= ((Ex) ? Ex[k] : 1) * dxj;
580 }
581 for (Int k = Ep[i]; k < Ep[i + 1]; k++)
582 {
583 grad[Ei[k]] -= ((Ex) ? Ex[k] : 1) * dxi;
584 }
585 grad[j] -= D[j] * dxj;
586 grad[i] -= D[i] * dxi;
587
588 // ---------------------------------------------------------------------
589 // the following 2 cases define the next j in the iteration:
590 // ---------------------------------------------------------------------
591
592 // Remove either i or j from the FreeSet. Note that it is possible for
593 // both x[i] and x[j] to reach their bounds at the same time. Only one
594 // is removed from the FreeSet; the other will be removed later.
595
596 if (bind1)
597 {
598 // j is bound.
599 // remove j from the FreeSet, and keep i. The FreeSet_list was
600 // FreeSet_list = [ .... i j ] becomes FreeSet_list = [ .... i ]
601 PR(("(b3):remove j = %ld from the FreeSet\n", j));
602 ASSERT(FreeSet_status[j] == 0);
603 FreeSet_status[j] = new_FreeSet_status;
604 ASSERT(FreeSet_status[j] != 0);
605 }
606 else
607 {
608 // i is bound.
609 // remove i from the FreeSet, and keep j. The FreeSet_list was
610 // FreeSet_list = [ .... i j ] becomes FreeSet_list = [ .... j ]
611 PR(("(b4):remove i = %ld from the FreeSet\n", i));
612 ASSERT(FreeSet_status[i] == 0);
613 FreeSet_status[i] = new_FreeSet_status;
614 ASSERT(FreeSet_status[i] != 0);
615 // shift j down by one in the list, thus discarding j.
616 ASSERT(FreeSet_list[nFreeSet - 2] == i);
617 FreeSet_list[nFreeSet - 2] = j;
618 }
619
620 // one fewer vertex in the FreeSet (i or j removed)
621 nFreeSet--;
622
623 DEBUG(FreeSet_dump("step 4", n, FreeSet_list, nFreeSet, FreeSet_status,
624 0, x));
625 DEBUG(QPcheckCom(graph, options, QP, 1, nFreeSet, b)); // check b
626 }
627
628 DEBUG(FreeSet_dump("wrapup", n, FreeSet_list, nFreeSet, FreeSet_status, 0,
629 x));
630
631 /* ---------------------------------------------------------------------- */
632 /* step 5: at most one free variable remaining */
633 /* ---------------------------------------------------------------------- */
634
635 ASSERT(nFreeSet == 0 || nFreeSet == 1);
636
637 PR(("Step 5: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo, b, hi,
638 b - lo, hi - b));
639
640 if (nFreeSet == 1) /* j is free, optimize over x [j] */
641 {
642 // j is the first and only item in the FreeSet
643 Int j = FreeSet_list[0];
644 PR(("ONE AND ONLY!! j = %ld x[j] %g\n", j, x[j]));
645
646 Int bind1 = 0;
647 double aj = (a) ? a[j] : 1;
648 double dxj = (hi - b) / aj;
649 PR(("dxj %g x[j] %g (1-x[j]): %g\n", dxj, x[j], 1 - x[j]));
650 if (dxj < 1. - x[j])
651 {
652 bind1 = 1;
653 }
654 else
655 {
656 dxj = 1. - x[j];
657 }
658
659 Int bind2 = 0;
660 double dxi = (lo - b) / aj;
661 PR(("dxi %g x[j] %g (-x[j]): %g\n", dxi, x[j], -x[j]));
662 if (dxi > -x[j])
663 {
664 bind2 = 1;
665 }
666 else
667 {
668 dxi = -x[j];
669 }
670
671 double c1 = (grad[j] - 0.5 * D[j] * dxj) * dxj;
672 double c2 = (grad[j] - 0.5 * D[j] * dxi) * dxi;
673 if (c1 <= c2) /* x [j] += dxj */
674 {
675 if (bind1)
676 {
677 PR(("bind1: xj changes from %g", x[j]));
678 x[j] += dxj;
679 PR((" to %g, b now at hi\n", x[j]));
680 ib = +1;
681 b = hi;
682 }
683 else
684 {
685 x[j] = 1.;
686 b += dxj * aj;
687 /// remove j from the FreeSet, which is now empty
688 PR(("(b5):remove j = %ld from FreeSet, now empty\n", j));
689 ASSERT(FreeSet_status[j] == 0);
690 FreeSet_status[j] = 1;
691 ASSERT(FreeSet_status[j] != 0);
692 nFreeSet--;
693 ASSERT(nFreeSet == 0);
694 }
695 }
696 else /* x [j] += dxi */
697 {
698 dxj = dxi;
699 if (bind2)
700 {
701 PR(("bind2: xj changes from %g", x[j]));
702 x[j] += dxj;
703 PR((" to %g, b now at lo\n", x[j]));
704 ib = -1;
705 b = lo;
706 }
707 else
708 {
709 x[j] = 0.;
710 b += dxj * aj;
711 /// remove j from the FreeSet, which is now empty
712 PR(("(b6):remove j = %ld from FreeSet, now empty\n", j));
713 ASSERT(FreeSet_status[j] == 0);
714 FreeSet_status[j] = -1;
715 ASSERT(FreeSet_status[j] != 0);
716 nFreeSet--;
717 ASSERT(nFreeSet == 0);
718 }
719 }
720
721 if (dxj != 0.)
722 {
723 for (Int p = Ep[j]; p < Ep[j + 1]; p++)
724 {
725 grad[Ei[p]] -= ((Ex) ? Ex[p] : 1) * dxj;
726 }
727 grad[j] -= D[j] * dxj;
728 }
729 }
730
731 /* ---------------------------------------------------------------------- */
732 // wrapup
733 /* ---------------------------------------------------------------------- */
734
735 PR(("QBboundary, done:\n"));
736 DEBUG(FreeSet_dump("QPBoundary: done ", n, FreeSet_list, nFreeSet,
737 FreeSet_status, 0, x));
738 ASSERT(nFreeSet == 0 || nFreeSet == 1);
739 PR(("Boundary done: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo, b,
740 hi, b - lo, hi - b));
741
742 QP->nFreeSet = nFreeSet;
743 QP->b = b;
744 QP->ib = ib;
745
746 // clear the marks from all the vertices
747 graph->clearMarkArray();
748
749 DEBUG(QPcheckCom(graph, options, QP, 1, nFreeSet, b)); // check b
750 PR(("----- QPBoundary end ]\n"));
751 }
752
753 } // end namespace Mongoose
754