1 /*
2 * Name: move.c
3 * Author: Pietro Belotti
4 * Purpose: move to next iterate and update all sparse data structures
5 *
6 * This code is published under the Eclipse Public License (EPL).
7 * See http://www.eclipse.org/legal/epl-v10.html
8 */
9
10 #include <stdlib.h>
11 #include <math.h>
12 #include <stdio.h>
13
14 #ifdef RTR_MPI
15 #include <mpi.h>
16 #endif
17
18 #include "sparse.h"
19 #include "linopt.h"
20 #include "rtr.h"
21
22 #ifdef RTR_USE_PRAGMAS
23
24 /*
25 * fast vector operations (use #pragmas)
26 */
27
build_dxk_block(int * ip,double * ic,double * dxk,double z)28 __inline void build_dxk_block (int *ip, double *ic, double *dxk, double z) {
29
30 register int i = CONCURRENT_FLOPS;
31
32 #pragma disjoint (*dxk, *ip, z, *ic)
33
34 for (; i>0; --i) dxk [ip [i]] += z * ic [i]; /* update variation of variable *ip */
35 }
36
37
build_dxk(int k,int * ip,double * ic,double * dxk,double z)38 __inline void build_dxk (int k, int *ip, double *ic, double *dxk, double z) {
39
40 register int i = k / CONCURRENT_FLOPS;
41
42 k %= CONCURRENT_FLOPS;
43
44 for (; i>0; --i) {
45
46 build_dxk_block (ip, ic, dxk, z);
47
48 ip += 10;
49 ic += 10;
50 }
51
52 /* there is something fishy here... why not call it with k%...? */
53 }
54
55 #endif
56
57
58 /*
59 * Update points and aux. vectors for a delta vector
60 */
61
update(sparseLP * lp,double * x,double * dx,double * b_Ax,char * sat,double * sum,double maxMod)62 int update (sparseLP *lp, double *x, double *dx, double *b_Ax,
63 char *sat, double *sum, double maxMod) {
64
65 int i, j, dSatd = 0;
66
67 int c = lp->c0;
68 int *vl = lp->vl;
69
70 int **colP = lp->vp;
71 double **colC = lp->vc;
72
73 double *ub = lp->ub;
74 double *lb = lp->lb;
75
76 int *vp;
77 double *vc;
78
79 double z;
80
81 for (i=0; i<c; i++, vl++, colC++, colP++, dx++, x++, ub++, lb++) {
82
83 if ((*dx > EPSILON) || (*dx < EPSILON)) {
84
85 vp = *colP;
86 vc = *colC;
87
88 z = *x;
89
90 if ((*x += maxMod * *dx) > *ub) {*x = *ub;}
91 else if (*x < *lb) {*x = *lb;}
92
93 z -= *x;
94
95 /* for all constraints containing this variable, update b_Ax */
96 /* and respective indicator */
97
98 #ifdef RTR_USE_PRAGMAS
99 #pragma execution_frequency(very_high)
100 #endif
101
102 for (j = *vl; j>0 ; --j) {
103
104 register double *pV = b_Ax + *vp;
105 register char *pS = sat + *vp++;
106 register double delta = z * *vc++;
107
108 if ((*pV) > 0) {
109
110 if ((delta > 0) || ((*pV) > - delta))
111 *sum += delta;
112 else {
113 *sum -= *pV;
114 if (!*pS) { dSatd++; *pS = SATD; } /* constraint is now fulfilled */
115 }
116 }
117 else {
118
119 if (delta > - (*pV)) {
120 *sum += (*pV + delta);
121 if (*pS) { dSatd--; *pS = UNSATD; } /* constraint is now violated */
122 }
123 }
124
125 *pV += delta; /* update b_Ax */
126 }
127
128 *dx = 0;
129 }
130 }
131
132 return dSatd;
133 }
134
135
136 /*
137 * create move direction and move to new x accordingly
138 */
139
move(sparseLP * lp,double * x,double * b_Ax,int * block,char * sat,double temp,double * sum)140 int move (sparseLP *lp, /* LP data */
141 double *x, /* current point */
142 double *b_Ax, /* current violation */
143 int *block, /* array with indices of selected constraints */
144 char *sat, /* indicator function of fulfilled constraints */
145 double temp, /* temperature */
146 double *sum /* total violation (updated here) */
147 ) {
148
149 static double *dx0 = NULL; /* move direction (global) */
150 static double *dxk = NULL; /* move direction (local) */
151
152 static char first = 1; /* is this the first call to the function? */
153
154 int i, c, r;
155
156 double *ub, *lb;
157
158 double z, maxMod;
159
160 /* free static malloc'ed arrays */
161
162 if (!lp) {
163
164 free (dx0);
165 free (dxk);
166 return 0;
167 }
168
169 /* first call, initialize local arrays */
170
171 c = lp->c0;
172 r = lp->rk;
173
174 if (first) {
175
176 first = 0;
177
178 dx0 = (double *) malloc (c * sizeof (double));
179 dxk = (double *) malloc (c * sizeof (double));
180
181 for (i=c; i>0; i--) *dxk++ = 0.0;
182
183 dxk -= c;
184 }
185
186 /* Compute dx based on b_Ax */
187
188 lb = lp -> lb;
189 ub = lp -> ub;
190
191 for (; *block >= 0; block++) {
192
193 int *ip = (lp->ip) [*block];
194 double *ic = (lp->ic) [*block];
195
196 z = b_Ax [*block];
197
198 if (z<0) z = 0;
199 else z = exp (- z / temp);
200
201 #ifdef RTR_USE_PRAGMAS
202 build_dxk ((lp->il) [*block], ip, ic, dxk, z);
203 #else
204 for (i = (lp->il) [*block]; i>0; i--, ip++)
205 dxk [*ip] += z * *ic++; /* update variation of variable *ip */
206 #endif
207 }
208
209 #ifdef RTR_MPI
210 MPI_Allreduce (dxk, dx0, c, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
211 #else
212 for (i=c; i>0; --i) *dx0++ = *dxk++;
213 dxk -= c;
214 dx0 -= c;
215 #endif
216
217 for (i=c; i>0; --i) *dxk++ = 0;
218 dxk -= c;
219
220 /* Project back into bounding box */
221
222 if (lp -> onedim) maxMod = lp -> stretch;
223 else maxMod = lp -> alpha;
224
225 /*
226 {
227 double mind = 1e40, maxd = -1e40, norm=0;
228
229 for (i=c;i>0;i--, dx0++) {
230 if (fabs (*dx0) > maxd) maxd = *dx0;
231 if (fabs (*dx0) < mind) mind = *dx0;
232 norm += *dx0 * *dx0;
233 }
234 norm = sqrt (norm);
235 dx0 -= c;
236 printf ("||dx0|| = %12.8f (%12.8f,%12.8f); ", norm, mind, maxd);
237 }
238 */
239
240 for (i=c; i>0; --i, dx0++, x++, ub++, lb++) {
241
242 if (*dx0 < - EPSILON) { /* move down, check if below lower bound */
243 if (*x <= *lb) *dx0 = 0; /* already on boundary */
244 else if ((z = (*lb - *x) / *dx0) < maxMod) maxMod = z;
245 }
246
247 else if (*dx0 > EPSILON) { /* move up, check if above upper bound */
248 if (*x >= *ub) *dx0 = 0; /* already on boundary */
249 else if ((z = (*ub - *x) / *dx0) < maxMod) maxMod = z;
250 }
251
252 /*
253 * Now maxMod contains a scaling factor for (x + maxMod * dx) to
254 * stay within bounding box
255 */
256 }
257
258 dx0 -= c;
259 x -= c;
260 ub -= c;
261 lb -= c;
262
263 /* If one-dimensional optimization is used, tune the next point so
264 * as to get the most fulfilling point in the segment
265 */
266
267 /* printf ("maxMod = %12.8f ", maxMod); */
268
269
270 if (lp -> onedim) { /* create sparse version of dx0 */
271
272 z = one_opt (lp, dx0, b_Ax, maxMod);
273
274 if (z > 0) maxMod = z;
275 }
276
277 /* printf ("-> %12.8f\n", maxMod); */
278
279 /* {
280 double summ = 0;
281 for (i=0; i<c; i++, dx0++) summ += *dx0 * *dx0;
282 printf ("||dx|| = %.6f\n", log (1e-15+fabs(maxMod * sqrt(summ))) / log (10));
283 dx0 -= c;
284 } */
285
286 /* Update b_Ax */
287
288 return update (lp, x, dx0, b_Ax, sat, sum, maxMod);
289 }
290