1 /* $Id: IdiSolve.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #include "CoinPragma.hpp"
7 #include <stdio.h>
8 #include <stdarg.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include "CoinHelperFunctions.hpp"
12 #include "Idiot.hpp"
13 #define FIT
14 #ifdef FIT
15 #define HISTORY 8
16 #else
17 #define HISTORY 7
18 #endif
19 #define NSOLVE HISTORY - 1
solveSmall(int nsolve,double ** aIn,double ** a,double * b)20 static void solveSmall(int nsolve, double **aIn, double **a, double *b)
21 {
22 int i, j;
23 /* copy */
24 for (i = 0; i < nsolve; i++) {
25 for (j = 0; j < nsolve; j++) {
26 a[i][j] = aIn[i][j];
27 }
28 }
29 for (i = 0; i < nsolve; i++) {
30 /* update using all previous */
31 double diagonal;
32 int j;
33 for (j = i; j < nsolve; j++) {
34 int k;
35 double value = a[i][j];
36 for (k = 0; k < i; k++) {
37 value -= a[k][i] * a[k][j];
38 }
39 a[i][j] = value;
40 }
41 diagonal = a[i][i];
42 if (diagonal < 1.0e-20) {
43 diagonal = 0.0;
44 } else {
45 diagonal = 1.0 / sqrt(diagonal);
46 }
47 a[i][i] = diagonal;
48 for (j = i + 1; j < nsolve; j++) {
49 a[i][j] *= diagonal;
50 }
51 }
52 /* update */
53 for (i = 0; i < nsolve; i++) {
54 int j;
55 double value = b[i];
56 for (j = 0; j < i; j++) {
57 value -= b[j] * a[j][i];
58 }
59 value *= a[i][i];
60 b[i] = value;
61 }
62 for (i = nsolve - 1; i >= 0; i--) {
63 int j;
64 double value = b[i];
65 for (j = i + 1; j < nsolve; j++) {
66 value -= b[j] * a[i][j];
67 }
68 value *= a[i][i];
69 b[i] = value;
70 }
71 }
72 IdiotResult
objval(int nrows,int ncols,double * rowsol,double * colsol,double * pi,double *,const double * cost,const double *,const double * rowupper,const double *,const double *,const double * elemnt,const int * row,const CoinBigIndex * columnStart,const int * length,int extraBlock,int * rowExtra,double * solExtra,double * elemExtra,double *,double * costExtra,double weight)73 Idiot::objval(int nrows, int ncols, double *rowsol, double *colsol,
74 double *pi, double * /*djs*/, const double *cost,
75 const double * /*rowlower*/,
76 const double *rowupper, const double * /*lower*/,
77 const double * /*upper*/, const double *elemnt,
78 const int *row, const CoinBigIndex *columnStart,
79 const int *length, int extraBlock, int *rowExtra,
80 double *solExtra, double *elemExtra, double * /*upperExtra*/,
81 double *costExtra, double weight)
82 {
83 IdiotResult result;
84 double objvalue = 0.0;
85 double sum1 = 0.0, sum2 = 0.0;
86 int i;
87 for (i = 0; i < nrows; i++) {
88 rowsol[i] = -rowupper[i];
89 }
90 for (i = 0; i < ncols; i++) {
91 CoinBigIndex j;
92 double value = colsol[i];
93 if (value) {
94 objvalue += value * cost[i];
95 if (elemnt) {
96 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
97 int irow = row[j];
98 rowsol[irow] += elemnt[j] * value;
99 }
100 } else {
101 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
102 int irow = row[j];
103 rowsol[irow] += value;
104 }
105 }
106 }
107 }
108 /* adjust to make as feasible as possible */
109 /* no */
110 if (extraBlock) {
111 for (i = 0; i < extraBlock; i++) {
112 double element = elemExtra[i];
113 int irow = rowExtra[i];
114 objvalue += solExtra[i] * costExtra[i];
115 rowsol[irow] += solExtra[i] * element;
116 }
117 }
118 for (i = 0; i < nrows; i++) {
119 double value = rowsol[i];
120 sum1 += fabs(value);
121 sum2 += value * value;
122 pi[i] = -2.0 * weight * value;
123 }
124 result.infeas = sum1;
125 result.objval = objvalue;
126 result.weighted = objvalue + weight * sum2;
127 result.dropThis = 0.0;
128 result.sumSquared = sum2;
129 return result;
130 }
131 IdiotResult
IdiSolve(int nrows,int ncols,double * COIN_RESTRICT rowsol,double * COIN_RESTRICT colsol,double * COIN_RESTRICT pi,double * COIN_RESTRICT djs,const double * COIN_RESTRICT origcost,double * COIN_RESTRICT rowlower,double * COIN_RESTRICT rowupper,const double * COIN_RESTRICT lower,const double * COIN_RESTRICT upper,const double * COIN_RESTRICT elemnt,const int * row,const CoinBigIndex * columnStart,const int * length,double * COIN_RESTRICT lambda,int maxIts,double mu,double drop,double maxmin,double offset,int strategy,double djTol,double djExit,double djFlag,CoinThreadRandom * randomNumberGenerator)132 Idiot::IdiSolve(
133 int nrows, int ncols, double *COIN_RESTRICT rowsol, double *COIN_RESTRICT colsol,
134 double *COIN_RESTRICT pi, double *COIN_RESTRICT djs, const double *COIN_RESTRICT origcost, double *COIN_RESTRICT rowlower,
135 double *COIN_RESTRICT rowupper, const double *COIN_RESTRICT lower,
136 const double *COIN_RESTRICT upper, const double *COIN_RESTRICT elemnt,
137 const int *row, const CoinBigIndex *columnStart,
138 const int *length, double *COIN_RESTRICT lambda,
139 int maxIts, double mu, double drop,
140 double maxmin, double offset,
141 int strategy, double djTol, double djExit, double djFlag,
142 CoinThreadRandom *randomNumberGenerator)
143 {
144 IdiotResult result;
145 int i, k, iter;
146 CoinBigIndex j;
147 double value = 0.0, objvalue = 0.0, weightedObj = 0.0;
148 double tolerance = 1.0e-8;
149 double *history[HISTORY + 1];
150 int ncolx;
151 int nChange;
152 int extraBlock = 0;
153 int *rowExtra = NULL;
154 double *COIN_RESTRICT solExtra = NULL;
155 double *COIN_RESTRICT elemExtra = NULL;
156 double *COIN_RESTRICT upperExtra = NULL;
157 double *COIN_RESTRICT costExtra = NULL;
158 double *COIN_RESTRICT useCostExtra = NULL;
159 double *COIN_RESTRICT saveExtra = NULL;
160 double *COIN_RESTRICT cost = NULL;
161 double saveValue = 1.0e30;
162 double saveOffset = offset;
163 double useOffset = offset;
164 /*#define NULLVECTOR*/
165 #ifndef NULLVECTOR
166 int nsolve = NSOLVE;
167 #else
168 int nsolve = NSOLVE + 1; /* allow for null vector */
169 #endif
170 int nflagged;
171 double *COIN_RESTRICT thetaX;
172 double *COIN_RESTRICT djX;
173 double *COIN_RESTRICT bX;
174 double *COIN_RESTRICT vX;
175 double **aX;
176 double **aworkX;
177 double **allsum;
178 double *COIN_RESTRICT saveSol = 0;
179 const double *COIN_RESTRICT useCost = cost;
180 double bestSol = 1.0e60;
181 double weight = 0.5 / mu;
182 char *statusSave = new char[2 * ncols];
183 char *statusWork = statusSave + ncols;
184 #define DJTEST 5
185 double djSave[DJTEST];
186 double largestDj = 0.0;
187 double smallestDj = 1.0e60;
188 double maxDj = 0.0;
189 int doFull = 0;
190 #define SAVEHISTORY 10
191 #define EVERY (2 * SAVEHISTORY)
192 #define AFTER SAVEHISTORY *(HISTORY + 1)
193 #define DROP 5
194 double after = AFTER;
195 double obj[DROP];
196 double kbad = 0, kgood = 0;
197 if (strategy & 128)
198 after = 999999; /* no acceleration at all */
199 for (i = 0; i < DROP; i++) {
200 obj[i] = 1.0e70;
201 }
202 //#define FOUR_GOES 2
203 #ifdef FOUR_GOES
204 double *COIN_RESTRICT pi2 = new double[3 * nrows];
205 double *COIN_RESTRICT rowsol2 = new double[3 * nrows];
206 double *COIN_RESTRICT piX[4];
207 double *COIN_RESTRICT rowsolX[4];
208 int startsX[2][5];
209 int nChangeX[4];
210 double maxDjX[4];
211 double objvalueX[4];
212 int nflaggedX[4];
213 piX[0] = pi;
214 piX[1] = pi2;
215 piX[2] = pi2 + nrows;
216 piX[3] = piX[2] + nrows;
217 rowsolX[0] = rowsol;
218 rowsolX[1] = rowsol2;
219 rowsolX[2] = rowsol2 + nrows;
220 rowsolX[3] = rowsolX[2] + nrows;
221 #endif
222 allsum = new double *[nsolve];
223 aX = new double *[nsolve];
224 aworkX = new double *[nsolve];
225 thetaX = new double[nsolve];
226 vX = new double[nsolve];
227 bX = new double[nsolve];
228 djX = new double[nsolve];
229 allsum[0] = pi;
230 for (i = 0; i < nsolve; i++) {
231 if (i)
232 allsum[i] = new double[nrows];
233 aX[i] = new double[nsolve];
234 aworkX[i] = new double[nsolve];
235 }
236 /* check = rows */
237 for (i = 0; i < nrows; i++) {
238 if (rowupper[i] - rowlower[i] > tolerance) {
239 extraBlock++;
240 }
241 }
242 cost = new double[ncols];
243 memset(rowsol, 0, nrows * sizeof(double));
244 for (i = 0; i < ncols; i++) {
245 CoinBigIndex j;
246 double value = origcost[i] * maxmin;
247 double value2 = colsol[i];
248 if (elemnt) {
249 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
250 int irow = row[j];
251 value += elemnt[j] * lambda[irow];
252 rowsol[irow] += elemnt[j] * value2;
253 }
254 } else {
255 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
256 int irow = row[j];
257 value += lambda[irow];
258 rowsol[irow] += value2;
259 }
260 }
261 cost[i] = value;
262 }
263 if (extraBlock) {
264 rowExtra = new int[extraBlock];
265 solExtra = new double[extraBlock];
266 elemExtra = new double[extraBlock];
267 upperExtra = new double[extraBlock];
268 costExtra = new double[extraBlock];
269 saveExtra = new double[extraBlock];
270 extraBlock = 0;
271 int nbad = 0;
272 for (i = 0; i < nrows; i++) {
273 if (rowupper[i] - rowlower[i] > tolerance) {
274 double smaller, difference;
275 double value;
276 saveExtra[extraBlock] = rowupper[i];
277 if (fabs(rowupper[i]) > fabs(rowlower[i])) {
278 smaller = rowlower[i];
279 value = -1.0;
280 } else {
281 smaller = rowupper[i];
282 value = 1.0;
283 }
284 if (fabs(smaller) > 1.0e10) {
285 if (!nbad)
286 COIN_DETAIL_PRINT(printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
287 i, smaller));
288 nbad++;
289 if (rowupper[i] < 0.0 || rowlower[i] > 0.0)
290 abort();
291 if (fabs(rowupper[i]) > fabs(rowlower[i])) {
292 rowlower[i] = -0.9e10;
293 smaller = rowlower[i];
294 value = -1.0;
295 } else {
296 rowupper[i] = 0.9e10;
297 saveExtra[extraBlock] = rowupper[i];
298 smaller = rowupper[i];
299 value = 1.0;
300 }
301 }
302 difference = rowupper[i] - rowlower[i];
303 difference = CoinMin(difference, 1.0e31);
304 rowupper[i] = smaller;
305 elemExtra[extraBlock] = value;
306 solExtra[extraBlock] = (rowupper[i] - rowsol[i]) / value;
307 if (solExtra[extraBlock] < 0.0)
308 solExtra[extraBlock] = 0.0;
309 if (solExtra[extraBlock] > difference)
310 solExtra[extraBlock] = difference;
311 costExtra[extraBlock] = lambda[i] * value;
312 upperExtra[extraBlock] = difference;
313 rowsol[i] += value * solExtra[extraBlock];
314 rowExtra[extraBlock++] = i;
315 }
316 }
317 if (nbad)
318 COIN_DETAIL_PRINT(printf("%d bad values - results may be wrong\n", nbad));
319 }
320 for (i = 0; i < nrows; i++) {
321 offset += lambda[i] * rowsol[i];
322 }
323 if ((strategy & 256) != 0) {
324 /* save best solution */
325 saveSol = new double[ncols];
326 CoinMemcpyN(colsol, ncols, saveSol);
327 if (extraBlock) {
328 useCostExtra = new double[extraBlock];
329 memset(useCostExtra, 0, extraBlock * sizeof(double));
330 }
331 useCost = origcost;
332 useOffset = saveOffset;
333 } else {
334 useCostExtra = costExtra;
335 useCost = cost;
336 useOffset = offset;
337 }
338 ncolx = ncols + extraBlock;
339 for (i = 0; i < HISTORY + 1; i++) {
340 history[i] = new double[ncolx];
341 }
342 for (i = 0; i < DJTEST; i++) {
343 djSave[i] = 1.0e30;
344 }
345 #ifndef OSI_IDIOT
346 int numberColumns = model_->numberColumns();
347 for (int i = 0; i < numberColumns; i++) {
348 if (model_->getColumnStatus(i) != ClpSimplex::isFixed)
349 statusSave[i] = 0;
350 else
351 statusSave[i] = 2;
352 }
353 memset(statusSave + numberColumns, 0, ncols - numberColumns);
354 if ((strategy_ & 131072) == 0) {
355 for (int i = 0; i < numberColumns; i++) {
356 if (model_->getColumnStatus(i) == ClpSimplex::isFixed) {
357 assert(colsol[i] < lower[i] + tolerance || colsol[i] > upper[i] - tolerance);
358 }
359 }
360 }
361 #else
362 for (i = 0; i < ncols; i++) {
363 if (upper[i] - lower[i]) {
364 statusSave[i] = 0;
365 } else {
366 statusSave[i] = 1;
367 }
368 }
369 #endif
370 // for two pass method
371 int start[2];
372 int stop[2];
373 int direction = -1;
374 start[0] = 0;
375 stop[0] = ncols;
376 start[1] = 0;
377 stop[1] = 0;
378 iter = 0;
379 for (; iter < maxIts; iter++) {
380 double sum1 = 0.0, sum2 = 0.0;
381 double lastObj = 1.0e70;
382 int good = 0, doScale = 0;
383 if (strategy & 16) {
384 int ii = iter / EVERY + 1;
385 ii = ii * EVERY;
386 if (iter > ii - HISTORY * 2 && (iter & 1) == 0) {
387 double *COIN_RESTRICT x = history[HISTORY - 1];
388 for (i = HISTORY - 1; i > 0; i--) {
389 history[i] = history[i - 1];
390 }
391 history[0] = x;
392 CoinMemcpyN(colsol, ncols, history[0]);
393 CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
394 }
395 }
396 if ((iter % SAVEHISTORY) == 0 || doFull) {
397 if ((strategy & 16) == 0) {
398 double *COIN_RESTRICT x = history[HISTORY - 1];
399 for (i = HISTORY - 1; i > 0; i--) {
400 history[i] = history[i - 1];
401 }
402 history[0] = x;
403 CoinMemcpyN(colsol, ncols, history[0]);
404 CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
405 }
406 }
407 /* start full try */
408 if ((iter % EVERY) == 0 || doFull) {
409 // for next pass
410 direction = -direction;
411 // randomize.
412 // The cast is to avoid gcc compiler warning
413 int kcol = static_cast< int >(ncols * randomNumberGenerator->randomDouble());
414 if (kcol == ncols)
415 kcol = ncols - 1;
416 if (direction > 0) {
417 start[0] = kcol;
418 stop[0] = ncols;
419 start[1] = 0;
420 stop[1] = kcol;
421 #ifdef FOUR_GOES
422 for (int itry = 0; itry < 2; itry++) {
423 int chunk = (stop[itry] - start[itry] + FOUR_GOES - 1) / FOUR_GOES;
424 startsX[itry][0] = start[itry];
425 for (int i = 1; i < 5; i++)
426 startsX[itry][i] = CoinMin(stop[itry], startsX[itry][i - 1] + chunk);
427 }
428 #endif
429 } else {
430 start[0] = kcol;
431 stop[0] = -1;
432 start[1] = ncols - 1;
433 stop[1] = kcol;
434 #ifdef FOUR_GOES
435 for (int itry = 0; itry < 2; itry++) {
436 int chunk = (start[itry] - stop[itry] + FOUR_GOES - 1) / FOUR_GOES;
437 startsX[itry][0] = start[itry];
438 for (int i = 1; i < 5; i++)
439 startsX[itry][i] = CoinMax(stop[itry], startsX[itry][i - 1] - chunk);
440 }
441 #endif
442 }
443 int itry = 0;
444 /*if ((strategy&16)==0) {
445 double * COIN_RESTRICT x=history[HISTORY-1];
446 for (i=HISTORY-1;i>0;i--) {
447 history[i]=history[i-1];
448 }
449 history[0]=x;
450 CoinMemcpyN(colsol,ncols,history[0]);
451 CoinMemcpyN(solExtra,extraBlock,history[0]+ncols);
452 }*/
453 while (!good) {
454 itry++;
455 #define MAXTRY 5
456 if (iter > after && doScale < 2 && itry < MAXTRY) {
457 /* now full one */
458 for (i = 0; i < nrows; i++) {
459 rowsol[i] = -rowupper[i];
460 }
461 sum2 = 0.0;
462 objvalue = 0.0;
463 memset(pi, 0, nrows * sizeof(double));
464 {
465 double *COIN_RESTRICT theta = thetaX;
466 double *COIN_RESTRICT dj = djX;
467 double *COIN_RESTRICT b = bX;
468 double **a = aX;
469 double **awork = aworkX;
470 double *COIN_RESTRICT v = vX;
471 double c;
472 #ifdef FIT
473 int ntot = 0, nsign = 0, ngood = 0, mgood[4] = { 0, 0, 0, 0 };
474 double diff1, diff2, val0, val1, val2, newValue;
475 CoinMemcpyN(colsol, ncols, history[HISTORY - 1]);
476 CoinMemcpyN(solExtra, extraBlock, history[HISTORY - 1] + ncols);
477 #endif
478 dj[0] = 0.0;
479 for (i = 1; i < nsolve; i++) {
480 dj[i] = 0.0;
481 memset(allsum[i], 0, nrows * sizeof(double));
482 }
483 for (i = 0; i < ncols; i++) {
484 double value2 = colsol[i];
485 if (value2 > lower[i] + tolerance) {
486 if (value2 < (upper[i] - tolerance)) {
487 int k;
488 objvalue += value2 * cost[i];
489 #ifdef FIT
490 ntot++;
491 val0 = history[0][i];
492 val1 = history[1][i];
493 val2 = history[2][i];
494 diff1 = val0 - val1;
495 diff2 = val1 - val2;
496 if (diff1 * diff2 >= 0.0) {
497 nsign++;
498 if (fabs(diff1) < fabs(diff2)) {
499 int ii = static_cast< int >(fabs(4.0 * diff1 / diff2));
500 if (ii == 4)
501 ii = 3;
502 mgood[ii]++;
503 ngood++;
504 }
505 if (fabs(diff1) < 0.75 * fabs(diff2)) {
506 newValue = val1 + (diff1 * diff2) / (diff2 - diff1);
507 } else {
508 newValue = val1 + 4.0 * diff1;
509 }
510 } else {
511 newValue = 0.333333333 * (val0 + val1 + val2);
512 }
513 if (newValue > upper[i] - tolerance) {
514 newValue = upper[i];
515 } else if (newValue < lower[i] + tolerance) {
516 newValue = lower[i];
517 }
518 history[HISTORY - 1][i] = newValue;
519 #endif
520 for (k = 0; k < HISTORY - 1; k++) {
521 value = history[k][i] - history[k + 1][i];
522 dj[k] += value * cost[i];
523 v[k] = value;
524 }
525 if (elemnt) {
526 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
527 int irow = row[j];
528 for (k = 0; k < HISTORY - 1; k++) {
529 allsum[k][irow] += elemnt[j] * v[k];
530 }
531 rowsol[irow] += elemnt[j] * value2;
532 }
533 } else {
534 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
535 int irow = row[j];
536 for (k = 0; k < HISTORY - 1; k++) {
537 allsum[k][irow] += v[k];
538 }
539 rowsol[irow] += value2;
540 }
541 }
542 } else {
543 /* at ub */
544 colsol[i] = upper[i];
545 value2 = colsol[i];
546 objvalue += value2 * cost[i];
547 if (elemnt) {
548 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
549 int irow = row[j];
550 rowsol[irow] += elemnt[j] * value2;
551 }
552 } else {
553 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
554 int irow = row[j];
555 rowsol[irow] += value2;
556 }
557 }
558 }
559 } else {
560 /* at lb */
561 if (value2) {
562 objvalue += value2 * cost[i];
563 if (elemnt) {
564 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
565 int irow = row[j];
566 rowsol[irow] += elemnt[j] * value2;
567 }
568 } else {
569 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
570 int irow = row[j];
571 rowsol[irow] += value2;
572 }
573 }
574 }
575 }
576 }
577 #ifdef FIT
578 /*printf("total %d, same sign %d, good %d %d %d %d %d\n",
579 ntot,nsign,ngood,mgood[0],mgood[1],mgood[2],mgood[3]);*/
580 #endif
581 if (extraBlock) {
582 for (i = 0; i < extraBlock; i++) {
583 double element = elemExtra[i];
584 int irow = rowExtra[i];
585 objvalue += solExtra[i] * costExtra[i];
586 if (solExtra[i] > tolerance
587 && solExtra[i] < (upperExtra[i] - tolerance)) {
588 double value2 = solExtra[i];
589 int k;
590 for (k = 0; k < HISTORY - 1; k++) {
591 value = history[k][i + ncols] - history[k + 1][i + ncols];
592 dj[k] += value * costExtra[i];
593 allsum[k][irow] += element * value;
594 }
595 rowsol[irow] += element * value2;
596 } else {
597 double value2 = solExtra[i];
598 double element = elemExtra[i];
599 int irow = rowExtra[i];
600 rowsol[irow] += element * value2;
601 }
602 }
603 }
604 #ifdef NULLVECTOR
605 if ((strategy & 64)) {
606 double djVal = dj[0];
607 for (i = 0; i < ncols - nrows; i++) {
608 double value2 = colsol[i];
609 if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
610 value = history[0][i] - history[1][i];
611 } else {
612 value = 0.0;
613 }
614 history[HISTORY][i] = value;
615 }
616 for (; i < ncols; i++) {
617 double value2 = colsol[i];
618 double delta;
619 int irow = i - (ncols - nrows);
620 double oldSum = allsum[0][irow];
621 ;
622 if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
623 delta = history[0][i] - history[1][i];
624 } else {
625 delta = 0.0;
626 }
627 djVal -= delta * cost[i];
628 oldSum -= delta;
629 delta = -oldSum;
630 djVal += delta * cost[i];
631 history[HISTORY][i] = delta;
632 }
633 dj[HISTORY - 1] = djVal;
634 djVal = 0.0;
635 for (i = 0; i < ncols; i++) {
636 double value2 = colsol[i];
637 if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance || i >= ncols - nrows) {
638 int k;
639 value = history[HISTORY][i];
640 djVal += value * cost[i];
641 for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
642 int irow = row[j];
643 allsum[nsolve - 1][irow] += value;
644 }
645 }
646 }
647 printf("djs %g %g\n", dj[HISTORY - 1], djVal);
648 }
649 #endif
650 for (i = 0; i < nsolve; i++) {
651 int j;
652 b[i] = 0.0;
653 for (j = 0; j < nsolve; j++) {
654 a[i][j] = 0.0;
655 }
656 }
657 c = 0.0;
658 for (i = 0; i < nrows; i++) {
659 double value = rowsol[i];
660 for (k = 0; k < nsolve; k++) {
661 v[k] = allsum[k][i];
662 b[k] += v[k] * value;
663 }
664 c += value * value;
665 for (k = 0; k < nsolve; k++) {
666 for (j = k; j < nsolve; j++) {
667 a[k][j] += v[k] * v[j];
668 }
669 }
670 }
671 sum2 = c;
672 if (itry == 1) {
673 lastObj = objvalue + weight * sum2;
674 }
675 for (k = 0; k < nsolve; k++) {
676 b[k] = -(weight * b[k] + 0.5 * dj[k]);
677 for (j = k; j < nsolve; j++) {
678 a[k][j] *= weight;
679 a[j][k] = a[k][j];
680 }
681 }
682 c *= weight;
683 for (k = 0; k < nsolve; k++) {
684 theta[k] = b[k];
685 }
686 solveSmall(nsolve, a, awork, theta);
687 if ((strategy & 64) != 0) {
688 value = 10.0;
689 for (k = 0; k < nsolve; k++) {
690 value = CoinMax(value, fabs(theta[k]));
691 }
692 if (value > 10.0 && ((logLevel_ & 4) != 0)) {
693 printf("theta %g %g %g\n", theta[0], theta[1], theta[2]);
694 }
695 value = 10.0 / value;
696 for (k = 0; k < nsolve; k++) {
697 theta[k] *= value;
698 }
699 }
700 for (i = 0; i < ncolx; i++) {
701 double valueh = 0.0;
702 for (k = 0; k < HISTORY - 1; k++) {
703 value = history[k][i] - history[k + 1][i];
704 valueh += value * theta[k];
705 }
706 #ifdef NULLVECTOR
707 value = history[HISTORY][i];
708 valueh += value * theta[HISTORY - 1];
709 #endif
710 history[HISTORY][i] = valueh;
711 }
712 }
713 #ifdef NULLVECTOR
714 if ((strategy & 64)) {
715 for (i = 0; i < ncols - nrows; i++) {
716 if (colsol[i] <= lower[i] + tolerance
717 || colsol[i] >= (upper[i] - tolerance)) {
718 history[HISTORY][i] = 0.0;
719 ;
720 }
721 }
722 tolerance = -tolerance; /* switch off test */
723 }
724 #endif
725 if (!doScale) {
726 for (i = 0; i < ncols; i++) {
727 if (colsol[i] > lower[i] + tolerance
728 && colsol[i] < (upper[i] - tolerance)) {
729 value = history[HISTORY][i];
730 colsol[i] += value;
731 if (colsol[i] < lower[i] + tolerance) {
732 colsol[i] = lower[i];
733 } else if (colsol[i] > upper[i] - tolerance) {
734 colsol[i] = upper[i];
735 }
736 }
737 }
738 if (extraBlock) {
739 for (i = 0; i < extraBlock; i++) {
740 if (solExtra[i] > tolerance
741 && solExtra[i] < (upperExtra[i] - tolerance)) {
742 value = history[HISTORY][i + ncols];
743 solExtra[i] += value;
744 if (solExtra[i] < 0.0) {
745 solExtra[i] = 0.0;
746 } else if (solExtra[i] > upperExtra[i]) {
747 solExtra[i] = upperExtra[i];
748 }
749 }
750 }
751 }
752 } else {
753 double theta = 1.0;
754 double saveTheta = theta;
755 for (i = 0; i < ncols; i++) {
756 if (colsol[i] > lower[i] + tolerance
757 && colsol[i] < (upper[i] - tolerance)) {
758 value = history[HISTORY][i];
759 if (value > 0) {
760 if (theta * value + colsol[i] > upper[i]) {
761 theta = (upper[i] - colsol[i]) / value;
762 }
763 } else if (value < 0) {
764 if (colsol[i] + theta * value < lower[i]) {
765 theta = (lower[i] - colsol[i]) / value;
766 }
767 }
768 }
769 }
770 if (extraBlock) {
771 for (i = 0; i < extraBlock; i++) {
772 if (solExtra[i] > tolerance
773 && solExtra[i] < (upperExtra[i] - tolerance)) {
774 value = history[HISTORY][i + ncols];
775 if (value > 0) {
776 if (theta * value + solExtra[i] > upperExtra[i]) {
777 theta = (upperExtra[i] - solExtra[i]) / value;
778 }
779 } else if (value < 0) {
780 if (solExtra[i] + theta * value < 0.0) {
781 theta = -solExtra[i] / value;
782 }
783 }
784 }
785 }
786 }
787 if ((iter % 100 == 0) && (logLevel_ & 8) != 0) {
788 if (theta < saveTheta) {
789 printf(" - modified theta %g\n", theta);
790 }
791 }
792 for (i = 0; i < ncols; i++) {
793 if (colsol[i] > lower[i] + tolerance
794 && colsol[i] < (upper[i] - tolerance)) {
795 value = history[HISTORY][i];
796 colsol[i] += value * theta;
797 }
798 }
799 if (extraBlock) {
800 for (i = 0; i < extraBlock; i++) {
801 if (solExtra[i] > tolerance
802 && solExtra[i] < (upperExtra[i] - tolerance)) {
803 value = history[HISTORY][i + ncols];
804 solExtra[i] += value * theta;
805 }
806 }
807 }
808 }
809 #ifdef NULLVECTOR
810 tolerance = fabs(tolerance); /* switch back on */
811 #endif
812 if ((iter % 100) == 0 && (logLevel_ & 8) != 0) {
813 printf("\n");
814 }
815 }
816 good = 1;
817 result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
818 rowlower, rowupper, lower, upper,
819 elemnt, row, columnStart, length, extraBlock, rowExtra,
820 solExtra, elemExtra, upperExtra, useCostExtra,
821 weight);
822 weightedObj = result.weighted;
823 if (!iter)
824 saveValue = weightedObj;
825 objvalue = result.objval;
826 sum1 = result.infeas;
827 if (saveSol) {
828 if (result.weighted < bestSol) {
829 COIN_DETAIL_PRINT(printf("%d %g better than %g\n", iter,
830 result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
831 bestSol = result.weighted;
832 CoinMemcpyN(colsol, ncols, saveSol);
833 }
834 }
835 #ifdef FITz
836 if (iter > after) {
837 IdiotResult result2;
838 double ww, oo, ss;
839 if (extraBlock)
840 abort();
841 result2 = objval(nrows, ncols, row2, sol2, pi2, djs, cost,
842 rowlower, rowupper, lower, upper,
843 elemnt, row, columnStart, extraBlock, rowExtra,
844 solExtra, elemExtra, upperExtra, costExtra,
845 weight);
846 ww = result2.weighted;
847 oo = result2.objval;
848 ss = result2.infeas;
849 printf("wobj %g obj %g inf %g last %g\n", ww, oo, ss, lastObj);
850 if (ww < weightedObj && ww < lastObj) {
851 printf(" taken");
852 ntaken++;
853 saving += weightedObj - ww;
854 weightedObj = ww;
855 objvalue = oo;
856 sum1 = ss;
857 CoinMemcpyN(row2, nrows, rowsol);
858 CoinMemcpyN(pi2, nrows, pi);
859 CoinMemcpyN(sol2, ncols, colsol);
860 result = objval(nrows, ncols, rowsol, colsol, pi, djs, cost,
861 rowlower, rowupper, lower, upper,
862 elemnt, row, columnStart, extraBlock, rowExtra,
863 solExtra, elemExtra, upperExtra, costExtra,
864 weight);
865 weightedObj = result.weighted;
866 objvalue = result.objval;
867 sum1 = result.infeas;
868 if (ww < weightedObj)
869 abort();
870 } else {
871 printf(" not taken");
872 nottaken++;
873 }
874 }
875 #endif
876 /*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/
877 if (weightedObj > lastObj + 1.0e-4 && itry < MAXTRY) {
878 if ((logLevel_ & 16) != 0 && doScale) {
879 printf("Weighted objective from %g to %g **** bad move\n",
880 lastObj, weightedObj);
881 }
882 if (doScale) {
883 good = 1;
884 }
885 if ((strategy & 3) == 1) {
886 good = 0;
887 if (weightedObj > lastObj + djExit) {
888 if ((logLevel_ & 16) != 0) {
889 printf("Weighted objective from %g to %g ?\n", lastObj, weightedObj);
890 }
891 CoinMemcpyN(history[0], ncols, colsol);
892 CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
893 good = 1;
894 }
895 } else if ((strategy & 3) == 2) {
896 if (weightedObj > lastObj + 0.1 * maxDj) {
897 CoinMemcpyN(history[0], ncols, colsol);
898 CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
899 doScale++;
900 good = 0;
901 }
902 } else if ((strategy & 3) == 3) {
903 if (weightedObj > lastObj + 0.001 * maxDj) {
904 /*doScale++;*/
905 good = 0;
906 }
907 }
908 }
909 }
910 if ((iter % checkFrequency_) == 0) {
911 double best = weightedObj;
912 double test = obj[0];
913 for (i = 1; i < DROP; i++) {
914 obj[i - 1] = obj[i];
915 if (best > obj[i])
916 best = obj[i];
917 }
918 obj[DROP - 1] = best;
919 if (test - best < drop && (strategy & 8) == 0) {
920 if ((logLevel_ & 8) != 0) {
921 printf("Exiting as drop in %d its is %g after %d iterations\n",
922 DROP * checkFrequency_, test - best, iter);
923 }
924 goto RETURN;
925 }
926 }
927 if ((iter % logFreq_) == 0) {
928 double piSum = 0.0;
929 for (i = 0; i < nrows; i++) {
930 piSum += (rowsol[i] + rowupper[i]) * pi[i];
931 }
932 if ((logLevel_ & 2) != 0) {
933 printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n",
934 iter, sum1, objvalue * maxmin - useOffset, weightedObj - useOffset,
935 piSum * maxmin - useOffset, maxDj);
936 }
937 }
938 CoinMemcpyN(statusSave, ncols, statusWork);
939 nflagged = 0;
940 }
941 nChange = 0;
942 doFull = 0;
943 maxDj = 0.0;
944 // go through forwards or backwards and starting at odd places
945 #ifdef FOUR_GOES
946 for (int i = 1; i < FOUR_GOES; i++) {
947 cilk_spawn memcpy(piX[i], pi, nrows * sizeof(double));
948 cilk_spawn memcpy(rowsolX[i], rowsol, nrows * sizeof(double));
949 }
950 for (int i = 0; i < FOUR_GOES; i++) {
951 nChangeX[i] = 0;
952 maxDjX[i] = 0.0;
953 objvalueX[i] = 0.0;
954 nflaggedX[i] = 0;
955 }
956 cilk_sync;
957 #endif
958 //printf("PASS\n");
959 #ifdef FOUR_GOES
960 cilk_for(int iPar = 0; iPar < FOUR_GOES; iPar++)
961 {
962 double *COIN_RESTRICT pi = piX[iPar];
963 double *COIN_RESTRICT rowsol = rowsolX[iPar];
964 for (int itry = 0; itry < 2; itry++) {
965 int istop;
966 int istart;
967 #if 0
968 int chunk = (start[itry]+stop[itry]+FOUR_GOES-1)/FOUR_GOES;
969 if (iPar == 0) {
970 istart=start[itry];
971 istop=(start[itry]+stop[itry])>>1;
972 } else {
973 istart=(start[itry]+stop[itry])>>1;
974 istop = stop[itry];
975 }
976 #endif
977 #if 0
978 printf("istart %d istop %d direction %d array %d %d new %d %d\n",
979 istart,istop,direction,start[itry],stop[itry],
980 startsX[itry][iPar],startsX[itry][iPar+1]);
981 #endif
982 istart = startsX[itry][iPar];
983 istop = startsX[itry][iPar + 1];
984 #else
985 for (int itry = 0; itry < 2; itry++) {
986 int istart = start[itry];
987 int istop = stop[itry];
988 #endif
989 for (int icol = istart; icol != istop; icol += direction) {
990 if (!statusWork[icol]) {
991 CoinBigIndex j;
992 double value = colsol[icol];
993 double djval = cost[icol];
994 double djval2, value2;
995 double theta, a, b, c;
996 if (elemnt) {
997 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
998 int irow = row[j];
999 djval -= elemnt[j] * pi[irow];
1000 }
1001 } else {
1002 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
1003 int irow = row[j];
1004 djval -= pi[irow];
1005 }
1006 }
1007 /*printf("xx iter %d seq %d djval %g value %g\n",
1008 iter,i,djval,value);*/
1009 if (djval > 1.0e-5) {
1010 value2 = (lower[icol] - value);
1011 } else {
1012 value2 = (upper[icol] - value);
1013 }
1014 djval2 = djval * value2;
1015 djval = fabs(djval);
1016 if (djval > djTol) {
1017 if (djval2 < -1.0e-4) {
1018 #ifndef FOUR_GOES
1019 nChange++;
1020 if (djval > maxDj)
1021 maxDj = djval;
1022 #else
1023 nChangeX[iPar]++;
1024 if (djval > maxDjX[iPar])
1025 maxDjX[iPar] = djval;
1026 #endif
1027 /*if (djval>3.55e6) {
1028 printf("big\n");
1029 }*/
1030 a = 0.0;
1031 b = 0.0;
1032 c = 0.0;
1033 djval2 = cost[icol];
1034 if (elemnt) {
1035 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
1036 int irow = row[j];
1037 double value = rowsol[irow];
1038 c += value * value;
1039 a += elemnt[j] * elemnt[j];
1040 b += value * elemnt[j];
1041 }
1042 } else {
1043 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
1044 int irow = row[j];
1045 double value = rowsol[irow];
1046 c += value * value;
1047 a += 1.0;
1048 b += value;
1049 }
1050 }
1051 a *= weight;
1052 b = b * weight + 0.5 * djval2;
1053 c *= weight;
1054 /* solve */
1055 theta = -b / a;
1056 #ifndef FOUR_GOES
1057 if ((strategy & 4) != 0) {
1058 double valuep, thetap;
1059 value2 = a * theta * theta + 2.0 * b * theta;
1060 thetap = 2.0 * theta;
1061 valuep = a * thetap * thetap + 2.0 * b * thetap;
1062 if (valuep < value2 + djTol) {
1063 theta = thetap;
1064 kgood++;
1065 } else {
1066 kbad++;
1067 }
1068 }
1069 #endif
1070 if (theta > 0.0) {
1071 if (theta < upper[icol] - colsol[icol]) {
1072 value2 = theta;
1073 } else {
1074 value2 = upper[icol] - colsol[icol];
1075 }
1076 } else {
1077 if (theta > lower[icol] - colsol[icol]) {
1078 value2 = theta;
1079 } else {
1080 value2 = lower[icol] - colsol[icol];
1081 }
1082 }
1083 colsol[icol] += value2;
1084 #ifndef FOUR_GOES
1085 objvalue += cost[icol] * value2;
1086 #else
1087 objvalueX[iPar] += cost[icol] * value2;
1088 #endif
1089 if (elemnt) {
1090 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
1091 int irow = row[j];
1092 double value;
1093 rowsol[irow] += elemnt[j] * value2;
1094 value = rowsol[irow];
1095 pi[irow] = -2.0 * weight * value;
1096 }
1097 } else {
1098 for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
1099 int irow = row[j];
1100 double value;
1101 rowsol[irow] += value2;
1102 value = rowsol[irow];
1103 pi[irow] = -2.0 * weight * value;
1104 }
1105 }
1106 } else {
1107 /* dj but at bound */
1108 if (djval > djFlag) {
1109 statusWork[icol] = 1;
1110 #ifndef FOUR_GOES
1111 nflagged++;
1112 #else
1113 nflaggedX[iPar]++;
1114 #endif
1115 }
1116 }
1117 }
1118 }
1119 }
1120 #ifdef FOUR_GOES
1121 }
1122 #endif
1123 }
1124 #ifdef FOUR_GOES
1125 for (int i = 0; i < FOUR_GOES; i++) {
1126 nChange += nChangeX[i];
1127 maxDj = CoinMax(maxDj, maxDjX[i]);
1128 objvalue += objvalueX[i];
1129 nflagged += nflaggedX[i];
1130 }
1131 cilk_for(int i = 0; i < nrows; i++)
1132 {
1133 #if FOUR_GOES == 2
1134 rowsol[i] = 0.5 * (rowsolX[0][i] + rowsolX[1][i]);
1135 pi[i] = 0.5 * (piX[0][i] + piX[1][i]);
1136 #elif FOUR_GOES == 3
1137 pi[i] = 0.33333333333333 * (piX[0][i] + piX[1][i] + piX[2][i]);
1138 rowsol[i] = 0.3333333333333 * (rowsolX[0][i] + rowsolX[1][i] + rowsolX[2][i]);
1139 #else
1140 pi[i] = 0.25 * (piX[0][i] + piX[1][i] + piX[2][i] + piX[3][i]);
1141 rowsol[i] = 0.25 * (rowsolX[0][i] + rowsolX[1][i] + rowsolX[2][i] + rowsolX[3][i]);
1142 #endif
1143 }
1144 #endif
1145 if (extraBlock) {
1146 for (int i = 0; i < extraBlock; i++) {
1147 double value = solExtra[i];
1148 double djval = costExtra[i];
1149 double djval2, value2;
1150 double element = elemExtra[i];
1151 double theta, a, b, c;
1152 int irow = rowExtra[i];
1153 djval -= element * pi[irow];
1154 /*printf("xxx iter %d extra %d djval %g value %g\n",
1155 iter,irow,djval,value);*/
1156 if (djval > 1.0e-5) {
1157 value2 = -value;
1158 } else {
1159 value2 = (upperExtra[i] - value);
1160 }
1161 djval2 = djval * value2;
1162 if (djval2 < -1.0e-4 && fabs(djval) > djTol) {
1163 nChange++;
1164 a = 0.0;
1165 b = 0.0;
1166 c = 0.0;
1167 djval2 = costExtra[i];
1168 value = rowsol[irow];
1169 c += value * value;
1170 a += element * element;
1171 b += element * value;
1172 a *= weight;
1173 b = b * weight + 0.5 * djval2;
1174 c *= weight;
1175 /* solve */
1176 theta = -b / a;
1177 if (theta > 0.0) {
1178 value2 = CoinMin(theta, upperExtra[i] - solExtra[i]);
1179 } else {
1180 value2 = CoinMax(theta, -solExtra[i]);
1181 }
1182 solExtra[i] += value2;
1183 rowsol[irow] += element * value2;
1184 value = rowsol[irow];
1185 pi[irow] = -2.0 * weight * value;
1186 }
1187 }
1188 }
1189 if ((iter % 10) == 2) {
1190 for (int i = DJTEST - 1; i > 0; i--) {
1191 djSave[i] = djSave[i - 1];
1192 }
1193 djSave[0] = maxDj;
1194 largestDj = CoinMax(largestDj, maxDj);
1195 smallestDj = CoinMin(smallestDj, maxDj);
1196 for (int i = DJTEST - 1; i > 0; i--) {
1197 maxDj += djSave[i];
1198 }
1199 maxDj = maxDj / static_cast< double >(DJTEST);
1200 if (maxDj < djExit && iter > 50) {
1201 //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
1202 break;
1203 }
1204 if (nChange < 100) {
1205 djTol *= 0.5;
1206 }
1207 }
1208 }
1209 RETURN:
1210 if (kgood || kbad) {
1211 COIN_DETAIL_PRINT(printf("%g good %g bad\n", kgood, kbad));
1212 }
1213 result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
1214 rowlower, rowupper, lower, upper,
1215 elemnt, row, columnStart, length, extraBlock, rowExtra,
1216 solExtra, elemExtra, upperExtra, useCostExtra,
1217 weight);
1218 result.djAtBeginning = largestDj;
1219 result.djAtEnd = smallestDj;
1220 result.dropThis = saveValue - result.weighted;
1221 if (saveSol) {
1222 if (result.weighted < bestSol) {
1223 bestSol = result.weighted;
1224 CoinMemcpyN(colsol, ncols, saveSol);
1225 } else {
1226 COIN_DETAIL_PRINT(printf("restoring previous - now %g best %g\n",
1227 result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
1228 }
1229 }
1230 if (saveSol) {
1231 if (extraBlock) {
1232 delete[] useCostExtra;
1233 }
1234 CoinMemcpyN(saveSol, ncols, colsol);
1235 delete[] saveSol;
1236 }
1237 for (i = 0; i < nsolve; i++) {
1238 if (i)
1239 delete[] allsum[i];
1240 delete[] aX[i];
1241 delete[] aworkX[i];
1242 }
1243 delete[] thetaX;
1244 delete[] djX;
1245 delete[] bX;
1246 delete[] vX;
1247 delete[] aX;
1248 delete[] aworkX;
1249 delete[] allsum;
1250 delete[] cost;
1251 #ifdef FOUR_GOES
1252 delete[] pi2;
1253 delete[] rowsol2;
1254 #endif
1255 for (i = 0; i < HISTORY + 1; i++) {
1256 delete[] history[i];
1257 }
1258 delete[] statusSave;
1259 /* do original costs objvalue*/
1260 result.objval = 0.0;
1261 for (i = 0; i < ncols; i++) {
1262 result.objval += colsol[i] * origcost[i];
1263 }
1264 if (extraBlock) {
1265 for (i = 0; i < extraBlock; i++) {
1266 int irow = rowExtra[i];
1267 rowupper[irow] = saveExtra[i];
1268 }
1269 delete[] rowExtra;
1270 delete[] solExtra;
1271 delete[] elemExtra;
1272 delete[] upperExtra;
1273 delete[] costExtra;
1274 delete[] saveExtra;
1275 }
1276 result.iteration = iter;
1277 result.objval -= saveOffset;
1278 result.weighted = result.objval + weight * result.sumSquared;
1279 return result;
1280 }
1281
1282 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1283 */
1284