1 /* ========================================================================== */
2 /* === Source/Mongoose_QPNapDown.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 /* === QPNapDown ============================================================ */
14 /* ========================================================================== */
15
16 /* Find x that minimizes ||x-y|| while satisfying the constraints
17 0 <= x <= 1, a'x = b.
18 The algorithm is described in the napsack comments.
19 It is assumed that the starting guess lambda for the dual multiplier is >=
20 the correct multiplier. Hence, lambda will be decreased. The slope of the
21 dual function, neglecting b, starts out smaller than b. We stop
22 when we reach b. We assume that a >= 0, so that as lambda decreases,
23 x_i (lambda) increases. Hence, the only bound variables that can become
24 free are those with x_i (lambda) <= 0 */
25
26 #include "Mongoose_QPNapDown.hpp"
27 #include "Mongoose_Debug.hpp"
28 #include "Mongoose_Internal.hpp"
29 #include "Mongoose_Logger.hpp"
30 #include "Mongoose_QPMaxHeap.hpp"
31
32 namespace Mongoose
33 {
34
QPNapDown(const double * x,const Int n,double lambda,const double * a,double b,double * breakpts,Int * bound_heap,Int * free_heap)35 double QPNapDown /* return lambda */
36 (const double *x, /* holds y on input, not modified */
37 const Int n, /* size of x */
38 double lambda, /* initial guess for the shift */
39 const double *a, /* input constraint vector */
40 double b, /* input constraint scalar */
41 double *breakpts, /* break points */
42 Int *bound_heap, /* work array */
43 Int *free_heap /* work array */
44 )
45 {
46 Int i, k, e, maxsteps, n_bound, n_free;
47 double ai, asum, a2sum, maxbound, maxfree, t;
48
49 maxbound = -INFINITY;
50 maxfree = -INFINITY;
51
52 /* -------------------------------------------------------------- */
53 /* construct the heaps */
54 /* -------------------------------------------------------------- */
55
56 n_bound = 0;
57 n_free = 0;
58 asum = 0.;
59 a2sum = 0.;
60
61 for (i = 0; i < n; i++)
62 {
63 ai = (a) ? a[i] : 1;
64 double xi = x[i] - ai * lambda;
65 if (xi < 0.)
66 {
67 n_bound++;
68 bound_heap[n_bound] = i;
69
70 t = x[i] / ai;
71 maxbound = std::max(maxbound, t);
72 breakpts[i] = t;
73 }
74 else if (xi < 1.)
75 {
76 n_free++;
77 free_heap[n_free] = i;
78 asum += x[i] * ai;
79 a2sum += ai * ai;
80 t = (x[i] - 1.) / ai;
81 maxfree = std::max(maxfree, t);
82 breakpts[i] = t;
83 }
84 else
85 {
86 asum += ai;
87 }
88 }
89
90 /*------------------------------------------------------------------- */
91 /* check to see if zero slope achieved without changing the free set */
92 /* remember that the slope must always be adjusted by b */
93 /*------------------------------------------------------------------- */
94
95 maxsteps = 2 * n + 1;
96 for (k = 1; k <= maxsteps; k++)
97 {
98 double new_break = std::max(maxfree, maxbound);
99 double s = asum - new_break * a2sum;
100 if ((s >= b) || (new_break == -INFINITY)) /* done */
101 {
102 if (a2sum != 0.)
103 {
104 lambda = (asum - b) / a2sum;
105 }
106 return lambda;
107 }
108 lambda = new_break;
109
110 if (k == 1)
111 {
112 QPMaxHeap_build(free_heap, n_free, breakpts);
113 QPMaxHeap_build(bound_heap, n_bound, breakpts);
114 }
115
116 /* -------------------------------------------------------------- */
117 /* update the heaps */
118 /* -------------------------------------------------------------- */
119
120 if (n_free > 0)
121 {
122 while (breakpts[e = free_heap[1]] >= lambda)
123 {
124 ai = (a) ? a[e] : 1;
125 a2sum -= ai * ai;
126 asum += ai * (1. - x[e]);
127 n_free = QPMaxHeap_delete(free_heap, n_free, breakpts);
128 if (n_free == 0)
129 {
130 a2sum = 0.;
131 break;
132 }
133 }
134 }
135
136 if (n_bound > 0)
137 {
138 while (breakpts[e = bound_heap[1]] >= lambda)
139 {
140 n_bound = QPMaxHeap_delete(bound_heap, n_bound, breakpts);
141 ai = (a) ? a[e] : 1;
142 a2sum += ai * ai;
143 asum += ai * x[e];
144 t = (x[e] - 1.) / ai;
145 breakpts[e] = t;
146 n_free = QPMaxHeap_add(e, free_heap, breakpts, n_free);
147 if (n_bound == 0)
148 break;
149 }
150 }
151
152 /*------------------------------------------------------------------- */
153 /* get the biggest entry in each heap */
154 /*------------------------------------------------------------------- */
155
156 maxfree = (n_free > 0 ? breakpts[free_heap[1]] : -INFINITY);
157 maxbound = (n_bound > 0 ? breakpts[bound_heap[1]] : -INFINITY);
158 }
159
160 /* This should never happen */
161 ASSERT(false);
162 lambda = 0.;
163 return lambda;
164 }
165
166 } // end namespace Mongoose
167