source: trunk/Clp/src/AbcMatrix.cpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 10 months ago

formatting

  • Property svn:keywords set to Id
File size: 128.8 KB
Line 
1/* $Id: AbcMatrix.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include <cstdio>
7#include "CoinPragma.hpp"
8#include "CoinIndexedVector.hpp"
9#include "CoinHelperFunctions.hpp"
10#include "CoinAbcHelperFunctions.hpp"
11#include "AbcSimplexFactorization.hpp"
12#include "AbcPrimalColumnDantzig.hpp"
13#include "AbcPrimalColumnSteepest.hpp"
14#include "CoinTime.hpp"
15
16#include "AbcSimplex.hpp"
17#include "AbcSimplexDual.hpp"
18// at end to get min/max!
19#include "AbcMatrix.hpp"
20#include "ClpMessage.hpp"
21#ifdef INTEL_MKL
22#include "mkl_spblas.h"
23#endif
24#if ABC_INSTRUMENT > 1
25extern int abcPricing[20];
26extern int abcPricingDense[20];
27#endif
28//=============================================================================
29
30//#############################################################################
31// Constructors / Destructor / Assignment
32//#############################################################################
33
34//-------------------------------------------------------------------
35// Default Constructor
36//-------------------------------------------------------------------
37AbcMatrix::AbcMatrix()
38  : matrix_(NULL)
39  , model_(NULL)
40  , rowStart_(NULL)
41  , element_(NULL)
42  , column_(NULL)
43  , numberColumnBlocks_(0)
44  , numberRowBlocks_(0)
45  ,
46#ifdef COUNT_COPY
47  countRealColumn_(NULL)
48  , countStartLarge_(NULL)
49  , countRow_(NULL)
50  , countElement_(NULL)
51  , smallestCount_(0)
52  , largestCount_(0)
53  ,
54#endif
55  startFraction_(0.0)
56  , endFraction_(1.0)
57  , savedBestDj_(0.0)
58  , originalWanted_(0)
59  , currentWanted_(0)
60  , savedBestSequence_(-1)
61  , minimumObjectsScan_(-1)
62  , minimumGoodReducedCosts_(-1)
63{
64}
65
66//-------------------------------------------------------------------
67// Copy constructor
68//-------------------------------------------------------------------
69AbcMatrix::AbcMatrix(const AbcMatrix &rhs)
70{
71#ifndef COIN_SPARSE_MATRIX
72  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -1, -1);
73#else
74  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
75#endif
76  model_ = rhs.model_;
77  rowStart_ = NULL;
78  element_ = NULL;
79  column_ = NULL;
80#ifdef COUNT_COPY
81  countRealColumn_ = NULL;
82  countStartLarge_ = NULL;
83  countRow_ = NULL;
84  countElement_ = NULL;
85#endif
86  numberColumnBlocks_ = rhs.numberColumnBlocks_;
87  CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_, numberColumnBlocks_ + 1);
88  numberRowBlocks_ = rhs.numberRowBlocks_;
89  if (numberRowBlocks_) {
90    assert(model_);
91    int numberRows = model_->numberRows();
92    int numberElements = matrix_->getNumElements();
93    memcpy(blockStart_, rhs.blockStart_, sizeof(blockStart_));
94    rowStart_ = CoinCopyOfArray(rhs.rowStart_, numberRows * (numberRowBlocks_ + 2));
95    element_ = CoinCopyOfArray(rhs.element_, numberElements);
96    column_ = CoinCopyOfArray(rhs.column_, numberElements);
97#ifdef COUNT_COPY
98    smallestCount_ = rhs.smallestCount_;
99    largestCount_ = rhs.largestCount_;
100    int numberColumns = model_->numberColumns();
101    countRealColumn_ = CoinCopyOfArray(rhs.countRealColumn_, numberColumns);
102    memcpy(countStart_, rhs.countStart_, reinterpret_cast< char * >(&countRealColumn_) - reinterpret_cast< char * >(countStart_));
103    int numberLarge = numberColumns - countStart_[MAX_COUNT];
104    countStartLarge_ = CoinCopyOfArray(rhs.countStartLarge_, numberLarge + 1);
105    numberElements = countStartLarge_[numberLarge];
106    countElement_ = CoinCopyOfArray(rhs.countElement_, numberElements);
107    countRow_ = CoinCopyOfArray(rhs.countRow_, numberElements);
108#endif
109  }
110}
111
112AbcMatrix::AbcMatrix(const CoinPackedMatrix &rhs)
113{
114#ifndef COIN_SPARSE_MATRIX
115  matrix_ = new CoinPackedMatrix(rhs, -1, -1);
116#else
117  matrix_ = new CoinPackedMatrix(rhs, -0, -0);
118#endif
119  matrix_->cleanMatrix();
120  model_ = NULL;
121  rowStart_ = NULL;
122  element_ = NULL;
123  column_ = NULL;
124#ifdef COUNT_COPY
125  countRealColumn_ = NULL;
126  countStartLarge_ = NULL;
127  countRow_ = NULL;
128  countElement_ = NULL;
129  smallestCount_ = 0;
130  largestCount_ = 0;
131#endif
132  numberColumnBlocks_ = 1;
133  startColumnBlock_[0] = 0;
134  startColumnBlock_[1] = 0;
135  numberRowBlocks_ = 0;
136  startFraction_ = 0;
137  endFraction_ = 1.0;
138  savedBestDj_ = 0;
139  originalWanted_ = 0;
140  currentWanted_ = 0;
141  savedBestSequence_ = -1;
142  minimumObjectsScan_ = -1;
143  minimumGoodReducedCosts_ = -1;
144}
145#ifdef ABC_SPRINT
146/* Subset constructor (without gaps). */
147AbcMatrix::AbcMatrix(const AbcMatrix &wholeMatrix,
148  int numberRows, const int *whichRows,
149  int numberColumns, const int *whichColumns)
150{
151#ifndef COIN_SPARSE_MATRIX
152  matrix_ = new CoinPackedMatrix(*wholeMatrix.matrix_, numberRows, whichRows,
153    numberColumns, whichColumns);
154#else
155  matrix_ = new CoinPackedMatrix(rhs, -0, -0);
156  abort();
157#endif
158  matrix_->cleanMatrix();
159  model_ = NULL;
160  rowStart_ = NULL;
161  element_ = NULL;
162  column_ = NULL;
163#ifdef COUNT_COPY
164  countRealColumn_ = NULL;
165  countStartLarge_ = NULL;
166  countRow_ = NULL;
167  countElement_ = NULL;
168  smallestCount_ = 0;
169  largestCount_ = 0;
170#endif
171  numberColumnBlocks_ = 1;
172  startColumnBlock_[0] = 0;
173  startColumnBlock_[1] = 0;
174  numberRowBlocks_ = 0;
175  startFraction_ = 0;
176  endFraction_ = 1.0;
177  savedBestDj_ = 0;
178  originalWanted_ = 0;
179  currentWanted_ = 0;
180  savedBestSequence_ = -1;
181  minimumObjectsScan_ = -1;
182  minimumGoodReducedCosts_ = -1;
183}
184#endif
185//-------------------------------------------------------------------
186// Destructor
187//-------------------------------------------------------------------
188AbcMatrix::~AbcMatrix()
189{
190  delete matrix_;
191  delete[] rowStart_;
192  delete[] element_;
193  delete[] column_;
194#ifdef COUNT_COPY
195  delete[] countRealColumn_;
196  delete[] countStartLarge_;
197  delete[] countRow_;
198  delete[] countElement_;
199#endif
200}
201
202//----------------------------------------------------------------
203// Assignment operator
204//-------------------------------------------------------------------
205AbcMatrix &
206AbcMatrix::operator=(const AbcMatrix &rhs)
207{
208  if (this != &rhs) {
209    delete matrix_;
210#ifndef COIN_SPARSE_MATRIX
211    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
212#else
213    matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
214#endif
215    model_ = rhs.model_;
216    delete[] rowStart_;
217    delete[] element_;
218    delete[] column_;
219#ifdef COUNT_COPY
220    delete[] countRealColumn_;
221    delete[] countStartLarge_;
222    delete[] countRow_;
223    delete[] countElement_;
224#endif
225    rowStart_ = NULL;
226    element_ = NULL;
227    column_ = NULL;
228#ifdef COUNT_COPY
229    countRealColumn_ = NULL;
230    countStartLarge_ = NULL;
231    countRow_ = NULL;
232    countElement_ = NULL;
233#endif
234    numberColumnBlocks_ = rhs.numberColumnBlocks_;
235    CoinAbcMemcpy(startColumnBlock_, rhs.startColumnBlock_, numberColumnBlocks_ + 1);
236    numberRowBlocks_ = rhs.numberRowBlocks_;
237    if (numberRowBlocks_) {
238      assert(model_);
239      int numberRows = model_->numberRows();
240      int numberElements = matrix_->getNumElements();
241      memcpy(blockStart_, rhs.blockStart_, sizeof(blockStart_));
242      rowStart_ = CoinCopyOfArray(rhs.rowStart_, numberRows * (numberRowBlocks_ + 2));
243      element_ = CoinCopyOfArray(rhs.element_, numberElements);
244      column_ = CoinCopyOfArray(rhs.column_, numberElements);
245#ifdef COUNT_COPY
246      smallestCount_ = rhs.smallestCount_;
247      largestCount_ = rhs.largestCount_;
248      int numberColumns = model_->numberColumns();
249      countRealColumn_ = CoinCopyOfArray(rhs.countRealColumn_, numberColumns);
250      memcpy(countStart_, rhs.countStart_, reinterpret_cast< char * >(&countRealColumn_) - reinterpret_cast< char * >(countStart_));
251      int numberLarge = numberColumns - countStart_[MAX_COUNT];
252      countStartLarge_ = CoinCopyOfArray(rhs.countStartLarge_, numberLarge + 1);
253      numberElements = countStartLarge_[numberLarge];
254      countElement_ = CoinCopyOfArray(rhs.countElement_, numberElements);
255      countRow_ = CoinCopyOfArray(rhs.countRow_, numberElements);
256#endif
257    }
258    startFraction_ = rhs.startFraction_;
259    endFraction_ = rhs.endFraction_;
260    savedBestDj_ = rhs.savedBestDj_;
261    originalWanted_ = rhs.originalWanted_;
262    currentWanted_ = rhs.currentWanted_;
263    savedBestSequence_ = rhs.savedBestSequence_;
264    minimumObjectsScan_ = rhs.minimumObjectsScan_;
265    minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
266  }
267  return *this;
268}
269// Sets model
270void AbcMatrix::setModel(AbcSimplex *model)
271{
272  model_ = model;
273  int numberColumns = model_->numberColumns();
274  bool needExtension = numberColumns > matrix_->getNumCols();
275  if (needExtension) {
276    CoinBigIndex lastElement = matrix_->getNumElements();
277    matrix_->reserve(numberColumns, lastElement, true);
278    CoinBigIndex *columnStart = matrix_->getMutableVectorStarts();
279    for (int i = numberColumns; i >= 0; i--) {
280      if (columnStart[i] == 0)
281        columnStart[i] = lastElement;
282      else
283        break;
284    }
285    assert(lastElement == columnStart[numberColumns]);
286  }
287}
288/* Returns a new matrix in reverse order without gaps */
289CoinPackedMatrix *
290AbcMatrix::reverseOrderedCopy() const
291{
292  CoinPackedMatrix *matrix = new CoinPackedMatrix();
293  matrix->setExtraGap(0.0);
294  matrix->setExtraMajor(0.0);
295  matrix->reverseOrderedCopyOf(*matrix_);
296  return matrix;
297}
298/// returns number of elements in column part of basis,
299CoinBigIndex
300AbcMatrix::countBasis(const int *whichColumn,
301  int &numberColumnBasic)
302{
303  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
304  CoinBigIndex numberElements = 0;
305  int numberRows = model_->numberRows();
306  // just count - can be over so ignore zero problem
307  for (int i = 0; i < numberColumnBasic; i++) {
308    int iColumn = whichColumn[i] - numberRows;
309    numberElements += columnLength[iColumn];
310  }
311  return numberElements;
312}
313void AbcMatrix::fillBasis(const int *COIN_RESTRICT whichColumn,
314  int &numberColumnBasic,
315  int *COIN_RESTRICT indexRowU,
316  int *COIN_RESTRICT start,
317  int *COIN_RESTRICT rowCount,
318  int *COIN_RESTRICT columnCount,
319  CoinFactorizationDouble *COIN_RESTRICT elementU)
320{
321  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
322  CoinBigIndex numberElements = start[0];
323  // fill
324  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
325  const int *COIN_RESTRICT row = matrix_->getIndices();
326  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
327  int numberRows = model_->numberRows();
328  for (int i = 0; i < numberColumnBasic; i++) {
329    int iColumn = whichColumn[i] - numberRows;
330    int length = columnLength[iColumn];
331    CoinBigIndex startThis = columnStart[iColumn];
332    columnCount[i] = length;
333    CoinBigIndex endThis = startThis + length;
334    for (CoinBigIndex j = startThis; j < endThis; j++) {
335      int iRow = row[j];
336      indexRowU[numberElements] = iRow;
337      rowCount[iRow]++;
338      assert(elementByColumn[j]);
339      elementU[numberElements++] = elementByColumn[j];
340    }
341    start[i + 1] = numberElements;
342  }
343}
344#ifdef ABC_LONG_FACTORIZATION
345void AbcMatrix::fillBasis(const int *COIN_RESTRICT whichColumn,
346  int &numberColumnBasic,
347  int *COIN_RESTRICT indexRowU,
348  int *COIN_RESTRICT start,
349  int *COIN_RESTRICT rowCount,
350  int *COIN_RESTRICT columnCount,
351  long double *COIN_RESTRICT elementU)
352{
353  const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths();
354  CoinBigIndex numberElements = start[0];
355  // fill
356  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
357  const int *COIN_RESTRICT row = matrix_->getIndices();
358  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
359  int numberRows = model_->numberRows();
360  for (int i = 0; i < numberColumnBasic; i++) {
361    int iColumn = whichColumn[i] - numberRows;
362    int length = columnLength[iColumn];
363    CoinBigIndex startThis = columnStart[iColumn];
364    columnCount[i] = length;
365    CoinBigIndex endThis = startThis + length;
366    for (CoinBigIndex j = startThis; j < endThis; j++) {
367      int iRow = row[j];
368      indexRowU[numberElements] = iRow;
369      rowCount[iRow]++;
370      assert(elementByColumn[j]);
371      elementU[numberElements++] = elementByColumn[j];
372    }
373    start[i + 1] = numberElements;
374  }
375}
376#endif
377#if 0 
378/// Move largest in column to beginning
379void
380AbcMatrix::moveLargestToStart()
381{
382  // get matrix data pointers
383  int *  COIN_RESTRICT row = matrix_->getMutableIndices();
384  const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts();
385  double *  COIN_RESTRICT elementByColumn = matrix_->getMutableElements();
386  int numberColumns=model_->numberColumns();
387  CoinBigIndex start = 0;
388  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
389    CoinBigIndex end = columnStart[iColumn+1];
390    double largest=0.0;
391    int position=-1;
392    for (CoinBigIndex j = start; j < end; j++) {
393      double value = fabs(elementByColumn[j]);
394      if (value>largest) {
395        largest=value;
396        position=j;
397      }
398    }
399    assert (position>=0); // ? empty column
400    if (position>start) {
401      double value=elementByColumn[start];
402      elementByColumn[start]=elementByColumn[position];
403      elementByColumn[position]=value;
404      int iRow=row[start];
405      row[start]=row[position];
406      row[position]=iRow;
407    }
408    start=end;
409  }
410}
411#endif
412// Creates row copy
413void AbcMatrix::createRowCopy()
414{
415#if ABC_PARALLEL
416  if (model_->parallelMode() == 0)
417#endif
418    numberRowBlocks_ = 1;
419#if ABC_PARALLEL
420  else
421    numberRowBlocks_ = CoinMin(NUMBER_ROW_BLOCKS, model_->numberCpus());
422#endif
423  int maximumRows = model_->maximumAbcNumberRows();
424  int numberRows = model_->numberRows();
425  int numberColumns = model_->numberColumns();
426  int numberElements = matrix_->getNumElements();
427  assert(!rowStart_);
428  char *whichBlock_ = new char[numberColumns];
429  rowStart_ = new CoinBigIndex[numberRows * (numberRowBlocks_ + 2)];
430  element_ = new double[numberElements];
431  column_ = new int[numberElements];
432  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
433  memset(blockStart_, 0, sizeof(blockStart_));
434  int ecount[10];
435  assert(numberRowBlocks_ < 16);
436  CoinAbcMemset0(ecount, 10);
437  // allocate to blocks (put a bit less in first as will be dealing with slacks) LATER
438  CoinBigIndex start = 0;
439  int block = 0;
440  CoinBigIndex work = (2 * numberColumns + matrix_->getNumElements() + numberRowBlocks_ - 1) / numberRowBlocks_;
441  CoinBigIndex thisWork = work;
442  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
443    CoinBigIndex end = columnStart[iColumn + 1];
444    assert(block >= 0 && block < numberRowBlocks_);
445    whichBlock_[iColumn] = static_cast< char >(block);
446    thisWork -= 2 + end - start;
447    ecount[block] += end - start;
448    start = end;
449    blockStart_[block]++;
450    if (thisWork <= 0) {
451      block++;
452      thisWork = work;
453    }
454  }
455#if 0
456  printf("Blocks ");
457  for (int i=0;i<numberRowBlocks_;i++)
458    printf("(%d %d) ",blockStart_[i],ecount[i]);
459  printf("\n");
460#endif
461  start = 0;
462  for (int i = 0; i < numberRowBlocks_; i++) {
463    int next = start + blockStart_[i];
464    blockStart_[i] = start;
465    start = next;
466  }
467  blockStart_[numberRowBlocks_] = start;
468  assert(start == numberColumns);
469  const int *COIN_RESTRICT row = matrix_->getIndices();
470  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
471  // counts
472  CoinAbcMemset0(rowStart_, numberRows * (numberRowBlocks_ + 2));
473  int *COIN_RESTRICT last = rowStart_ + numberRows * (numberRowBlocks_ + 1);
474  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
475    int block = whichBlock_[iColumn];
476    blockStart_[block]++;
477    int base = (block + 1) * numberRows;
478    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
479      int iRow = row[j];
480      rowStart_[base + iRow]++;
481      last[iRow]++;
482    }
483  }
484  // back to real starts
485  for (int iBlock = numberRowBlocks_ - 1; iBlock >= 0; iBlock--)
486    blockStart_[iBlock + 1] = blockStart_[iBlock];
487  blockStart_[0] = 0;
488  CoinBigIndex put = 0;
489  for (int iRow = 1; iRow < numberRows; iRow++) {
490    int n = last[iRow - 1];
491    last[iRow - 1] = put;
492    put += n;
493    rowStart_[iRow] = put;
494  }
495  int n = last[numberRows - 1];
496  last[numberRows - 1] = put;
497  put += n;
498  assert(put == numberElements);
499  //last[numberRows-1]=put;
500  // starts
501  int base = 0;
502  for (int iBlock = 0; iBlock < numberRowBlocks_; iBlock++) {
503    for (int iRow = 0; iRow < numberRows; iRow++) {
504      rowStart_[base + numberRows + iRow] += rowStart_[base + iRow];
505    }
506    base += numberRows;
507  }
508  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
509    int block = whichBlock_[iColumn];
510    int where = blockStart_[block];
511    blockStart_[block]++;
512    int base = block * numberRows;
513    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
514      int iRow = row[j];
515      CoinBigIndex put = rowStart_[base + iRow];
516      rowStart_[base + iRow]++;
517      column_[put] = where + maximumRows;
518      element_[put] = elementByColumn[j];
519    }
520  }
521  // back to real starts etc
522  base = numberRows * numberRowBlocks_;
523  for (int iBlock = numberRowBlocks_ - 1; iBlock >= 0; iBlock--) {
524    blockStart_[iBlock + 1] = blockStart_[iBlock] + maximumRows;
525    CoinAbcMemcpy(rowStart_ + base, rowStart_ + base - numberRows, numberRows);
526    base -= numberRows;
527  }
528  blockStart_[0] = 0; //maximumRows;
529  delete[] whichBlock_;
530  // and move
531  CoinAbcMemcpy(rowStart_, last, numberRows);
532  // All in useful
533  CoinAbcMemcpy(rowStart_ + (numberRowBlocks_ + 1) * numberRows,
534    rowStart_ + (numberRowBlocks_)*numberRows, numberRows);
535#ifdef COUNT_COPY
536  // now blocked by element count
537  countRealColumn_ = new int[numberColumns];
538  int counts[2 * MAX_COUNT];
539  memset(counts, 0, sizeof(counts));
540  //memset(countFirst_,0,sizeof(countFirst_));
541  int numberLarge = 0;
542  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
543    int n = columnStart[iColumn + 1] - columnStart[iColumn];
544    if (n < MAX_COUNT)
545      counts[n]++;
546    else
547      numberLarge++;
548  }
549  CoinBigIndex numberExtra = 3; // for alignment
550#define SMALL_COUNT 1
551  for (int i = 0; i < MAX_COUNT; i++) {
552    int n = counts[i];
553    if (n >= SMALL_COUNT) {
554      n &= 3;
555      int extra = (4 - n) & 3;
556      numberExtra += i * extra;
557    } else {
558      // treat as large
559      numberLarge += n;
560    }
561  }
562  countElement_ = new double[numberElements + numberExtra];
563  countRow_ = new int[numberElements + numberExtra];
564  countStartLarge_ = new CoinBigIndex[numberLarge + 1];
565  countStartLarge_[numberLarge] = numberElements + numberExtra;
566  //return;
567  CoinInt64 xx = reinterpret_cast< CoinInt64 >(countElement_);
568  int iBottom = static_cast< int >(xx & 31);
569  int offset = iBottom >> 3;
570  CoinBigIndex firstElementLarge = 0;
571  if (offset)
572    firstElementLarge = 4 - offset;
573  //countStart_[0]=firstElementLarge;
574  int positionLarge = 0;
575  smallestCount_ = 0;
576  largestCount_ = 0;
577  for (int i = 0; i < MAX_COUNT; i++) {
578    int n = counts[i];
579    countFirst_[i] = positionLarge;
580    countStart_[i] = firstElementLarge;
581    if (n >= SMALL_COUNT) {
582      counts[i + MAX_COUNT] = 1;
583      if (smallestCount_ == 0)
584        smallestCount_ = i;
585      largestCount_ = i;
586      positionLarge += n;
587      firstElementLarge += n * i;
588      n &= 3;
589      int extra = (4 - n) & 3;
590      firstElementLarge += i * extra;
591    }
592    counts[i] = 0;
593  }
594  largestCount_++;
595  countFirst_[MAX_COUNT] = positionLarge;
596  countStart_[MAX_COUNT] = firstElementLarge;
597  numberLarge = 0;
598  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
599    CoinBigIndex start = columnStart[iColumn];
600    int n = columnStart[iColumn + 1] - start;
601    CoinBigIndex put;
602    int position;
603    if (n < MAX_COUNT && counts[MAX_COUNT + n] != 0) {
604      int iWhich = counts[n];
605      position = countFirst_[n] + iWhich;
606      int iBlock4 = iWhich & (~3);
607      iWhich &= 3;
608      put = countStart_[n] + iBlock4 * n;
609      put += iWhich;
610      counts[n]++;
611      for (int i = 0; i < n; i++) {
612        countRow_[put] = row[start + i];
613        countElement_[put] = elementByColumn[start + i];
614        put += 4;
615      }
616    } else {
617      position = positionLarge + numberLarge;
618      put = firstElementLarge;
619      countStartLarge_[numberLarge] = put;
620      firstElementLarge += n;
621      numberLarge++;
622      CoinAbcMemcpy(countRow_ + put, row + start, n);
623      CoinAbcMemcpy(countElement_ + put, elementByColumn + start, n);
624    }
625    countRealColumn_[position] = iColumn + maximumRows;
626  }
627  countStartLarge_[numberLarge] = firstElementLarge;
628  assert(firstElementLarge <= numberElements + numberExtra);
629#endif
630}
631// Make all useful
632void AbcMatrix::makeAllUseful(CoinIndexedVector & /*spare*/)
633{
634  int numberRows = model_->numberRows();
635  CoinBigIndex *COIN_RESTRICT rowEnd = rowStart_ + numberRows;
636  const CoinBigIndex *COIN_RESTRICT rowReallyEnd = rowStart_ + 2 * numberRows;
637  for (int iRow = 0; iRow < numberRows; iRow++) {
638    rowEnd[iRow] = rowReallyEnd[iRow];
639  }
640}
641#define SQRT_ARRAY
642// Creates scales for column copy (rowCopy in model may be modified)
643void AbcMatrix::scale(int numberAlreadyScaled)
644{
645  int numberRows = model_->numberRows();
646  bool inBranchAndBound = (model_->specialOptions(), 0x1000000) != 0;
647  bool doScaling = numberAlreadyScaled >= 0;
648  if (!doScaling)
649    numberAlreadyScaled = 0;
650  if (numberAlreadyScaled == numberRows)
651    return; // no need to do anything
652  int numberColumns = model_->numberColumns();
653  double *COIN_RESTRICT rowScale = model_->rowScale2();
654  double *COIN_RESTRICT inverseRowScale = model_->inverseRowScale2();
655  double *COIN_RESTRICT columnScale = model_->columnScale2();
656  double *COIN_RESTRICT inverseColumnScale = model_->inverseColumnScale2();
657  // we are going to mark bits we are interested in
658  int whichArray = model_->getAvailableArrayPublic();
659  char *COIN_RESTRICT usefulColumn = reinterpret_cast< char * >(model_->usefulArray(whichArray)->getIndices());
660  memset(usefulColumn, 1, numberColumns);
661  const double *COIN_RESTRICT rowLower = model_->rowLower();
662  const double *COIN_RESTRICT rowUpper = model_->rowUpper();
663  const double *COIN_RESTRICT columnLower = model_->columnLower();
664  const double *COIN_RESTRICT columnUpper = model_->columnUpper();
665  //#define LEAVE_FIXED
666  // mark empty and fixed columns
667  // get matrix data pointers
668  const int *COIN_RESTRICT row = matrix_->getIndices();
669  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
670  double *COIN_RESTRICT elementByColumn = matrix_->getMutableElements();
671  CoinPackedMatrix *COIN_RESTRICT rowCopy = reverseOrderedCopy();
672  const int *column = rowCopy->getIndices();
673  const CoinBigIndex *COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
674  const double *COIN_RESTRICT element = rowCopy->getElements();
675  assert(numberAlreadyScaled >= 0 && numberAlreadyScaled < numberRows);
676#ifndef LEAVE_FIXED
677  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
678    if (columnUpper[iColumn] == columnLower[iColumn])
679      usefulColumn[iColumn] = 0;
680  }
681#endif
682  double overallLargest = -1.0e-20;
683  double overallSmallest = 1.0e20;
684  if (!numberAlreadyScaled) {
685    // may be redundant
686    CoinFillN(rowScale, numberRows, 1.0);
687    CoinFillN(columnScale, numberColumns, 1.0);
688    CoinBigIndex start = 0;
689    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
690      CoinBigIndex end = columnStart[iColumn + 1];
691      if (usefulColumn[iColumn]) {
692        if (end > start) {
693          for (CoinBigIndex j = start; j < end; j++) {
694            double value = fabs(elementByColumn[j]);
695            overallLargest = CoinMax(overallLargest, value);
696            overallSmallest = CoinMin(overallSmallest, value);
697          }
698        } else {
699          usefulColumn[iColumn] = 0;
700        }
701      }
702      start = end;
703    }
704  } else {
705    CoinBigIndex start = rowStart[numberAlreadyScaled];
706    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
707      rowScale[iRow] = 1.0;
708      CoinBigIndex end = rowStart[iRow + 1];
709      for (CoinBigIndex j = start; j < end; j++) {
710        int iColumn = column[j];
711        if (usefulColumn[iColumn]) {
712          double value = fabs(elementByColumn[j]) * columnScale[iColumn];
713          overallLargest = CoinMax(overallLargest, value);
714          overallSmallest = CoinMin(overallSmallest, value);
715        }
716      }
717    }
718  }
719  if ((overallSmallest >= 0.5 && overallLargest <= 2.0) || !doScaling) {
720    //printf("no scaling\n");
721    delete rowCopy;
722    model_->clearArraysPublic(whichArray);
723    CoinFillN(inverseRowScale + numberAlreadyScaled, numberRows - numberAlreadyScaled, 1.0);
724    if (!numberAlreadyScaled)
725      CoinFillN(inverseColumnScale, numberColumns, 1.0);
726    //moveLargestToStart();
727    return;
728  }
729  // need to scale
730  double largest;
731  double smallest;
732  int scalingMethod = model_->scalingFlag();
733  if (scalingMethod == 4) {
734    // As auto
735    scalingMethod = 3;
736  } else if (scalingMethod == 5) {
737    // As geometric
738    scalingMethod = 2;
739  }
740  double savedOverallRatio = 0.0;
741  double tolerance = 5.0 * model_->primalTolerance();
742  bool finished = false;
743  // if scalingMethod 3 then may change
744  bool extraDetails = (model_->logLevel() > 2);
745  bool secondTime = false;
746  while (!finished) {
747    int numberPass = !numberAlreadyScaled ? 3 : 1;
748    overallLargest = -1.0e-20;
749    overallSmallest = 1.0e20;
750    if (!secondTime) {
751      secondTime = true;
752    } else {
753      CoinFillN(rowScale, numberRows, 1.0);
754      CoinFillN(columnScale, numberColumns, 1.0);
755    }
756    if (scalingMethod == 1 || scalingMethod == 3) {
757      // Maximum in each row
758      for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
759        largest = 1.0e-10;
760        for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
761          int iColumn = column[j];
762          if (usefulColumn[iColumn]) {
763            double value = fabs(element[j]);
764            largest = CoinMax(largest, value);
765            assert(largest < 1.0e40);
766          }
767        }
768        rowScale[iRow] = 1.0 / largest;
769#ifdef COIN_DEVELOP
770        if (extraDetails) {
771          overallLargest = CoinMax(overallLargest, largest);
772          overallSmallest = CoinMin(overallSmallest, largest);
773        }
774#endif
775      }
776    } else {
777      while (numberPass) {
778        overallLargest = 0.0;
779        overallSmallest = 1.0e50;
780        numberPass--;
781        // Geometric mean on row scales
782        for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
783          largest = 1.0e-50;
784          smallest = 1.0e50;
785          for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) {
786            int iColumn = column[j];
787            if (usefulColumn[iColumn]) {
788              double value = fabs(element[j]);
789              value *= columnScale[iColumn];
790              largest = CoinMax(largest, value);
791              smallest = CoinMin(smallest, value);
792            }
793          }
794#ifdef SQRT_ARRAY
795          rowScale[iRow] = smallest * largest;
796#else
797          rowScale[iRow] = 1.0 / sqrt(smallest * largest);
798#endif
799          if (extraDetails) {
800            overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
801            overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
802          }
803        }
804        if (model_->scalingFlag() == 5)
805          break; // just scale rows
806#ifdef SQRT_ARRAY
807        CoinAbcInverseSqrts(rowScale, numberRows);
808#endif
809        if (!inBranchAndBound)
810          model_->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model_->messagesPointer())
811            << overallSmallest
812            << overallLargest
813            << CoinMessageEol;
814        // skip last column round
815        if (numberPass == 1)
816          break;
817        // Geometric mean on column scales
818        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
819          if (usefulColumn[iColumn]) {
820            largest = 1.0e-50;
821            smallest = 1.0e50;
822            for (CoinBigIndex j = columnStart[iColumn];
823                 j < columnStart[iColumn + 1]; j++) {
824              int iRow = row[j];
825              double value = fabs(elementByColumn[j]);
826              value *= rowScale[iRow];
827              largest = CoinMax(largest, value);
828              smallest = CoinMin(smallest, value);
829            }
830#ifdef SQRT_ARRAY
831            columnScale[iColumn] = smallest * largest;
832#else
833            columnScale[iColumn] = 1.0 / sqrt(smallest * largest);
834#endif
835          }
836        }
837#ifdef SQRT_ARRAY
838        CoinAbcInverseSqrts(columnScale, numberColumns);
839#endif
840      }
841    }
842    // If ranges will make horrid then scale
843    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
844      double difference = rowUpper[iRow] - rowLower[iRow];
845      double scaledDifference = difference * rowScale[iRow];
846      if (scaledDifference > tolerance && scaledDifference < 1.0e-4) {
847        // make gap larger
848        rowScale[iRow] *= 1.0e-4 / scaledDifference;
849        rowScale[iRow] = CoinMax(1.0e-10, CoinMin(1.0e10, rowScale[iRow]));
850        //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
851        // scaledDifference,difference*rowScale[iRow]);
852      }
853    }
854    // final pass to scale columns so largest is reasonable
855    // See what smallest will be if largest is 1.0
856    if (model_->scalingFlag() != 5) {
857      overallSmallest = 1.0e50;
858      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
859        if (usefulColumn[iColumn]) {
860          largest = 1.0e-20;
861          smallest = 1.0e50;
862          for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
863            int iRow = row[j];
864            double value = fabs(elementByColumn[j] * rowScale[iRow]);
865            largest = CoinMax(largest, value);
866            smallest = CoinMin(smallest, value);
867          }
868          if (overallSmallest * largest > smallest)
869            overallSmallest = smallest / largest;
870        }
871      }
872    }
873    if (scalingMethod == 1 || scalingMethod == 2) {
874      finished = true;
875    } else if (savedOverallRatio == 0.0 && scalingMethod != 4) {
876      savedOverallRatio = overallSmallest;
877      scalingMethod = 4;
878    } else {
879      assert(scalingMethod == 4);
880      if (overallSmallest > 2.0 * savedOverallRatio) {
881        finished = true; // geometric was better
882        if (model_->scalingFlag() == 4) {
883          // If in Branch and bound change
884          if ((model_->specialOptions() & 1024) != 0) {
885            model_->scaling(2);
886          }
887        }
888      } else {
889        scalingMethod = 1; // redo equilibrium
890        if (model_->scalingFlag() == 4) {
891          // If in Branch and bound change
892          if ((model_->specialOptions() & 1024) != 0) {
893            model_->scaling(1);
894          }
895        }
896      }
897#if 0
898      if (extraDetails) {
899        if (finished)
900          printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
901                 savedOverallRatio, overallSmallest);
902        else
903          printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
904                 savedOverallRatio, overallSmallest);
905      }
906#endif
907    }
908  }
909  //#define RANDOMIZE
910#ifdef RANDOMIZE
911  // randomize by up to 10%
912  for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
913    double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
914    rowScale[iRow] *= (1.0 + 0.1 * value);
915  }
916#endif
917  overallLargest = 1.0;
918  if (overallSmallest < 1.0e-1)
919    overallLargest = 1.0 / sqrt(overallSmallest);
920  overallLargest = CoinMin(100.0, overallLargest);
921  overallSmallest = 1.0e50;
922  char *COIN_RESTRICT usedRow = reinterpret_cast< char * >(inverseRowScale);
923  memset(usedRow, 0, numberRows);
924  //printf("scaling %d\n",model_->scalingFlag());
925  if (model_->scalingFlag() != 5) {
926    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
927      if (usefulColumn[iColumn]) {
928        largest = 1.0e-20;
929        smallest = 1.0e50;
930        for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
931          int iRow = row[j];
932          usedRow[iRow] = 1;
933          double value = fabs(elementByColumn[j] * rowScale[iRow]);
934          largest = CoinMax(largest, value);
935          smallest = CoinMin(smallest, value);
936        }
937        columnScale[iColumn] = overallLargest / largest;
938        //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
939#ifdef RANDOMIZE
940        if (!numberAlreadyScaled) {
941          double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
942          columnScale[iColumn] *= (1.0 + 0.1 * value);
943        }
944#endif
945        double difference = columnUpper[iColumn] - columnLower[iColumn];
946        if (difference < 1.0e-5 * columnScale[iColumn]) {
947          // make gap larger
948          columnScale[iColumn] = difference / 1.0e-5;
949          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
950          // scaledDifference,difference*columnScale[iColumn]);
951        }
952        double value = smallest * columnScale[iColumn];
953        if (overallSmallest > value)
954          overallSmallest = value;
955        //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
956      } else {
957        assert(columnScale[iColumn] == 1.0);
958        //columnScale[iColumn]=1.0;
959      }
960    }
961    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
962      if (!usedRow[iRow]) {
963        rowScale[iRow] = 1.0;
964      }
965    }
966  }
967  if (!inBranchAndBound)
968    model_->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model_->messagesPointer())
969      << overallSmallest
970      << overallLargest
971      << CoinMessageEol;
972  if (overallSmallest < 1.0e-13) {
973    // Change factorization zero tolerance
974    double newTolerance = CoinMax(1.0e-15 * (overallSmallest / 1.0e-13),
975      1.0e-18);
976    if (model_->factorization()->zeroTolerance() > newTolerance)
977      model_->factorization()->zeroTolerance(newTolerance);
978    newTolerance = CoinMax(overallSmallest * 0.5, 1.0e-18);
979    model_->setZeroTolerance(newTolerance);
980#ifndef NDEBUG
981    assert(newTolerance < 0.0); // just so we can fix
982#endif
983  }
984  // make copy (could do faster by using previous values)
985  // could just do partial
986  CoinAbcReciprocal(inverseRowScale + numberAlreadyScaled, numberRows - numberAlreadyScaled,
987    rowScale + numberAlreadyScaled);
988  if (!numberAlreadyScaled)
989    CoinAbcReciprocal(inverseColumnScale, numberColumns, columnScale);
990  // Do scaled copy //NO and move largest to start
991  CoinBigIndex start = 0;
992  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
993    double scale = columnScale[iColumn];
994    CoinBigIndex end = columnStart[iColumn + 1];
995    for (CoinBigIndex j = start; j < end; j++) {
996      double value = elementByColumn[j];
997      int iRow = row[j];
998      if (iRow >= numberAlreadyScaled) {
999        value *= scale * rowScale[iRow];
1000        elementByColumn[j] = value;
1001      }
1002    }
1003    start = end;
1004  }
1005  delete rowCopy;
1006#if 0
1007  if (model_->rowCopy()) {
1008    // need to replace row by row
1009    CoinPackedMatrix * rowCopy = NULL;
1010    //static_cast< AbcMatrix*>(model_->rowCopy());
1011    double * element = rowCopy->getMutableElements();
1012    const int * column = rowCopy->getIndices();
1013    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1014    // scale row copy
1015    for (iRow = 0; iRow < numberRows; iRow++) {
1016      CoinBigIndex j;
1017      double scale = rowScale[iRow];
1018      double * elementsInThisRow = element + rowStart[iRow];
1019      const int * columnsInThisRow = column + rowStart[iRow];
1020      int number = rowStart[iRow+1] - rowStart[iRow];
1021      assert (number <= numberColumns);
1022      for (j = 0; j < number; j++) {
1023        int iColumn = columnsInThisRow[j];
1024        elementsInThisRow[j] *= scale * columnScale[iColumn];
1025      }
1026    }
1027  }
1028#endif
1029  model_->clearArraysPublic(whichArray);
1030}
1031/* Return <code>y + A * scalar *x</code> in <code>y</code>.
1032   @pre <code>x</code> must be of size <code>numColumns()</code>
1033   @pre <code>y</code> must be of size <code>numRows()</code> */
1034//scaled versions
1035void AbcMatrix::timesModifyExcludingSlacks(double scalar,
1036  const double *x, double *y) const
1037{
1038  int numberTotal = model_->numberTotal();
1039  int maximumRows = model_->maximumAbcNumberRows();
1040  // get matrix data pointers
1041  const int *COIN_RESTRICT row = matrix_->getIndices();
1042  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
1043  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
1044  for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1045    double value = x[iColumn];
1046    if (value) {
1047      CoinBigIndex start = columnStart[iColumn];
1048      CoinBigIndex end = columnStart[iColumn + 1];
1049      value *= scalar;
1050      for (CoinBigIndex j = start; j < end; j++) {
1051        int iRow = row[j];
1052        y[iRow] += value * elementByColumn[j];
1053      }
1054    }
1055  }
1056}
1057/* Return <code>y + A * scalar(+-1) *x</code> in <code>y</code>.
1058   @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
1059   @pre <code>y</code> must be of size <code>numRows()</code> */
1060void AbcMatrix::timesModifyIncludingSlacks(double scalar,
1061  const double *x, double *y) const
1062{
1063  int numberRows = model_->numberRows();
1064  int numberTotal = model_->numberTotal();
1065  int maximumRows = model_->maximumAbcNumberRows();
1066  // makes no sense for x==y??
1067  assert(x != y);
1068  // For now just by column and assumes already scaled (use reallyScale)
1069  // get matrix data pointers
1070  const int *COIN_RESTRICT row = matrix_->getIndices();
1071  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
1072  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
1073  if (scalar == 1.0) {
1074    // add
1075    if (x != y) {
1076      for (int i = 0; i < numberRows; i++)
1077        y[i] += x[i];
1078    } else {
1079      for (int i = 0; i < numberRows; i++)
1080        y[i] += y[i];
1081    }
1082    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1083      double value = x[iColumn];
1084      if (value) {
1085        CoinBigIndex start = columnStart[iColumn];
1086        CoinBigIndex next = columnStart[iColumn + 1];
1087        for (CoinBigIndex j = start; j < next; j++) {
1088          int jRow = row[j];
1089          y[jRow] += value * elementByColumn[j];
1090        }
1091      }
1092    }
1093  } else {
1094    // subtract
1095    if (x != y) {
1096      for (int i = 0; i < numberRows; i++)
1097        y[i] -= x[i];
1098    } else {
1099      for (int i = 0; i < numberRows; i++)
1100        y[i] = 0.0;
1101    }
1102    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1103      double value = x[iColumn];
1104      if (value) {
1105        CoinBigIndex start = columnStart[iColumn];
1106        CoinBigIndex next = columnStart[iColumn + 1];
1107        for (CoinBigIndex j = start; j < next; j++) {
1108          int jRow = row[j];
1109          y[jRow] -= value * elementByColumn[j];
1110        }
1111      }
1112    }
1113  }
1114}
1115/* Return A * scalar(+-1) *x</code> in <code>y</code>.
1116   @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
1117   @pre <code>y</code> must be of size <code>numRows()</code> */
1118void AbcMatrix::timesIncludingSlacks(double scalar,
1119  const double *x, double *y) const
1120{
1121  int numberRows = model_->numberRows();
1122  int numberTotal = model_->numberTotal();
1123  int maximumRows = model_->maximumAbcNumberRows();
1124  // For now just by column and assumes already scaled (use reallyScale)
1125  // get matrix data pointers
1126  const int *COIN_RESTRICT row = matrix_->getIndices();
1127  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
1128  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
1129  if (scalar == 1.0) {
1130    // add
1131    if (x != y) {
1132      for (int i = 0; i < numberRows; i++)
1133        y[i] = x[i];
1134    }
1135    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1136      double value = x[iColumn];
1137      if (value) {
1138        CoinBigIndex start = columnStart[iColumn];
1139        CoinBigIndex next = columnStart[iColumn + 1];
1140        for (CoinBigIndex j = start; j < next; j++) {
1141          int jRow = row[j];
1142          y[jRow] += value * elementByColumn[j];
1143        }
1144      }
1145    }
1146  } else {
1147    // subtract
1148    if (x != y) {
1149      for (int i = 0; i < numberRows; i++)
1150        y[i] = -x[i];
1151    } else {
1152      for (int i = 0; i < numberRows; i++)
1153        y[i] = -y[i];
1154    }
1155    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1156      double value = x[iColumn];
1157      if (value) {
1158        CoinBigIndex start = columnStart[iColumn];
1159        CoinBigIndex next = columnStart[iColumn + 1];
1160        for (CoinBigIndex j = start; j < next; j++) {
1161          int jRow = row[j];
1162          y[jRow] -= value * elementByColumn[j];
1163        }
1164      }
1165    }
1166  }
1167}
1168extern double parallelDual4(AbcSimplexDual *);
1169static double firstPass(AbcSimplex *model, int iBlock,
1170  double upperTheta, int &freeSequence,
1171  CoinPartitionedVector &tableauRow,
1172  CoinPartitionedVector &candidateList)
1173{
1174  int *COIN_RESTRICT index = tableauRow.getIndices();
1175  double *COIN_RESTRICT array = tableauRow.denseVector();
1176  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
1177  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
1178  const double *COIN_RESTRICT abcDj = model->djRegion();
1179  const unsigned char *COIN_RESTRICT internalStatus = model->internalStatus();
1180  // do first pass to get possibles
1181  double bestPossible = 0.0;
1182  // We can also see if infeasible or pivoting on free
1183  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1184  double acceptablePivot = model->currentAcceptablePivot();
1185  double dualT = -model->currentDualTolerance();
1186  // fixed will have been taken out by now
1187  const double multiplier[] = { 1.0, -1.0 };
1188  freeSequence = -1;
1189  int firstIn = model->abcMatrix()->blockStart(iBlock);
1190  int numberNonZero = tableauRow.getNumElements(iBlock) + firstIn;
1191  int numberRemaining = firstIn;
1192  //first=tableauRow.getNumElements();
1193  // could pass in last and if numberNonZero==last-firstIn scan as well
1194  if (model->ordinaryVariables()) {
1195    for (int i = firstIn; i < numberNonZero; i++) {
1196      int iSequence = index[i];
1197      double tableauValue = array[i];
1198      unsigned char iStatus = internalStatus[iSequence] & 7;
1199      double mult = multiplier[iStatus];
1200      double alpha = tableauValue * mult;
1201      double oldValue = abcDj[iSequence] * mult;
1202      double value = oldValue - tentativeTheta * alpha;
1203      if (value < dualT) {
1204        bestPossible = CoinMax(bestPossible, alpha);
1205        value = oldValue - upperTheta * alpha;
1206        if (value < dualT && alpha >= acceptablePivot) {
1207          upperTheta = (oldValue - dualT) / alpha;
1208        }
1209        // add to list
1210        arrayCandidate[numberRemaining] = alpha;
1211        indexCandidate[numberRemaining++] = iSequence;
1212      }
1213    }
1214  } else {
1215    double badFree = 0.0;
1216    double freeAlpha = model->currentAcceptablePivot();
1217    int freeSequenceIn = model->freeSequenceIn();
1218    double currentDualTolerance = model->currentDualTolerance();
1219    for (int i = firstIn; i < numberNonZero; i++) {
1220      int iSequence = index[i];
1221      double tableauValue = array[i];
1222      unsigned char iStatus = internalStatus[iSequence] & 7;
1223      if ((iStatus & 4) == 0) {
1224        double mult = multiplier[iStatus];
1225        double alpha = tableauValue * mult;
1226        double oldValue = abcDj[iSequence] * mult;
1227        double value = oldValue - tentativeTheta * alpha;
1228        if (value < dualT) {
1229          bestPossible = CoinMax(bestPossible, alpha);
1230          value = oldValue - upperTheta * alpha;
1231          if (value < dualT && alpha >= acceptablePivot) {
1232            upperTheta = (oldValue - dualT) / alpha;
1233          }
1234          // add to list
1235          arrayCandidate[numberRemaining] = alpha;
1236          indexCandidate[numberRemaining++] = iSequence;
1237        }
1238      } else {
1239        bool keep;
1240        bestPossible = CoinMax(bestPossible, fabs(tableauValue));
1241        double oldValue = abcDj[iSequence];
1242        // If free has to be very large - should come in via dualRow
1243        //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
1244        //break;
1245        if (oldValue > currentDualTolerance) {
1246          keep = true;
1247        } else if (oldValue < -currentDualTolerance) {
1248          keep = true;
1249        } else {
1250          if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
1251            keep = true;
1252          } else {
1253            keep = false;
1254            badFree = CoinMax(badFree, fabs(tableauValue));
1255          }
1256        }
1257        if (keep) {
1258#ifdef PAN
1259          if (model->fakeSuperBasic(iSequence) >= 0) {
1260#endif
1261            if (iSequence == freeSequenceIn)
1262              tableauValue = COIN_DBL_MAX;
1263            // free - choose largest
1264            if (fabs(tableauValue) > fabs(freeAlpha)) {
1265              freeAlpha = tableauValue;
1266              freeSequence = iSequence;
1267            }
1268#ifdef PAN
1269          }
1270#endif
1271        }
1272      }
1273    }
1274  }
1275  //firstInX=numberNonZero-firstIn;
1276  //lastInX=-1;//numberRemaining-lastInX;
1277  tableauRow.setNumElementsPartition(iBlock, numberNonZero - firstIn);
1278  candidateList.setNumElementsPartition(iBlock, numberRemaining - firstIn);
1279  return upperTheta;
1280}
1281// gets sorted tableau row and a possible value of theta
1282double
1283AbcMatrix::dualColumn1Row(int iBlock, double upperTheta, int &freeSequence,
1284  const CoinIndexedVector &update,
1285  CoinPartitionedVector &tableauRow,
1286  CoinPartitionedVector &candidateList) const
1287{
1288  int maximumRows = model_->maximumAbcNumberRows();
1289  int number = update.getNumElements();
1290  const double *COIN_RESTRICT pi = update.denseVector();
1291  const int *COIN_RESTRICT piIndex = update.getIndices();
1292  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
1293  int numberRows = model_->numberRows();
1294  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
1295  // count down
1296  int nColumns;
1297  int firstIn = blockStart_[iBlock];
1298  int first = firstIn;
1299  if (!first)
1300    first = maximumRows;
1301  int last = blockStart_[iBlock + 1];
1302  nColumns = last - first;
1303  int target = nColumns;
1304  rowStart += iBlock * numberRows;
1305  rowEnd = rowStart + numberRows;
1306  for (int i = 0; i < number; i++) {
1307    int iRow = piIndex[i];
1308    CoinBigIndex end = rowEnd[iRow];
1309    target -= end - rowStart[iRow];
1310    if (target < 0)
1311      break;
1312  }
1313  if (target > 0) {
1314    //printf("going to few %d ops %d\n",number,nColumns-target);
1315    return dualColumn1RowFew(iBlock, upperTheta, freeSequence,
1316      update, tableauRow, candidateList);
1317  }
1318  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
1319  int *COIN_RESTRICT index = tableauRow.getIndices();
1320  double *COIN_RESTRICT array = tableauRow.denseVector();
1321  const double *COIN_RESTRICT element = element_;
1322  const int *COIN_RESTRICT column = column_;
1323  for (int i = 0; i < number; i++) {
1324    int iRow = piIndex[i];
1325    double piValue = pi[iRow];
1326    CoinBigIndex end = rowEnd[iRow];
1327    for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
1328      int iColumn = column[j];
1329      assert(iColumn >= first && iColumn < last);
1330      double value = element[j];
1331      array[iColumn] += piValue * value;
1332    }
1333  }
1334  int numberNonZero;
1335  int numberRemaining;
1336  if (iBlock == 0) {
1337#if 1
1338    upperTheta = static_cast< AbcSimplexDual * >(model_)->dualColumn1A();
1339    numberNonZero = tableauRow.getNumElements(0);
1340    numberRemaining = candidateList.getNumElements(0);
1341#else
1342    numberNonZero = 0;
1343    for (int i = 0; i < number; i++) {
1344      int iRow = piIndex[i];
1345      unsigned char type = internalStatus[iRow];
1346      if ((type & 7) < 6) {
1347        index[numberNonZero] = iRow;
1348        double piValue = pi[iRow];
1349        array[numberNonZero++] = piValue;
1350      }
1351    }
1352    numberRemaining = 0;
1353#endif
1354  } else {
1355    numberNonZero = firstIn;
1356    numberRemaining = firstIn;
1357  }
1358  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
1359  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
1360  //printf("first %d last %d firstIn %d lastIn %d\n",
1361  //     first,last,firstIn,lastIn);
1362  const double *COIN_RESTRICT abcDj = model_->djRegion();
1363  // do first pass to get possibles
1364  double bestPossible = 0.0;
1365  // We can also see if infeasible or pivoting on free
1366  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1367  double acceptablePivot = model_->currentAcceptablePivot();
1368  double dualT = -model_->currentDualTolerance();
1369  const double multiplier[] = { 1.0, -1.0 };
1370  double zeroTolerance = model_->zeroTolerance();
1371  freeSequence = -1;
1372  if (model_->ordinaryVariables()) {
1373    for (int iSequence = first; iSequence < last; iSequence++) {
1374      double tableauValue = array[iSequence];
1375      if (tableauValue) {
1376        array[iSequence] = 0.0;
1377        if (fabs(tableauValue) > zeroTolerance) {
1378          unsigned char iStatus = internalStatus[iSequence] & 7;
1379          if (iStatus < 4) {
1380            index[numberNonZero] = iSequence;
1381            array[numberNonZero++] = tableauValue;
1382            double mult = multiplier[iStatus];
1383            double alpha = tableauValue * mult;
1384            double oldValue = abcDj[iSequence] * mult;
1385            double value = oldValue - tentativeTheta * alpha;
1386            if (value < dualT) {
1387              bestPossible = CoinMax(bestPossible, alpha);
1388              value = oldValue - upperTheta * alpha;
1389              if (value < dualT && alpha >= acceptablePivot) {
1390                upperTheta = (oldValue - dualT) / alpha;
1391              }
1392              // add to list
1393              arrayCandidate[numberRemaining] = alpha;
1394              indexCandidate[numberRemaining++] = iSequence;
1395            }
1396          }
1397        }
1398      }
1399    }
1400  } else {
1401    double badFree = 0.0;
1402    double freeAlpha = model_->currentAcceptablePivot();
1403    int freeSequenceIn = model_->freeSequenceIn();
1404    //printf("block %d freeSequence %d acceptable %g\n",iBlock,freeSequenceIn,freeAlpha);
1405    double currentDualTolerance = model_->currentDualTolerance();
1406    for (int iSequence = first; iSequence < last; iSequence++) {
1407      double tableauValue = array[iSequence];
1408      if (tableauValue) {
1409        array[iSequence] = 0.0;
1410        if (fabs(tableauValue) > zeroTolerance) {
1411          unsigned char iStatus = internalStatus[iSequence] & 7;
1412          if (iStatus < 6) {
1413            if ((iStatus & 4) == 0) {
1414              index[numberNonZero] = iSequence;
1415              array[numberNonZero++] = tableauValue;
1416              double mult = multiplier[iStatus];
1417              double alpha = tableauValue * mult;
1418              double oldValue = abcDj[iSequence] * mult;
1419              double value = oldValue - tentativeTheta * alpha;
1420              if (value < dualT) {
1421                bestPossible = CoinMax(bestPossible, alpha);
1422                value = oldValue - upperTheta * alpha;
1423                if (value < dualT && alpha >= acceptablePivot) {
1424                  upperTheta = (oldValue - dualT) / alpha;
1425                }
1426                // add to list
1427                arrayCandidate[numberRemaining] = alpha;
1428                indexCandidate[numberRemaining++] = iSequence;
1429              }
1430            } else {
1431              bool keep;
1432              index[numberNonZero] = iSequence;
1433              array[numberNonZero++] = tableauValue;
1434              bestPossible = CoinMax(bestPossible, fabs(tableauValue));
1435              double oldValue = abcDj[iSequence];
1436              // If free has to be very large - should come in via dualRow
1437              //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
1438              //break;
1439              // may be fake super basic
1440              if (oldValue > currentDualTolerance) {
1441                keep = true;
1442              } else if (oldValue < -currentDualTolerance) {
1443                keep = true;
1444              } else {
1445                if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
1446                  keep = true;
1447                } else {
1448                  keep = false;
1449                  badFree = CoinMax(badFree, fabs(tableauValue));
1450                }
1451              }
1452#if 0
1453              if (iSequence==freeSequenceIn)
1454                assert (keep);
1455#endif
1456              if (keep) {
1457#ifdef PAN
1458                if (model_->fakeSuperBasic(iSequence) >= 0) {
1459#endif
1460                  if (iSequence == freeSequenceIn)
1461                    tableauValue = COIN_DBL_MAX;
1462                  // free - choose largest
1463                  if (fabs(tableauValue) > fabs(freeAlpha)) {
1464                    freeAlpha = tableauValue;
1465                    freeSequence = iSequence;
1466                  }
1467#ifdef PAN
1468                }
1469#endif
1470              }
1471            }
1472          }
1473        }
1474      }
1475    }
1476  }
1477#if 0
1478  if (model_->freeSequenceIn()>=first&&model_->freeSequenceIn()<last)
1479    assert (freeSequence==model_->freeSequenceIn());
1480  extern int xxInfo[6][8];
1481  xxInfo[0][iBlock]=first;
1482  xxInfo[1][iBlock]=last;
1483  xxInfo[2][iBlock]=firstIn;
1484  xxInfo[3][iBlock]=numberNonZero-firstIn;
1485  xxInfo[4][iBlock]=numberRemaining-firstIn;
1486#endif
1487  tableauRow.setNumElementsPartition(iBlock, numberNonZero - firstIn);
1488  candidateList.setNumElementsPartition(iBlock, numberRemaining - firstIn);
1489  return upperTheta;
1490}
1491// gets sorted tableau row and a possible value of theta
1492double
1493AbcMatrix::dualColumn1Row2(double upperTheta, int &freeSequence,
1494  const CoinIndexedVector &update,
1495  CoinPartitionedVector &tableauRow,
1496  CoinPartitionedVector &candidateList) const
1497{
1498  //int first=model_->maximumAbcNumberRows();
1499  assert(update.getNumElements() == 2);
1500  const double *COIN_RESTRICT pi = update.denseVector();
1501  const int *COIN_RESTRICT piIndex = update.getIndices();
1502  int *COIN_RESTRICT index = tableauRow.getIndices();
1503  double *COIN_RESTRICT array = tableauRow.denseVector();
1504  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
1505  int numberRows = model_->numberRows();
1506  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
1507  const double *COIN_RESTRICT element = element_;
1508  const int *COIN_RESTRICT column = column_;
1509  int iRow0 = piIndex[0];
1510  int iRow1 = piIndex[1];
1511  CoinBigIndex end0 = rowEnd[iRow0];
1512  CoinBigIndex end1 = rowEnd[iRow1];
1513  if (end0 - rowStart[iRow0] > end1 - rowStart[iRow1]) {
1514    int temp = iRow0;
1515    iRow0 = iRow1;
1516    iRow1 = temp;
1517  }
1518  CoinBigIndex start = rowStart[iRow0];
1519  CoinBigIndex end = rowEnd[iRow0];
1520  double piValue = pi[iRow0];
1521  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
1522  int numberNonZero;
1523  numberNonZero = tableauRow.getNumElements(0);
1524  int n = numberNonZero;
1525  for (CoinBigIndex j = start; j < end; j++) {
1526    int iSequence = column[j];
1527    double value = element[j];
1528    arrayCandidate[iSequence] = piValue * value;
1529    index[n++] = iSequence;
1530  }
1531  start = rowStart[iRow1];
1532  end = rowEnd[iRow1];
1533  piValue = pi[iRow1];
1534  for (CoinBigIndex j = start; j < end; j++) {
1535    int iSequence = column[j];
1536    double value = element[j];
1537    value *= piValue;
1538    if (!arrayCandidate[iSequence]) {
1539      arrayCandidate[iSequence] = value;
1540      index[n++] = iSequence;
1541    } else {
1542      arrayCandidate[iSequence] += value;
1543    }
1544  }
1545  // pack down
1546  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
1547  double zeroTolerance = model_->zeroTolerance();
1548  while (numberNonZero < n) {
1549    int iSequence = index[numberNonZero];
1550    double value = arrayCandidate[iSequence];
1551    arrayCandidate[iSequence] = 0.0;
1552    unsigned char iStatus = internalStatus[iSequence] & 7;
1553    if (fabs(value) > zeroTolerance && iStatus < 6) {
1554      index[numberNonZero] = iSequence;
1555      array[numberNonZero++] = value;
1556    } else {
1557      // kill
1558      n--;
1559      index[numberNonZero] = index[n];
1560    }
1561  }
1562  tableauRow.setNumElementsPartition(0, numberNonZero);
1563  return firstPass(model_, 0, upperTheta, freeSequence,
1564    tableauRow,
1565    candidateList);
1566}
1567//static int ixxxxxx=1;
1568// gets sorted tableau row and a possible value of theta
1569double
1570AbcMatrix::dualColumn1RowFew(int iBlock, double upperTheta, int &freeSequence,
1571  const CoinIndexedVector &update,
1572  CoinPartitionedVector &tableauRow,
1573  CoinPartitionedVector &candidateList) const
1574{
1575  //int first=model_->maximumAbcNumberRows();
1576  int number = update.getNumElements();
1577  const double *COIN_RESTRICT pi = update.denseVector();
1578  const int *COIN_RESTRICT piIndex = update.getIndices();
1579  int *COIN_RESTRICT index = tableauRow.getIndices();
1580  double *COIN_RESTRICT array = tableauRow.denseVector();
1581  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
1582  int numberRows = model_->numberRows();
1583  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
1584  const double *COIN_RESTRICT element = element_;
1585  const int *COIN_RESTRICT column = column_;
1586  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
1587  int numberNonZero;
1588  assert(iBlock >= 0);
1589  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
1590  int firstIn = blockStart_[iBlock];
1591  if (iBlock == 0) {
1592    numberNonZero = 0;
1593    for (int i = 0; i < number; i++) {
1594      int iRow = piIndex[i];
1595      unsigned char type = internalStatus[iRow];
1596      if ((type & 7) < 6) {
1597        index[numberNonZero] = iRow;
1598        double piValue = pi[iRow];
1599        array[numberNonZero++] = piValue;
1600      }
1601    }
1602  } else {
1603    numberNonZero = firstIn;
1604  }
1605  int n = numberNonZero;
1606  rowStart += iBlock * numberRows;
1607  rowEnd = rowStart + numberRows;
1608  for (int i = 0; i < number; i++) {
1609    int iRow = piIndex[i];
1610    double piValue = pi[iRow];
1611    CoinBigIndex end = rowEnd[iRow];
1612    for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
1613      int iColumn = column[j];
1614      double value = element[j] * piValue;
1615      double oldValue = arrayCandidate[iColumn];
1616      value += oldValue;
1617      if (!oldValue) {
1618        arrayCandidate[iColumn] = value;
1619        index[n++] = iColumn;
1620      } else if (value) {
1621        arrayCandidate[iColumn] = value;
1622      } else {
1623        arrayCandidate[iColumn] = COIN_INDEXED_REALLY_TINY_ELEMENT;
1624      }
1625    }
1626  }
1627  // pack down
1628  double zeroTolerance = model_->zeroTolerance();
1629  while (numberNonZero < n) {
1630    int iSequence = index[numberNonZero];
1631    double value = arrayCandidate[iSequence];
1632    arrayCandidate[iSequence] = 0.0;
1633    unsigned char iStatus = internalStatus[iSequence] & 7;
1634    if (fabs(value) > zeroTolerance && iStatus < 6) {
1635      index[numberNonZero] = iSequence;
1636      array[numberNonZero++] = value;
1637    } else {
1638      // kill
1639      n--;
1640      index[numberNonZero] = index[n];
1641    }
1642  }
1643  tableauRow.setNumElementsPartition(iBlock, numberNonZero - firstIn);
1644  upperTheta = firstPass(model_, iBlock, upperTheta, freeSequence,
1645    tableauRow,
1646    candidateList);
1647  return upperTheta;
1648}
1649// gets sorted tableau row and a possible value of theta
1650double
1651AbcMatrix::dualColumn1Row1(double upperTheta, int &freeSequence,
1652  const CoinIndexedVector &update,
1653  CoinPartitionedVector &tableauRow,
1654  CoinPartitionedVector &candidateList) const
1655{
1656  assert(update.getNumElements() == 1);
1657  int iRow = update.getIndices()[0];
1658  double piValue = update.denseVector()[iRow];
1659  int *COIN_RESTRICT index = tableauRow.getIndices();
1660  double *COIN_RESTRICT array = tableauRow.denseVector();
1661  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_;
1662  int numberRows = model_->numberRows();
1663  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows * numberRowBlocks_;
1664  CoinBigIndex start = rowStart[iRow];
1665  CoinBigIndex end = rowEnd[iRow];
1666  const double *COIN_RESTRICT element = element_;
1667  const int *COIN_RESTRICT column = column_;
1668  int numberNonZero;
1669  int numberRemaining;
1670  numberNonZero = tableauRow.getNumElements(0);
1671  numberRemaining = candidateList.getNumElements(0);
1672  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
1673  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
1674  const double *COIN_RESTRICT abcDj = model_->djRegion();
1675  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
1676  // do first pass to get possibles
1677  double bestPossible = 0.0;
1678  // We can also see if infeasible or pivoting on free
1679  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1680  double acceptablePivot = model_->currentAcceptablePivot();
1681  double dualT = -model_->currentDualTolerance();
1682  const double multiplier[] = { 1.0, -1.0 };
1683  double zeroTolerance = model_->zeroTolerance();
1684  freeSequence = -1;
1685  if (model_->ordinaryVariables()) {
1686    for (CoinBigIndex j = start; j < end; j++) {
1687      int iSequence = column[j];
1688      double value = element[j];
1689      double tableauValue = piValue * value;
1690      if (fabs(tableauValue) > zeroTolerance) {
1691        unsigned char iStatus = internalStatus[iSequence] & 7;
1692        if (iStatus < 4) {
1693          index[numberNonZero] = iSequence;
1694          array[numberNonZero++] = tableauValue;
1695          double mult = multiplier[iStatus];
1696          double alpha = tableauValue * mult;
1697          double oldValue = abcDj[iSequence] * mult;
1698          double value = oldValue - tentativeTheta * alpha;
1699          if (value < dualT) {
1700            bestPossible = CoinMax(bestPossible, alpha);
1701            value = oldValue - upperTheta * alpha;
1702            if (value < dualT && alpha >= acceptablePivot) {
1703              upperTheta = (oldValue - dualT) / alpha;
1704            }
1705            // add to list
1706            arrayCandidate[numberRemaining] = alpha;
1707            indexCandidate[numberRemaining++] = iSequence;
1708          }
1709        }
1710      }
1711    }
1712  } else {
1713    double badFree = 0.0;
1714    double freeAlpha = model_->currentAcceptablePivot();
1715    int freeSequenceIn = model_->freeSequenceIn();
1716    double currentDualTolerance = model_->currentDualTolerance();
1717    for (CoinBigIndex j = start; j < end; j++) {
1718      int iSequence = column[j];
1719      double value = element[j];
1720      double tableauValue = piValue * value;
1721      if (fabs(tableauValue) > zeroTolerance) {
1722        unsigned char iStatus = internalStatus[iSequence] & 7;
1723        if (iStatus < 6) {
1724          if ((iStatus & 4) == 0) {
1725            index[numberNonZero] = iSequence;
1726            array[numberNonZero++] = tableauValue;
1727            double mult = multiplier[iStatus];
1728            double alpha = tableauValue * mult;
1729            double oldValue = abcDj[iSequence] * mult;
1730            double value = oldValue - tentativeTheta * alpha;
1731            if (value < dualT) {
1732              bestPossible = CoinMax(bestPossible, alpha);
1733              value = oldValue - upperTheta * alpha;
1734              if (value < dualT && alpha >= acceptablePivot) {
1735                upperTheta = (oldValue - dualT) / alpha;
1736              }
1737              // add to list
1738              arrayCandidate[numberRemaining] = alpha;
1739              indexCandidate[numberRemaining++] = iSequence;
1740            }
1741          } else {
1742            bool keep;
1743            index[numberNonZero] = iSequence;
1744            array[numberNonZero++] = tableauValue;
1745            bestPossible = CoinMax(bestPossible, fabs(tableauValue));
1746            double oldValue = abcDj[iSequence];
1747            // If free has to be very large - should come in via dualRow
1748            //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
1749            //break;
1750            if (oldValue > currentDualTolerance) {
1751              keep = true;
1752            } else if (oldValue < -currentDualTolerance) {
1753              keep = true;
1754            } else {
1755              if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
1756                keep = true;
1757              } else {
1758                keep = false;
1759                badFree = CoinMax(badFree, fabs(tableauValue));
1760              }
1761            }
1762            if (keep) {
1763#ifdef PAN
1764              if (model_->fakeSuperBasic(iSequence) >= 0) {
1765#endif
1766                if (iSequence == freeSequenceIn)
1767                  tableauValue = COIN_DBL_MAX;
1768                // free - choose largest
1769                if (fabs(tableauValue) > fabs(freeAlpha)) {
1770                  freeAlpha = tableauValue;
1771                  freeSequence = iSequence;
1772                }
1773#ifdef PAN
1774              }
1775#endif
1776            }
1777          }
1778        }
1779      }
1780    }
1781  }
1782  tableauRow.setNumElementsPartition(0, numberNonZero);
1783  candidateList.setNumElementsPartition(0, numberRemaining);
1784  return upperTheta;
1785}
1786//#define PARALLEL2
1787#ifdef PARALLEL2
1788#undef cilk_for
1789#undef cilk_spawn
1790#undef cilk_sync
1791#include <cilk/cilk.h>
1792#include <cilk/cilk_api.h>
1793#endif
1794#if 0
1795static void compact(int numberBlocks,CoinIndexedVector * vector,const int * starts,const int * lengths)
1796{
1797  int numberNonZeroIn=vector->getNumElements();
1798  int * index = vector->getIndices();
1799  double * array = vector->denseVector();
1800  CoinAbcCompact(numberBlocks,numberNonZeroIn,
1801                 array,starts,lengths);
1802  int numberNonZero = CoinAbcCompact(numberBlocks,numberNonZeroIn,
1803                                     index,starts,lengths);
1804  vector->setNumElements(numberNonZero);
1805}
1806static void compactBoth(int numberBlocks,CoinIndexedVector * vector1,CoinIndexedVector * vector2,
1807                        const int * starts,const int * lengths1, const int * lengths2)
1808{
1809  cilk_spawn compact(numberBlocks,vector1,starts,lengths1);
1810  compact(numberBlocks,vector2,starts,lengths2);
1811  cilk_sync;
1812}
1813#endif
1814void AbcMatrix::rebalance() const
1815{
1816  int maximumRows = model_->maximumAbcNumberRows();
1817  int numberTotal = model_->numberTotal();
1818  /* rebalance
1819     For non-vector version
1820     each basic etc column is 1
1821     each real column is 5+2*nel
1822     each basic slack is 0
1823     each real slack is 3
1824   */
1825#if ABC_PARALLEL
1826  int howOften = CoinMax(model_->factorization()->maximumPivots(), 200);
1827  if ((model_->numberIterations() % howOften) == 0 || !startColumnBlock_[1]) {
1828    int numberCpus = model_->numberCpus();
1829    if (numberCpus > 1) {
1830      const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
1831      const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
1832      int numberRows = model_->numberRows();
1833      int total = 0;
1834      for (int iSequence = 0; iSequence < numberRows; iSequence++) {
1835        unsigned char iStatus = internalStatus[iSequence] & 7;
1836        if (iStatus < 4)
1837          total += 3;
1838      }
1839      int totalSlacks = total;
1840      CoinBigIndex end = columnStart[maximumRows];
1841      for (int iSequence = maximumRows; iSequence < numberTotal; iSequence++) {
1842        CoinBigIndex start = end;
1843        end = columnStart[iSequence + 1];
1844        unsigned char iStatus = internalStatus[iSequence] & 7;
1845        if (iStatus < 4)
1846          total += 5 + 2 * (end - start);
1847        else
1848          total += 1;
1849      }
1850      int chunk = (total + numberCpus - 1) / numberCpus;
1851      total = totalSlacks;
1852      int iCpu = 0;
1853      startColumnBlock_[0] = 0;
1854      end = columnStart[maximumRows];
1855      for (int iSequence = maximumRows; iSequence < numberTotal; iSequence++) {
1856        CoinBigIndex start = end;
1857        end = columnStart[iSequence + 1];
1858        unsigned char iStatus = internalStatus[iSequence] & 7;
1859        if (iStatus < 4)
1860          total += 5 + 2 * (end - start);
1861        else
1862          total += 1;
1863        if (total > chunk) {
1864          iCpu++;
1865          total = 0;
1866          startColumnBlock_[iCpu] = iSequence;
1867        }
1868      }
1869      assert(iCpu < numberCpus);
1870      iCpu++;
1871      for (; iCpu <= numberCpus; iCpu++)
1872        startColumnBlock_[iCpu] = numberTotal;
1873      numberColumnBlocks_ = numberCpus;
1874    } else {
1875      numberColumnBlocks_ = 1;
1876      startColumnBlock_[0] = 0;
1877      startColumnBlock_[1] = numberTotal;
1878    }
1879  }
1880#else
1881  startColumnBlock_[0] = 0;
1882  startColumnBlock_[1] = numberTotal;
1883#endif
1884}
1885double
1886AbcMatrix::dualColumn1(const CoinIndexedVector &update,
1887  CoinPartitionedVector &tableauRow, CoinPartitionedVector &candidateList) const
1888{
1889  int numberTotal = model_->numberTotal();
1890  // rebalance
1891  rebalance();
1892  tableauRow.setPackedMode(true);
1893  candidateList.setPackedMode(true);
1894  int number = update.getNumElements();
1895  double upperTheta;
1896  if (rowStart_ && number < 3) {
1897#if ABC_INSTRUMENT > 1
1898    {
1899      int n = update.getNumElements();
1900      abcPricing[n < 19 ? n : 19]++;
1901    }
1902#endif
1903    // always do serially
1904    // do slacks first
1905    int starts[2];
1906    starts[0] = 0;
1907    starts[1] = numberTotal;
1908    tableauRow.setPartitions(1, starts);
1909    candidateList.setPartitions(1, starts);
1910    upperTheta = static_cast< AbcSimplexDual * >(model_)->dualColumn1A();
1911    //assert (upperTheta>-100*model_->dualTolerance()||model_->sequenceIn()>=0||model_->lastFirstFree()>=0);
1912    int freeSequence = -1;
1913    // worth using row copy
1914    assert(number);
1915    if (number == 2) {
1916      upperTheta = dualColumn1Row2(upperTheta, freeSequence, update, tableauRow, candidateList);
1917    } else {
1918      upperTheta = dualColumn1Row1(upperTheta, freeSequence, update, tableauRow, candidateList);
1919    }
1920    if (freeSequence >= 0) {
1921      int numberNonZero = tableauRow.getNumElements(0);
1922      const int *COIN_RESTRICT index = tableauRow.getIndices();
1923      const double *COIN_RESTRICT array = tableauRow.denseVector();
1924      // search for free coming in
1925      double freeAlpha = 0.0;
1926      int bestSequence = model_->sequenceIn();
1927      if (bestSequence >= 0)
1928        freeAlpha = model_->alpha();
1929      index = tableauRow.getIndices();
1930      array = tableauRow.denseVector();
1931      // free variable - search
1932      for (int k = 0; k < numberNonZero; k++) {
1933        if (freeSequence == index[k]) {
1934          if (fabs(freeAlpha) < fabs(array[k])) {
1935            freeAlpha = array[k];
1936          }
1937          break;
1938        }
1939      }
1940      if (model_->sequenceIn() < 0 || fabs(freeAlpha) > fabs(model_->alpha())) {
1941        double oldValue = model_->djRegion()[freeSequence];
1942        model_->setSequenceIn(freeSequence);
1943        model_->setAlpha(freeAlpha);
1944        model_->setTheta(oldValue / freeAlpha);
1945      }
1946    }
1947  } else {
1948    // three or more
1949    // need to do better job on dividing up (but wait until vector or by row)
1950    upperTheta = parallelDual4(static_cast< AbcSimplexDual * >(model_));
1951  }
1952  //tableauRow.compact();
1953  //candidateList.compact();
1954#if 0 //ndef NDEBUG
1955  model_->checkArrays();
1956#endif
1957  candidateList.computeNumberElements();
1958  int numberRemaining = candidateList.getNumElements();
1959  if (!numberRemaining && model_->sequenceIn() < 0) {
1960    return COIN_DBL_MAX; // Looks infeasible
1961  } else {
1962    return upperTheta;
1963  }
1964}
1965#define _mm256_broadcast_sd(x) static_cast< __m256d >(__builtin_ia32_vbroadcastsd256(x))
1966#define _mm256_load_pd(x) *(__m256d *)(x)
1967#define _mm256_store_pd (s, x) * ((__m256d *)s) = x
1968void AbcMatrix::dualColumn1Part(int iBlock, int &sequenceIn, double &upperThetaResult,
1969  const CoinIndexedVector &update,
1970  CoinPartitionedVector &tableauRow, CoinPartitionedVector &candidateList) const
1971{
1972  double upperTheta = upperThetaResult;
1973#if 0
1974  double time0=CoinCpuTime();
1975#endif
1976  int maximumRows = model_->maximumAbcNumberRows();
1977  int firstIn = startColumnBlock_[iBlock];
1978  int last = startColumnBlock_[iBlock + 1];
1979  int numberNonZero;
1980  int numberRemaining;
1981  int first;
1982  if (firstIn == 0) {
1983    upperTheta = static_cast< AbcSimplexDual * >(model_)->dualColumn1A();
1984    numberNonZero = tableauRow.getNumElements(0);
1985    numberRemaining = candidateList.getNumElements(0);
1986    first = maximumRows;
1987  } else {
1988    first = firstIn;
1989    numberNonZero = first;
1990    numberRemaining = first;
1991  }
1992  sequenceIn = -1;
1993  // get matrix data pointers
1994  const int *COIN_RESTRICT row = matrix_->getIndices();
1995  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
1996  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
1997  const double *COIN_RESTRICT pi = update.denseVector();
1998  //const int *  COIN_RESTRICT piIndex = update.getIndices();
1999  int *COIN_RESTRICT index = tableauRow.getIndices();
2000  double *COIN_RESTRICT array = tableauRow.denseVector();
2001  // pivot elements
2002  double *COIN_RESTRICT arrayCandidate = candidateList.denseVector();
2003  // indices
2004  int *COIN_RESTRICT indexCandidate = candidateList.getIndices();
2005  const double *COIN_RESTRICT abcDj = model_->djRegion();
2006  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
2007  // do first pass to get possibles
2008  double bestPossible = 0.0;
2009  // We can also see if infeasible or pivoting on free
2010  double tentativeTheta = 1.0e25; // try with this much smaller as guess
2011  double acceptablePivot = model_->currentAcceptablePivot();
2012  double dualT = -model_->currentDualTolerance();
2013  const double multiplier[] = { 1.0, -1.0 };
2014  double zeroTolerance = model_->zeroTolerance();
2015  if (model_->ordinaryVariables()) {
2016#ifdef COUNT_COPY
2017    if (first > maximumRows || last < model_->numberTotal() || false) {
2018#endif
2019#if 1
2020      for (int iSequence = first; iSequence < last; iSequence++) {
2021        unsigned char iStatus = internalStatus[iSequence] & 7;
2022        if (iStatus < 4) {
2023          CoinBigIndex start = columnStart[iSequence];
2024          CoinBigIndex next = columnStart[iSequence + 1];
2025          double tableauValue = 0.0;
2026          for (CoinBigIndex j = start; j < next; j++) {
2027            int jRow = row[j];
2028            tableauValue += pi[jRow] * elementByColumn[j];
2029          }
2030          if (fabs(tableauValue) > zeroTolerance) {
2031#else
2032    cilk_for(int iSequence = first; iSequence < last; iSequence++)
2033    {
2034      unsigned char iStatus = internalStatus[iSequence] & 7;
2035      if (iStatus < 4) {
2036        CoinBigIndex start = columnStart[iSequence];
2037        CoinBigIndex next = columnStart[iSequence + 1];
2038        double tableauValue = 0.0;
2039        for (CoinBigIndex j = start; j < next; j++) {
2040          int jRow = row[j];
2041          tableauValue += pi[jRow] * elementByColumn[j];
2042        }
2043        array[iSequence] = tableauValue;
2044      }
2045    }
2046    //printf("time %.6g\n",CoinCpuTime()-time0);
2047    for (int iSequence = first; iSequence < last; iSequence++) {
2048      double tableauValue = array[iSequence];
2049      if (tableauValue) {
2050        array[iSequence] = 0.0;
2051        if (fabs(tableauValue) > zeroTolerance) {
2052          unsigned char iStatus = internalStatus[iSequence] & 7;
2053#endif
2054            index[numberNonZero] = iSequence;
2055            array[numberNonZero++] = tableauValue;
2056            double mult = multiplier[iStatus];
2057            double alpha = tableauValue * mult;
2058            double oldValue = abcDj[iSequence] * mult;
2059            double value = oldValue - tentativeTheta * alpha;
2060            if (value < dualT) {
2061              bestPossible = CoinMax(bestPossible, alpha);
2062              value = oldValue - upperTheta * alpha;
2063              if (value < dualT && alpha >= acceptablePivot) {
2064                upperTheta = (oldValue - dualT) / alpha;
2065              }
2066              // add to list
2067              arrayCandidate[numberRemaining] = alpha;
2068              indexCandidate[numberRemaining++] = iSequence;
2069            }
2070          }
2071        }
2072      }
2073#ifdef COUNT_COPY
2074    } else {
2075      const double *COIN_RESTRICT element = countElement_;
2076      const int *COIN_RESTRICT row = countRow_;
2077      double temp[4] __attribute__((aligned(32)));
2078      memset(temp, 0, sizeof(temp));
2079      for (int iCount = smallestCount_; iCount < largestCount_; iCount++) {
2080        int iStart = countFirst_[iCount];
2081        int number = countFirst_[iCount + 1] - iStart;
2082        if (!number)
2083          continue;
2084        const int *COIN_RESTRICT blockRow = row + countStart_[iCount];
2085        const double *COIN_RESTRICT blockElement = element + countStart_[iCount];
2086        const int *COIN_RESTRICT sequence = countRealColumn_ + iStart;
2087        // really need to sort and swap to avoid tests
2088        int numberBlocks = number >> 2;
2089        number &= 3;
2090        for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
2091          double tableauValue0 = 0.0;
2092          double tableauValue1 = 0.0;
2093          double tableauValue2 = 0.0;
2094          double tableauValue3 = 0.0;
2095          for (int j = 0; j < iCount; j++) {
2096            tableauValue0 += pi[blockRow[0]] * blockElement[0];
2097            tableauValue1 += pi[blockRow[1]] * blockElement[1];
2098            tableauValue2 += pi[blockRow[2]] * blockElement[2];
2099            tableauValue3 += pi[blockRow[3]] * blockElement[3];
2100            blockRow += 4;
2101            blockElement += 4;
2102          }
2103          if (fabs(tableauValue0) > zeroTolerance) {
2104            int iSequence = sequence[0];
2105            unsigned char iStatus = internalStatus[iSequence] & 7;
2106            if (iStatus < 4) {
2107              index[numberNonZero] = iSequence;
2108              array[numberNonZero++] = tableauValue0;
2109              double mult = multiplier[iStatus];
2110              double alpha = tableauValue0 * mult;
2111              double oldValue = abcDj[iSequence] * mult;
2112              double value = oldValue - tentativeTheta * alpha;
2113              if (value < dualT) {
2114                bestPossible = CoinMax(bestPossible, alpha);
2115                value = oldValue - upperTheta * alpha;
2116                if (value < dualT && alpha >= acceptablePivot) {
2117                  upperTheta = (oldValue - dualT) / alpha;
2118                }
2119                // add to list
2120                arrayCandidate[numberRemaining] = alpha;
2121                indexCandidate[numberRemaining++] = iSequence;
2122              }
2123            }
2124          }
2125          if (fabs(tableauValue1) > zeroTolerance) {
2126            int iSequence = sequence[1];
2127            unsigned char iStatus = internalStatus[iSequence] & 7;
2128            if (iStatus < 4) {
2129              index[numberNonZero] = iSequence;
2130              array[numberNonZero++] = tableauValue1;
2131              double mult = multiplier[iStatus];
2132              double alpha = tableauValue1 * mult;
2133              double oldValue = abcDj[iSequence] * mult;
2134              double value = oldValue - tentativeTheta * alpha;
2135              if (value < dualT) {
2136                bestPossible = CoinMax(bestPossible, alpha);
2137                value = oldValue - upperTheta * alpha;
2138                if (value < dualT && alpha >= acceptablePivot) {
2139                  upperTheta = (oldValue - dualT) / alpha;
2140                }
2141                // add to list
2142                arrayCandidate[numberRemaining] = alpha;
2143                indexCandidate[numberRemaining++] = iSequence;
2144              }
2145            }
2146          }
2147          if (fabs(tableauValue2) > zeroTolerance) {
2148            int iSequence = sequence[2];
2149            unsigned char iStatus = internalStatus[iSequence] & 7;
2150            if (iStatus < 4) {
2151              index[numberNonZero] = iSequence;
2152              array[numberNonZero++] = tableauValue2;
2153              double mult = multiplier[iStatus];
2154              double alpha = tableauValue2 * mult;
2155              double oldValue = abcDj[iSequence] * mult;
2156              double value = oldValue - tentativeTheta * alpha;
2157              if (value < dualT) {
2158                bestPossible = CoinMax(bestPossible, alpha);
2159                value = oldValue - upperTheta * alpha;
2160                if (value < dualT && alpha >= acceptablePivot) {
2161                  upperTheta = (oldValue - dualT) / alpha;
2162                }
2163                // add to list
2164                arrayCandidate[numberRemaining] = alpha;
2165                indexCandidate[numberRemaining++] = iSequence;
2166              }
2167            }
2168          }
2169          if (fabs(tableauValue3) > zeroTolerance) {
2170            int iSequence = sequence[3];
2171            unsigned char iStatus = internalStatus[iSequence] & 7;
2172            if (iStatus < 4) {
2173              index[numberNonZero] = iSequence;
2174              array[numberNonZero++] = tableauValue3;
2175              double mult = multiplier[iStatus];
2176              double alpha = tableauValue3 * mult;
2177              double oldValue = abcDj[iSequence] * mult;
2178              double value = oldValue - tentativeTheta * alpha;
2179              if (value < dualT) {
2180                bestPossible = CoinMax(bestPossible, alpha);
2181                value = oldValue - upperTheta * alpha;
2182                if (value < dualT && alpha >= acceptablePivot) {
2183                  upperTheta = (oldValue - dualT) / alpha;
2184                }
2185                // add to list
2186                arrayCandidate[numberRemaining] = alpha;
2187                indexCandidate[numberRemaining++] = iSequence;
2188              }
2189            }
2190          }
2191          sequence += 4;
2192        }
2193        for (int i = 0; i < number; i++) {
2194          int iSequence = sequence[i];
2195          unsigned char iStatus = internalStatus[iSequence] & 7;
2196          if (iStatus < 4) {
2197            double tableauValue = 0.0;
2198            for (int j = 0; j < iCount; j++) {
2199              int iRow = blockRow[4 * j];
2200              tableauValue += pi[iRow] * blockElement[4 * j];
2201            }
2202            if (fabs(tableauValue) > zeroTolerance) {
2203              index[numberNonZero] = iSequence;
2204              array[numberNonZero++] = tableauValue;
2205              double mult = multiplier[iStatus];
2206              double alpha = tableauValue * mult;
2207              double oldValue = abcDj[iSequence] * mult;
2208              double value = oldValue - tentativeTheta * alpha;
2209              if (value < dualT) {
2210                bestPossible = CoinMax(bestPossible, alpha);
2211                value = oldValue - upperTheta * alpha;
2212                if (value < dualT && alpha >= acceptablePivot) {
2213                  upperTheta = (oldValue - dualT) / alpha;
2214                }
2215                // add to list
2216                arrayCandidate[numberRemaining] = alpha;
2217                indexCandidate[numberRemaining++] = iSequence;
2218              }
2219            }
2220          }
2221          blockRow++;
2222          blockElement++;
2223        }
2224      }
2225      int numberColumns = model_->numberColumns();
2226      // large ones
2227      const CoinBigIndex *COIN_RESTRICT largeStart = countStartLarge_ - countFirst_[MAX_COUNT];
2228      for (int i = countFirst_[MAX_COUNT]; i < numberColumns; i++) {
2229        int iSequence = countRealColumn_[i];
2230        unsigned char iStatus = internalStatus[iSequence] & 7;
2231        if (iStatus < 4) {
2232          CoinBigIndex start = largeStart[i];
2233          CoinBigIndex next = largeStart[i + 1];
2234          double tableauValue = 0.0;
2235          for (CoinBigIndex j = start; j < next; j++) {
2236            int jRow = row[j];
2237            tableauValue += pi[jRow] * element[j];
2238          }
2239          if (fabs(tableauValue) > zeroTolerance) {
2240            index[numberNonZero] = iSequence;
2241            array[numberNonZero++] = tableauValue;
2242            double mult = multiplier[iStatus];
2243            double alpha = tableauValue * mult;
2244            double oldValue = abcDj[iSequence] * mult;
2245            double value = oldValue - tentativeTheta * alpha;
2246            if (value < dualT) {
2247              bestPossible = CoinMax(bestPossible, alpha);
2248              value = oldValue - upperTheta * alpha;
2249              if (value < dualT && alpha >= acceptablePivot) {
2250                upperTheta = (oldValue - dualT) / alpha;
2251              }
2252              // add to list
2253              arrayCandidate[numberRemaining] = alpha;
2254              indexCandidate[numberRemaining++] = iSequence;
2255            }
2256          }
2257        }
2258      }
2259    }
2260#endif
2261  } else {
2262    double badFree = 0.0;
2263    double freePivot = model_->currentAcceptablePivot();
2264    int freeSequenceIn = model_->freeSequenceIn();
2265    double currentDualTolerance = model_->currentDualTolerance();
2266    for (int iSequence = first; iSequence < last; iSequence++) {
2267      unsigned char iStatus = internalStatus[iSequence] & 7;
2268      if (iStatus < 6) {
2269        CoinBigIndex start = columnStart[iSequence];
2270        CoinBigIndex next = columnStart[iSequence + 1];
2271        double tableauValue = 0.0;
2272        for (CoinBigIndex j = start; j < next; j++) {
2273          int jRow = row[j];
2274          tableauValue += pi[jRow] * elementByColumn[j];
2275        }
2276        if (fabs(tableauValue) > zeroTolerance) {
2277          if ((iStatus & 4) == 0) {
2278            index[numberNonZero] = iSequence;
2279            array[numberNonZero++] = tableauValue;
2280            double mult = multiplier[iStatus];
2281            double alpha = tableauValue * mult;
2282            double oldValue = abcDj[iSequence] * mult;
2283            double value = oldValue - tentativeTheta * alpha;
2284            if (value < dualT) {
2285              bestPossible = CoinMax(bestPossible, alpha);
2286              value = oldValue - upperTheta * alpha;
2287              if (value < dualT && alpha >= acceptablePivot) {
2288                upperTheta = (oldValue - dualT) / alpha;
2289              }
2290              // add to list
2291              arrayCandidate[numberRemaining] = alpha;
2292              indexCandidate[numberRemaining++] = iSequence;
2293            }
2294          } else {
2295            bool keep;
2296            index[numberNonZero] = iSequence;
2297            array[numberNonZero++] = tableauValue;
2298            bestPossible = CoinMax(bestPossible, fabs(tableauValue));
2299            double oldValue = abcDj[iSequence];
2300            // If free has to be very large - should come in via dualRow
2301            //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
2302            //break;
2303            if (oldValue > currentDualTolerance) {
2304              keep = true;
2305            } else if (oldValue < -currentDualTolerance) {
2306              keep = true;
2307            } else {
2308              if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
2309                keep = true;
2310              } else {
2311                keep = false;
2312                badFree = CoinMax(badFree, fabs(tableauValue));
2313              }
2314            }
2315            if (keep) {
2316#ifdef PAN
2317              if (model_->fakeSuperBasic(iSequence) >= 0) {
2318#endif
2319                if (iSequence == freeSequenceIn)
2320                  tableauValue = COIN_DBL_MAX;
2321                // free - choose largest
2322                if (fabs(tableauValue) > fabs(freePivot)) {
2323                  freePivot = tableauValue;
2324                  sequenceIn = iSequence;
2325                }
2326#ifdef PAN
2327              }
2328#endif
2329            }
2330          }
2331        }
2332      }
2333    }
2334  }
2335  // adjust
2336  numberNonZero -= firstIn;
2337  numberRemaining -= firstIn;
2338  tableauRow.setNumElementsPartition(iBlock, numberNonZero);
2339  candidateList.setNumElementsPartition(iBlock, numberRemaining);
2340  if (!numberRemaining && model_->sequenceIn() < 0) {
2341
2342    upperThetaResult = COIN_DBL_MAX; // Looks infeasible
2343  } else {
2344    upperThetaResult = upperTheta;
2345  }
2346}
2347// gets tableau row - returns number of slacks in block
2348int AbcMatrix::primalColumnRow(int iBlock, bool doByRow, const CoinIndexedVector &updates,
2349  CoinPartitionedVector &tableauRow) const
2350{
2351  int maximumRows = model_->maximumAbcNumberRows();
2352  int first = tableauRow.startPartition(iBlock);
2353  int last = tableauRow.startPartition(iBlock + 1);
2354  if (doByRow) {
2355    assert(first == blockStart_[iBlock]);
2356    assert(last == blockStart_[iBlock + 1]);
2357  } else {
2358    assert(first == startColumnBlock_[iBlock]);
2359    assert(last == startColumnBlock_[iBlock + 1]);
2360  }
2361  int numberSlacks = updates.getNumElements();
2362  int numberNonZero = 0;
2363  const double *COIN_RESTRICT pi = updates.denseVector();
2364  const int *COIN_RESTRICT piIndex = updates.getIndices();
2365  int *COIN_RESTRICT index = tableauRow.getIndices() + first;
2366  double *COIN_RESTRICT array = tableauRow.denseVector() + first;
2367  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
2368  if (!iBlock) {
2369    numberNonZero = numberSlacks;
2370    for (int i = 0; i < numberSlacks; i++) {
2371      int iRow = piIndex[i];
2372      index[i] = iRow;
2373      array[i] = pi[iRow]; // ? what if small or basic
2374    }
2375    first = maximumRows;
2376  }
2377  double zeroTolerance = model_->zeroTolerance();
2378  if (doByRow) {
2379    int numberRows = model_->numberRows();
2380    const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_ + iBlock * numberRows;
2381    const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows;
2382    const double *COIN_RESTRICT element = element_;
2383    const int *COIN_RESTRICT column = column_;
2384    if (numberSlacks > 1) {
2385      double *COIN_RESTRICT arrayTemp = tableauRow.denseVector();
2386      for (int i = 0; i < numberSlacks; i++) {
2387        int iRow = piIndex[i];
2388        double piValue = pi[iRow];
2389        CoinBigIndex end = rowEnd[iRow];
2390        for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
2391          int iColumn = column[j];
2392          arrayTemp[iColumn] += element[j] * piValue;
2393        }
2394      }
2395      for (int iSequence = first; iSequence < last; iSequence++) {
2396        double tableauValue = arrayTemp[iSequence];
2397        if (tableauValue) {
2398          arrayTemp[iSequence] = 0.0;
2399          if (fabs(tableauValue) > zeroTolerance) {
2400            unsigned char iStatus = internalStatus[iSequence] & 7;
2401            if (iStatus < 6) {
2402              index[numberNonZero] = iSequence;
2403              array[numberNonZero++] = tableauValue;
2404            }
2405          }
2406        }
2407      }
2408    } else {
2409      int iRow = piIndex[0];
2410      double piValue = pi[iRow];
2411      CoinBigIndex end = rowEnd[iRow];
2412      for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
2413        int iSequence = column[j];
2414        double tableauValue = element[j] * piValue;
2415        if (fabs(tableauValue) > zeroTolerance) {
2416          unsigned char iStatus = internalStatus[iSequence] & 7;
2417          if (iStatus < 6) {
2418            index[numberNonZero] = iSequence;
2419            array[numberNonZero++] = tableauValue;
2420          }
2421        }
2422      }
2423    }
2424  } else {
2425    // get matrix data pointers
2426    const int *COIN_RESTRICT row = matrix_->getIndices();
2427    const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
2428    const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
2429    const double *COIN_RESTRICT pi = updates.denseVector();
2430    for (int iSequence = first; iSequence < last; iSequence++) {
2431      unsigned char iStatus = internalStatus[iSequence] & 7;
2432      if (iStatus < 6) {
2433        CoinBigIndex start = columnStart[iSequence];
2434        CoinBigIndex next = columnStart[iSequence + 1];
2435        double tableauValue = 0.0;
2436        for (CoinBigIndex j = start; j < next; j++) {
2437          int jRow = row[j];
2438          tableauValue += pi[jRow] * elementByColumn[j];
2439        }
2440        if (fabs(tableauValue) > zeroTolerance) {
2441          index[numberNonZero] = iSequence;
2442          array[numberNonZero++] = tableauValue;
2443        }
2444      }
2445    }
2446  }
2447  tableauRow.setNumElementsPartition(iBlock, numberNonZero);
2448  return numberSlacks;
2449}
2450// Get sequenceIn when Dantzig
2451int AbcMatrix::pivotColumnDantzig(const CoinIndexedVector &updates,
2452  CoinPartitionedVector &spare) const
2453{
2454  int maximumRows = model_->maximumAbcNumberRows();
2455  // rebalance
2456  rebalance();
2457  spare.setPackedMode(true);
2458  bool useRowCopy = false;
2459  if (rowStart_) {
2460    // see if worth using row copy
2461    int number = updates.getNumElements();
2462    assert(number);
2463    useRowCopy = (number << 2) < maximumRows;
2464  }
2465  const int *starts;
2466  if (useRowCopy)
2467    starts = blockStart();
2468  else
2469    starts = startColumnBlock();
2470#if ABC_PARALLEL
2471#define NUMBER_BLOCKS NUMBER_ROW_BLOCKS
2472  int numberBlocks = CoinMin(NUMBER_BLOCKS, model_->numberCpus());
2473#else
2474#define NUMBER_BLOCKS 1
2475  int numberBlocks = NUMBER_BLOCKS;
2476#endif
2477#if ABC_PARALLEL
2478  if (model_->parallelMode() == 0)
2479    numberBlocks = 1;
2480#endif
2481  spare.setPartitions(numberBlocks, starts);
2482  int which[NUMBER_BLOCKS];
2483  double best[NUMBER_BLOCKS];
2484  for (int i = 0; i < numberBlocks - 1; i++)
2485    which[i] = cilk_spawn pivotColumnDantzig(i, useRowCopy, updates, spare, best[i]);
2486  which[numberBlocks - 1] = pivotColumnDantzig(numberBlocks - 1, useRowCopy, updates,
2487    spare, best[numberBlocks - 1]);
2488  cilk_sync;
2489  int bestSequence = -1;
2490  double bestValue = model_->dualTolerance();
2491  for (int i = 0; i < numberBlocks; i++) {
2492    if (best[i] > bestValue) {
2493      bestValue = best[i];
2494      bestSequence = which[i];
2495    }
2496  }
2497  return bestSequence;
2498}
2499// Get sequenceIn when Dantzig (One block)
2500int AbcMatrix::pivotColumnDantzig(int iBlock, bool doByRow, const CoinIndexedVector &updates,
2501  CoinPartitionedVector &spare,
2502  double &bestValue) const
2503{
2504  // could rewrite to do more directly
2505  int numberSlacks = primalColumnRow(iBlock, doByRow, updates, spare);
2506  double *COIN_RESTRICT reducedCost = model_->djRegion();
2507  int first = spare.startPartition(iBlock);
2508  int last = spare.startPartition(iBlock + 1);
2509  int *COIN_RESTRICT index = spare.getIndices() + first;
2510  double *COIN_RESTRICT array = spare.denseVector() + first;
2511  int numberNonZero = spare.getNumElements(iBlock);
2512  int bestSequence = -1;
2513  double bestDj = model_->dualTolerance();
2514  double bestFreeDj = model_->dualTolerance();
2515  int bestFreeSequence = -1;
2516  // redo LOTS so signs for slacks and columns same
2517  if (!first) {
2518    first = model_->maximumAbcNumberRows();
2519    for (int i = 0; i < numberSlacks; i++) {
2520      int iSequence = index[i];
2521      double value = reducedCost[iSequence];
2522      value += array[i];
2523      array[i] = 0.0;
2524      reducedCost[iSequence] = value;
2525    }
2526#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
2527#define CLP_PRIMAL_SLACK_MULTIPLIER 1.0
2528#endif
2529    int numberRows = model_->numberRows();
2530    for (int iSequence = 0; iSequence < numberRows; iSequence++) {
2531      // check flagged variable
2532      if (!model_->flagged(iSequence)) {
2533        double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER;
2534        AbcSimplex::Status status = model_->getInternalStatus(iSequence);
2535
2536        switch (status) {
2537
2538        case AbcSimplex::basic:
2539        case AbcSimplex::isFixed:
2540          break;
2541        case AbcSimplex::isFree:
2542        case AbcSimplex::superBasic:
2543          if (fabs(value) > bestFreeDj) {
2544            bestFreeDj = fabs(value);
2545            bestFreeSequence = iSequence;
2546          }
2547          break;
2548        case AbcSimplex::atUpperBound:
2549          if (value > bestDj) {
2550            bestDj = value;
2551            bestSequence = iSequence;
2552          }
2553          break;
2554        case AbcSimplex::atLowerBound:
2555          if (value < -bestDj) {
2556            bestDj = -value;
2557            bestSequence = iSequence;
2558          }
2559        }
2560      }
2561    }
2562  }
2563  for (int i = numberSlacks; i < numberNonZero; i++) {
2564    int iSequence = index[i];
2565    double value = reducedCost[iSequence];
2566    value += array[i];
2567    array[i] = 0.0;
2568    reducedCost[iSequence] = value;
2569  }
2570  for (int iSequence = first; iSequence < last; iSequence++) {
2571    // check flagged variable
2572    if (!model_->flagged(iSequence)) {
2573      double value = reducedCost[iSequence];
2574      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
2575
2576      switch (status) {
2577
2578      case AbcSimplex::basic:
2579      case AbcSimplex::isFixed:
2580        break;
2581      case AbcSimplex::isFree:
2582      case AbcSimplex::superBasic:
2583        if (fabs(value) > bestFreeDj) {
2584          bestFreeDj = fabs(value);
2585          bestFreeSequence = iSequence;
2586        }
2587        break;
2588      case AbcSimplex::atUpperBound:
2589        if (value > bestDj) {
2590          bestDj = value;
2591          bestSequence = iSequence;
2592        }
2593        break;
2594      case AbcSimplex::atLowerBound:
2595        if (value < -bestDj) {
2596          bestDj = -value;
2597          bestSequence = iSequence;
2598        }
2599      }
2600    }
2601  }
2602  spare.setNumElementsPartition(iBlock, 0);
2603  // bias towards free
2604  if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj) {
2605    bestSequence = bestFreeSequence;
2606    bestDj = 10.0 * bestFreeDj;
2607  }
2608  bestValue = bestDj;
2609  return bestSequence;
2610}
2611// gets tableau row and dj row - returns number of slacks in block
2612int AbcMatrix::primalColumnRowAndDjs(int iBlock, const CoinIndexedVector &updateTableau,
2613  const CoinIndexedVector &updateDjs,
2614  CoinPartitionedVector &tableauRow) const
2615{
2616  int maximumRows = model_->maximumAbcNumberRows();
2617  int first = tableauRow.startPartition(iBlock);
2618  int last = tableauRow.startPartition(iBlock + 1);
2619  assert(first == startColumnBlock_[iBlock]);
2620  assert(last == startColumnBlock_[iBlock + 1]);
2621  int numberSlacks = updateTableau.getNumElements();
2622  int numberNonZero = 0;
2623  const double *COIN_RESTRICT piTableau = updateTableau.denseVector();
2624  const double *COIN_RESTRICT pi = updateDjs.denseVector();
2625  int *COIN_RESTRICT index = tableauRow.getIndices() + first;
2626  double *COIN_RESTRICT array = tableauRow.denseVector() + first;
2627  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
2628  double *COIN_RESTRICT reducedCost = model_->djRegion();
2629  if (!iBlock) {
2630    const int *COIN_RESTRICT piIndexTableau = updateTableau.getIndices();
2631    numberNonZero = numberSlacks;
2632    for (int i = 0; i < numberSlacks; i++) {
2633      int iRow = piIndexTableau[i];
2634      index[i] = iRow;
2635      array[i] = piTableau[iRow]; // ? what if small or basic
2636    }
2637    const int *COIN_RESTRICT piIndex = updateDjs.getIndices();
2638    int number = updateDjs.getNumElements();
2639    for (int i = 0; i < number; i++) {
2640      int iRow = piIndex[i];
2641      reducedCost[iRow] -= pi[iRow]; // sign?
2642    }
2643    first = maximumRows;
2644  }
2645  double zeroTolerance = model_->zeroTolerance();
2646  // get matrix data pointers
2647  const int *COIN_RESTRICT row = matrix_->getIndices();
2648  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
2649  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
2650  for (int iSequence = first; iSequence < last; iSequence++) {
2651    unsigned char iStatus = internalStatus[iSequence] & 7;
2652    if (iStatus < 6) {
2653      CoinBigIndex start = columnStart[iSequence];
2654      CoinBigIndex next = columnStart[iSequence + 1];
2655      double tableauValue = 0.0;
2656      double djValue = reducedCost[iSequence];
2657      for (CoinBigIndex j = start; j < next; j++) {
2658        int jRow = row[j];
2659        tableauValue += piTableau[jRow] * elementByColumn[j];
2660        djValue -= pi[jRow] * elementByColumn[j]; // sign?
2661      }
2662      reducedCost[iSequence] = djValue;
2663      if (fabs(tableauValue) > zeroTolerance) {
2664        index[numberNonZero] = iSequence;
2665        array[numberNonZero++] = tableauValue;
2666      }
2667    }
2668  }
2669  tableauRow.setNumElementsPartition(iBlock, numberNonZero);
2670  return numberSlacks;
2671}
2672/* does steepest edge double or triple update
2673   If scaleFactor!=0 then use with tableau row to update djs
2674   otherwise use updateForDjs
2675*/
2676int AbcMatrix::primalColumnDouble(int iBlock, CoinPartitionedVector &updateForTableauRow,
2677  CoinPartitionedVector &updateForDjs,
2678  const CoinIndexedVector &updateForWeights,
2679  CoinPartitionedVector &spareColumn1,
2680  double *infeasibilities,
2681  double referenceIn, double devex,
2682  // Array for exact devex to say what is in reference framework
2683  unsigned int *reference,
2684  double *weights, double scaleFactor) const
2685{
2686  int maximumRows = model_->maximumAbcNumberRows();
2687  int first = startColumnBlock_[iBlock];
2688  int last = startColumnBlock_[iBlock + 1];
2689  const double *COIN_RESTRICT piTableau = updateForTableauRow.denseVector();
2690  const double *COIN_RESTRICT pi = updateForDjs.denseVector();
2691  const double *COIN_RESTRICT piWeights = updateForWeights.denseVector();
2692  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
2693  double *COIN_RESTRICT reducedCost = model_->djRegion();
2694  double tolerance = model_->currentDualTolerance();
2695  // use spareColumn to track new infeasibilities
2696  int *COIN_RESTRICT newInfeas = spareColumn1.getIndices() + first;
2697  int numberNewInfeas = 0;
2698  // we can't really trust infeasibilities if there is dual error
2699  // this coding has to mimic coding in checkDualSolution
2700  double error = CoinMin(1.0e-2, model_->largestDualError());
2701  // allow tolerance at least slightly bigger than standard
2702  tolerance = tolerance + error;
2703  double mult[2] = { 1.0, -1.0 };
2704  bool updateWeights = devex != 0.0;
2705#define isReference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0)
2706  // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor*
2707  if (!iBlock) {
2708    int numberSlacks = updateForTableauRow.getNumElements();
2709    const int *COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
2710    if (!scaleFactor) {
2711      const int *COIN_RESTRICT piIndex = updateForDjs.getIndices();
2712      int number = updateForDjs.getNumElements();
2713      for (int i = 0; i < number; i++) {
2714        int iRow = piIndex[i];
2715        int iStatus = internalStatus[iRow] & 7;
2716        double value = reducedCost[iRow];
2717        value += pi[iRow];
2718        if (iStatus < 4) {
2719          reducedCost[iRow] = value;
2720          value *= mult[iStatus];
2721          if (value < 0.0) {
2722            if (!infeasibilities[iRow])
2723              newInfeas[numberNewInfeas++] = iRow;
2724#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
2725            infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER;
2726#else
2727            infeasibilities[iRow] = value * value;
2728#endif
2729          } else {
2730            if (infeasibilities[iRow])
2731              infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
2732          }
2733        }
2734      }
2735    } else if (scaleFactor != COIN_DBL_MAX) {
2736      for (int i = 0; i < numberSlacks; i++) {
2737        int iRow = piIndexTableau[i];
2738        int iStatus = internalStatus[iRow] & 7;
2739        if (iStatus < 4) {
2740          double value = reducedCost[iRow];
2741          value += scaleFactor * piTableau[iRow]; // sign?
2742          reducedCost[iRow] = value;
2743          value *= mult[iStatus];
2744          if (value < 0.0) {
2745            if (!infeasibilities[iRow])
2746              newInfeas[numberNewInfeas++] = iRow;
2747#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
2748            infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER;
2749#else
2750            infeasibilities[iRow] = value * value;
2751#endif
2752          } else {
2753            if (infeasibilities[iRow])
2754              infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
2755          }
2756        }
2757      }
2758    }
2759    if (updateWeights) {
2760      for (int i = 0; i < numberSlacks; i++) {
2761        int iRow = piIndexTableau[i];
2762        double modification = piWeights[iRow];
2763        double thisWeight = weights[iRow];
2764        double pivot = piTableau[iRow];
2765        double pivotSquared = pivot * pivot;
2766        thisWeight += pivotSquared * devex - pivot * modification;
2767        if (thisWeight < DEVEX_TRY_NORM) {
2768          if (referenceIn < 0.0) {
2769            // steepest
2770            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2771          } else {
2772            // exact
2773            thisWeight = referenceIn * pivotSquared;
2774            if (isReference(iRow))
2775              thisWeight += 1.0;
2776            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2777          }
2778        }
2779        weights[iRow] = thisWeight;
2780      }
2781    }
2782    first = maximumRows;
2783  }
2784  double zeroTolerance = model_->zeroTolerance();
2785  double freeTolerance = FREE_ACCEPT * tolerance;
2786  ;
2787  // get matrix data pointers
2788  const int *COIN_RESTRICT row = matrix_->getIndices();
2789  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
2790  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
2791  bool updateDjs = scaleFactor != COIN_DBL_MAX;
2792  for (int iSequence = first; iSequence < last; iSequence++) {
2793    unsigned char iStatus = internalStatus[iSequence] & 7;
2794    if (iStatus < 6) {
2795      CoinBigIndex start = columnStart[iSequence];
2796      CoinBigIndex next = columnStart[iSequence + 1];
2797      double tableauValue = 0.0;
2798      double djValue = reducedCost[iSequence];
2799      if (!scaleFactor) {
2800        for (CoinBigIndex j = start; j < next; j++) {
2801          int jRow = row[j];
2802          tableauValue += piTableau[jRow] * elementByColumn[j];
2803          djValue += pi[jRow] * elementByColumn[j];
2804        }
2805      } else {
2806        for (CoinBigIndex j = start; j < next; j++) {
2807          int jRow = row[j];
2808          tableauValue += piTableau[jRow] * elementByColumn[j];
2809        }
2810        djValue += tableauValue * scaleFactor; // sign?
2811      }
2812      if (updateDjs) {
2813        reducedCost[iSequence] = djValue;
2814        if (iStatus < 4) {
2815          djValue *= mult[iStatus];
2816          if (djValue < 0.0) {
2817            if (!infeasibilities[iSequence])
2818              newInfeas[numberNewInfeas++] = iSequence;
2819            infeasibilities[iSequence] = djValue * djValue;
2820          } else {
2821            if (infeasibilities[iSequence])
2822              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
2823          }
2824        } else {
2825          if (fabs(djValue) > freeTolerance) {
2826            // we are going to bias towards free (but only if reasonable)
2827            djValue *= FREE_BIAS;
2828            if (!infeasibilities[iSequence])
2829              newInfeas[numberNewInfeas++] = iSequence;
2830            // store square in list
2831            infeasibilities[iSequence] = djValue * djValue;
2832          } else {
2833            if (infeasibilities[iSequence])
2834              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
2835          }
2836        }
2837      }
2838      if (fabs(tableauValue) > zeroTolerance && updateWeights) {
2839        double modification = 0.0;
2840        for (CoinBigIndex j = start; j < next; j++) {
2841          int jRow = row[j];
2842          modification += piWeights[jRow] * elementByColumn[j];
2843        }
2844        double thisWeight = weights[iSequence];
2845        double pivot = tableauValue;
2846        double pivotSquared = pivot * pivot;
2847        thisWeight += pivotSquared * devex - pivot * modification;
2848        if (thisWeight < DEVEX_TRY_NORM) {
2849          if (referenceIn < 0.0) {
2850            // steepest
2851            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2852          } else {
2853            // exact
2854            thisWeight = referenceIn * pivotSquared;
2855            if (isReference(iSequence))
2856              thisWeight += 1.0;
2857            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2858          }
2859        }
2860        weights[iSequence] = thisWeight;
2861      }
2862    }
2863  }
2864  spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas);
2865  int bestSequence = -1;
2866#if 0
2867  if (!iBlock)
2868    first=0;
2869  // not at present - maybe later
2870  double bestDj=tolerance*tolerance;
2871  for (int iSequence=first;iSequence<last;iSequence++) {
2872    if (infeasibilities[iSequence]>bestDj*weights[iSequence]) {
2873      int iStatus=internalStatus[iSequence]&7;
2874      assert (iStatus<6);
2875      bestSequence=iSequence;
2876      bestDj=infeasibilities[iSequence]/weights[iSequence];
2877    }
2878  }
2879#endif
2880  return bestSequence;
2881}
2882/* does steepest edge double or triple update
2883   If scaleFactor!=0 then use with tableau row to update djs
2884   otherwise use updateForDjs
2885*/
2886int AbcMatrix::primalColumnSparseDouble(int iBlock, CoinPartitionedVector &updateForTableauRow,
2887  CoinPartitionedVector &updateForDjs,
2888  const CoinIndexedVector &updateForWeights,
2889  CoinPartitionedVector &spareColumn1,
2890  double *infeasibilities,
2891  double referenceIn, double devex,
2892  // Array for exact devex to say what is in reference framework
2893  unsigned int *reference,
2894  double *weights, double scaleFactor) const
2895{
2896  int maximumRows = model_->maximumAbcNumberRows();
2897  int first = blockStart_[iBlock];
2898  //int last=blockStart_[iBlock+1];
2899  const double *COIN_RESTRICT piTableau = updateForTableauRow.denseVector();
2900  //const double *  COIN_RESTRICT pi=updateForDjs.denseVector();
2901  const double *COIN_RESTRICT piWeights = updateForWeights.denseVector();
2902  const unsigned char *COIN_RESTRICT internalStatus = model_->internalStatus();
2903  double *COIN_RESTRICT reducedCost = model_->djRegion();
2904  double tolerance = model_->currentDualTolerance();
2905  // use spareColumn to track new infeasibilities
2906  int *COIN_RESTRICT newInfeas = spareColumn1.getIndices() + first;
2907  int numberNewInfeas = 0;
2908  // we can't really trust infeasibilities if there is dual error
2909  // this coding has to mimic coding in checkDualSolution
2910  double error = CoinMin(1.0e-2, model_->largestDualError());
2911  // allow tolerance at least slightly bigger than standard
2912  tolerance = tolerance + error;
2913  double mult[2] = { 1.0, -1.0 };
2914  bool updateWeights = devex != 0.0;
2915  int numberSlacks = updateForTableauRow.getNumElements();
2916  const int *COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
2917#define isReference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0)
2918  // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor*
2919  assert(scaleFactor);
2920  if (!iBlock) {
2921    if (scaleFactor != COIN_DBL_MAX) {
2922      for (int i = 0; i < numberSlacks; i++) {
2923        int iRow = piIndexTableau[i];
2924        int iStatus = internalStatus[iRow] & 7;
2925        if (iStatus < 4) {
2926          double value = reducedCost[iRow];
2927          value += scaleFactor * piTableau[iRow]; // sign?
2928          reducedCost[iRow] = value;
2929          value *= mult[iStatus];
2930          if (value < 0.0) {
2931            if (!infeasibilities[iRow])
2932              newInfeas[numberNewInfeas++] = iRow;
2933#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
2934            infeasibilities[iRow] = value * value * CLP_PRIMAL_SLACK_MULTIPLIER;
2935#else
2936            infeasibilities[iRow] = value * value;
2937#endif
2938          } else {
2939            if (infeasibilities[iRow])
2940              infeasibilities[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
2941          }
2942        }
2943      }
2944    }
2945    if (updateWeights) {
2946      for (int i = 0; i < numberSlacks; i++) {
2947        int iRow = piIndexTableau[i];
2948        double modification = piWeights[iRow];
2949        double thisWeight = weights[iRow];
2950        double pivot = piTableau[iRow];
2951        double pivotSquared = pivot * pivot;
2952        thisWeight += pivotSquared * devex - pivot * modification;
2953        if (thisWeight < DEVEX_TRY_NORM) {
2954          if (referenceIn < 0.0) {
2955            // steepest
2956            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2957          } else {
2958            // exact
2959            thisWeight = referenceIn * pivotSquared;
2960            if (isReference(iRow))
2961              thisWeight += 1.0;
2962            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2963          }
2964        }
2965        weights[iRow] = thisWeight;
2966      }
2967    }
2968    first = maximumRows;
2969  }
2970  double zeroTolerance = model_->zeroTolerance();
2971  double freeTolerance = FREE_ACCEPT * tolerance;
2972  ;
2973  // get matrix data pointers
2974  const int *COIN_RESTRICT row = matrix_->getIndices();
2975  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
2976  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
2977  // get tableau row
2978  int *COIN_RESTRICT index = updateForTableauRow.getIndices() + first;
2979  double *COIN_RESTRICT array = updateForTableauRow.denseVector();
2980  int numberRows = model_->numberRows();
2981  const CoinBigIndex *COIN_RESTRICT rowStart = rowStart_ + iBlock * numberRows;
2982  const CoinBigIndex *COIN_RESTRICT rowEnd = rowStart + numberRows;
2983  const double *COIN_RESTRICT element = element_;
2984  const int *COIN_RESTRICT column = column_;
2985  int numberInRow = 0;
2986  if (numberSlacks > 2) {
2987    for (int i = 0; i < numberSlacks; i++) {
2988      int iRow = piIndexTableau[i];
2989      double piValue = piTableau[iRow];
2990      CoinBigIndex end = rowEnd[iRow];
2991      for (CoinBigIndex j = rowStart[iRow]; j < end; j++) {
2992        int iSequence = column[j];
2993        double value = element[j] * piValue;
2994        double oldValue = array[iSequence];
2995        value += oldValue;
2996        if (!oldValue) {
2997          array[iSequence] = value;
2998          index[numberInRow++] = iSequence;
2999        } else if (value) {
3000          array[iSequence] = value;
3001        } else {
3002          array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
3003        }
3004      }
3005    }
3006  } else if (numberSlacks == 2) {
3007    int iRow0 = piIndexTableau[0];
3008    int iRow1 = piIndexTableau[1];
3009    CoinBigIndex end0 = rowEnd[iRow0];
3010    CoinBigIndex end1 = rowEnd[iRow1];
3011    if (end0 - rowStart[iRow0] > end1 - rowStart[iRow1]) {
3012      int temp = iRow0;
3013      iRow0 = iRow1;
3014      iRow1 = temp;
3015    }
3016    CoinBigIndex start = rowStart[iRow0];
3017    CoinBigIndex end = rowEnd[iRow0];
3018    double piValue = piTableau[iRow0];
3019    for (CoinBigIndex j = start; j < end; j++) {
3020      int iSequence = column[j];
3021      double value = element[j];
3022      array[iSequence] = piValue * value;
3023      index[numberInRow++] = iSequence;
3024    }
3025    start = rowStart[iRow1];
3026    end = rowEnd[iRow1];
3027    piValue = piTableau[iRow1];
3028    for (CoinBigIndex j = start; j < end; j++) {
3029      int iSequence = column[j];
3030      double value = element[j];
3031      value *= piValue;
3032      if (!array[iSequence]) {
3033        array[iSequence] = value;
3034        index[numberInRow++] = iSequence;
3035      } else {
3036        value += array[iSequence];
3037        if (value)
3038          array[iSequence] = value;
3039        else
3040          array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
3041      }
3042    }
3043  } else {
3044    int iRow0 = piIndexTableau[0];
3045    CoinBigIndex start = rowStart[iRow0];
3046    CoinBigIndex end = rowEnd[iRow0];
3047    double piValue = piTableau[iRow0];
3048    for (CoinBigIndex j = start; j < end; j++) {
3049      int iSequence = column[j];
3050      double value = element[j];
3051      array[iSequence] = piValue * value;
3052      index[numberInRow++] = iSequence;
3053    }
3054  }
3055  bool updateDjs = scaleFactor != COIN_DBL_MAX;
3056  for (int iLook = 0; iLook < numberInRow; iLook++) {
3057    int iSequence = index[iLook];
3058    unsigned char iStatus = internalStatus[iSequence] & 7;
3059    if (iStatus < 6) {
3060      double tableauValue = array[iSequence];
3061      array[iSequence] = 0.0;
3062      double djValue = reducedCost[iSequence];
3063      djValue += tableauValue * scaleFactor; // sign?
3064      if (updateDjs) {
3065        reducedCost[iSequence] = djValue;
3066        if (iStatus < 4) {
3067          djValue *= mult[iStatus];
3068          if (djValue < 0.0) {
3069            if (!infeasibilities[iSequence])
3070              newInfeas[numberNewInfeas++] = iSequence;
3071            infeasibilities[iSequence] = djValue * djValue;
3072          } else {
3073            if (infeasibilities[iSequence])
3074              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
3075          }
3076        } else {
3077          if (fabs(djValue) > freeTolerance) {
3078            // we are going to bias towards free (but only if reasonable)
3079            djValue *= FREE_BIAS;
3080            if (!infeasibilities[iSequence])
3081              newInfeas[numberNewInfeas++] = iSequence;
3082            // store square in list
3083            infeasibilities[iSequence] = djValue * djValue;
3084          } else {
3085            if (infeasibilities[iSequence])
3086              infeasibilities[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
3087          }
3088        }
3089      }
3090      if (fabs(tableauValue) > zeroTolerance && updateWeights) {
3091        CoinBigIndex start = columnStart[iSequence];
3092        CoinBigIndex next = columnStart[iSequence + 1];
3093        double modification = 0.0;
3094        for (CoinBigIndex j = start; j < next; j++) {
3095          int jRow = row[j];
3096          modification += piWeights[jRow] * elementByColumn[j];
3097        }
3098        double thisWeight = weights[iSequence];
3099        double pivot = tableauValue;
3100        double pivotSquared = pivot * pivot;
3101        thisWeight += pivotSquared * devex - pivot * modification;
3102        if (thisWeight < DEVEX_TRY_NORM) {
3103          if (referenceIn < 0.0) {
3104            // steepest
3105            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
3106          } else {
3107            // exact
3108            thisWeight = referenceIn * pivotSquared;
3109            if (isReference(iSequence))
3110              thisWeight += 1.0;
3111            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
3112          }
3113        }
3114        weights[iSequence] = thisWeight;
3115      }
3116    } else {
3117      array[iSequence] = 0.0;
3118    }
3119  }
3120  spareColumn1.setTempNumElementsPartition(iBlock, numberNewInfeas);
3121  int bestSequence = -1;
3122#if 0
3123  if (!iBlock)
3124    first=0;
3125  // not at present - maybe later
3126  double bestDj=tolerance*tolerance;
3127  for (int iSequence=first;iSequence<last;iSequence++) {
3128    if (infeasibilities[iSequence]>bestDj*weights[iSequence]) {
3129      int iStatus=internalStatus[iSequence]&7;
3130      assert (iStatus<6);
3131      bestSequence=iSequence;
3132      bestDj=infeasibilities[iSequence]/weights[iSequence];
3133    }
3134  }
3135#endif
3136  return bestSequence;
3137}
3138#if 0
3139/* Chooses best weighted dj
3140 */
3141int
3142AbcMatrix::chooseBestDj(int iBlock,const CoinIndexedVector & infeasibilities,
3143                        const double * weights) const
3144{
3145  return 0;
3146}
3147#endif
3148/* does steepest edge double or triple update
3149   If scaleFactor!=0 then use with tableau row to update djs
3150   otherwise use updateForDjs
3151   Returns best
3152*/
3153int AbcMatrix::primalColumnDouble(CoinPartitionedVector &updateForTableauRow,
3154  CoinPartitionedVector &updateForDjs,
3155  const CoinIndexedVector &updateForWeights,
3156  CoinPartitionedVector &spareColumn1,
3157  CoinIndexedVector &infeasible,
3158  double referenceIn, double devex,
3159  // Array for exact devex to say what is in reference framework
3160  unsigned int *reference,
3161  double *weights, double scaleFactor) const
3162{
3163  //int maximumRows = model_->maximumAbcNumberRows();
3164  // rebalance
3165  rebalance();
3166#ifdef PRICE_IN_ABC_MATRIX
3167  int which[NUMBER_BLOCKS];
3168#endif
3169  double *infeasibilities = infeasible.denseVector();
3170  int bestSequence = -1;
3171  // see if worth using row copy
3172  int maximumRows = model_->maximumAbcNumberRows();
3173  int number = updateForTableauRow.getNumElements();
3174#ifdef GCC_4_9
3175  assert(number);
3176#else
3177  if (!number) {
3178    printf("Null tableau row!\n");
3179  }
3180#endif
3181  bool useRowCopy = (gotRowCopy() && (number << 2) < maximumRows);
3182  //useRowCopy=false;
3183  if (!scaleFactor)
3184    useRowCopy = false; // look again later
3185  const int *starts;
3186  int numberBlocks;
3187  if (useRowCopy) {
3188    starts = blockStart_;
3189    numberBlocks = numberRowBlocks_;
3190  } else {
3191    starts = startColumnBlock_;
3192    numberBlocks = numberColumnBlocks_;
3193  }
3194  if (useRowCopy) {
3195    for (int i = 0; i < numberBlocks; i++)
3196#ifdef PRICE_IN_ABC_MATRIX
3197      which[i] =
3198#endif
3199        cilk_spawn primalColumnSparseDouble(i, updateForTableauRow, updateForDjs, updateForWeights,
3200          spareColumn1,
3201          infeasibilities, referenceIn, devex, reference, weights, scaleFactor);
3202    cilk_sync;
3203  } else {
3204    for (int i = 0; i < numberBlocks; i++)
3205#ifdef PRICE_IN_ABC_MATRIX
3206      which[i] =
3207#endif
3208        cilk_spawn primalColumnDouble(i, updateForTableauRow, updateForDjs, updateForWeights,
3209          spareColumn1,
3210          infeasibilities, referenceIn, devex, reference, weights, scaleFactor);
3211    cilk_sync;
3212  }
3213#ifdef PRICE_IN_ABC_MATRIX
3214  double bestValue = model_->dualTolerance();
3215  int sequenceIn[8] = { -1, -1, -1, -1, -1, -1, -1, -1 };
3216#endif
3217  assert(numberColumnBlocks_ <= 8);
3218#ifdef PRICE_IN_ABC_MATRIX
3219  double weightedDj[8];
3220  int nPossible = 0;
3221#endif
3222  //const double * reducedCost = model_->djRegion();
3223  // use spareColumn to track new infeasibilities
3224  int numberInfeas = infeasible.getNumElements();
3225  int *COIN_RESTRICT infeas = infeasible.getIndices();
3226  const int *COIN_RESTRICT newInfeasAll = spareColumn1.getIndices();
3227  for (int i = 0; i < numberColumnBlocks_; i++) {
3228    const int *COIN_RESTRICT newInfeas = newInfeasAll + starts[i];
3229    int numberNewInfeas = spareColumn1.getNumElements(i);
3230    spareColumn1.setTempNumElementsPartition(i, 0);
3231    for (int j = 0; j < numberNewInfeas; j++)
3232      infeas[numberInfeas++] = newInfeas[j];
3233#ifdef PRICE_IN_ABC_MATRIX
3234    if (which[i] >= 0) {
3235      int iSequence = which[i];
3236      double value = fabs(infeasibilities[iSequence] / weights[iSequence]);
3237      if (value > bestValue) {
3238        bestValue = value;
3239        bestSequence = which[i];
3240      }
3241      weightedDj[nPossible] = -value;
3242      sequenceIn[nPossible++] = iSequence;
3243    }
3244#endif
3245  }
3246  infeasible.setNumElements(numberInfeas);
3247#ifdef PRICE_IN_ABC_MATRIX
3248  CoinSort_2(weightedDj, weightedDj + nPossible, sequenceIn);
3249  //model_->setMultipleSequenceIn(sequenceIn);
3250#endif
3251  return bestSequence;
3252}
3253// Partial pricing
3254void AbcMatrix::partialPricing(double startFraction, double endFraction,
3255  int &bestSequence, int &numberWanted)
3256{
3257  numberWanted = currentWanted_;
3258  int maximumRows = model_->maximumAbcNumberRows();
3259  // get matrix data pointers
3260  const int *COIN_RESTRICT row = matrix_->getIndices();
3261  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
3262  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
3263  int numberColumns = model_->numberColumns();
3264  int start = static_cast< int >(startFraction * numberColumns);
3265  int end = CoinMin(static_cast< int >(endFraction * numberColumns + 1), numberColumns);
3266  // adjust
3267  start += maximumRows;
3268  end += maximumRows;
3269  double tolerance = model_->currentDualTolerance();
3270  double *reducedCost = model_->djRegion();
3271  const double *duals = model_->dualRowSolution();
3272  const double *cost = model_->costRegion();
3273  double bestDj;
3274  if (bestSequence >= 0)
3275    bestDj = fabs(reducedCost[bestSequence]);
3276  else
3277    bestDj = tolerance;
3278  int sequenceOut = model_->sequenceOut();
3279  int saveSequence = bestSequence;
3280  int lastScan = minimumObjectsScan_ < 0 ? end : start + minimumObjectsScan_;
3281  int minNeg = minimumGoodReducedCosts_ == -1 ? numberWanted : minimumGoodReducedCosts_;
3282  for (int iSequence = start; iSequence < end; iSequence++) {
3283    if (iSequence != sequenceOut) {
3284      double value;
3285      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
3286      switch (status) {
3287      case AbcSimplex::basic:
3288      case AbcSimplex::isFixed:
3289        break;
3290      case AbcSimplex::isFree:
3291      case AbcSimplex::superBasic:
3292        value = cost[iSequence];
3293        for (CoinBigIndex j = columnStart[iSequence];
3294             j < columnStart[iSequence + 1]; j++) {
3295          int jRow = row[j];
3296          value -= duals[jRow] * elementByColumn[j];
3297        }
3298        value = fabs(value);
3299        if (value > FREE_ACCEPT * tolerance) {
3300          numberWanted--;
3301          // we are going to bias towards free (but only if reasonable)
3302          value *= FREE_BIAS;
3303          if (value > bestDj) {
3304            // check flagged variable and correct dj
3305            if (!model_->flagged(iSequence)) {
3306              bestDj = value;
3307              bestSequence = iSequence;
3308            } else {
3309              // just to make sure we don't exit before got something
3310              numberWanted++;
3311            }
3312          }
3313        }
3314        break;
3315      case AbcSimplex::atUpperBound:
3316        value = cost[iSequence];
3317        // scaled
3318        for (CoinBigIndex j = columnStart[iSequence];
3319             j < columnStart[iSequence + 1]; j++) {
3320          int jRow = row[j];
3321          value -= duals[jRow] * elementByColumn[j];
3322        }
3323        if (value > tolerance) {
3324          numberWanted--;
3325          if (value > bestDj) {
3326            // check flagged variable and correct dj
3327            if (!model_->flagged(iSequence)) {
3328              bestDj = value;
3329              bestSequence = iSequence;
3330            } else {
3331              // just to make sure we don't exit before got something
3332              numberWanted++;
3333            }
3334          }
3335        }
3336        break;
3337      case AbcSimplex::atLowerBound:
3338        value = cost[iSequence];
3339        for (CoinBigIndex j = columnStart[iSequence];
3340             j < columnStart[iSequence + 1]; j++) {
3341          int jRow = row[j];
3342          value -= duals[jRow] * elementByColumn[j];
3343        }
3344        value = -value;
3345        if (value > tolerance) {
3346          numberWanted--;
3347          if (value > bestDj) {
3348            // check flagged variable and correct dj
3349            if (!model_->flagged(iSequence)) {
3350              bestDj = value;
3351              bestSequence = iSequence;
3352            } else {
3353              // just to make sure we don't exit before got something
3354              numberWanted++;
3355            }
3356          }
3357        }
3358        break;
3359      }
3360    }
3361    if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
3362      // give up
3363      break;
3364    }
3365    if (!numberWanted)
3366      break;
3367  }
3368  if (bestSequence != saveSequence) {
3369    // recompute dj
3370    double value = cost[bestSequence];
3371    for (CoinBigIndex j = columnStart[bestSequence];
3372         j < columnStart[bestSequence + 1]; j++) {
3373      int jRow = row[j];
3374      value -= duals[jRow] * elementByColumn[j];
3375    }
3376    reducedCost[bestSequence] = value;
3377    savedBestSequence_ = bestSequence;
3378    savedBestDj_ = reducedCost[savedBestSequence_];
3379  }
3380  currentWanted_ = numberWanted;
3381}
3382/* Return <code>x *A</code> in <code>z</code> but
3383   just for indices Already in z.
3384   Note - z always packed mode */
3385void AbcMatrix::subsetTransposeTimes(const CoinIndexedVector &x,
3386  CoinIndexedVector &z) const
3387{
3388  int maximumRows = model_->maximumAbcNumberRows();
3389  // get matrix data pointers
3390  const int *COIN_RESTRICT row = matrix_->getIndices();
3391  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
3392  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
3393  const double *COIN_RESTRICT other = x.denseVector();
3394  int numberNonZero = z.getNumElements();
3395  int *COIN_RESTRICT index = z.getIndices();
3396  double *COIN_RESTRICT array = z.denseVector();
3397  int numberRows = model_->numberRows();
3398  for (int i = 0; i < numberNonZero; i++) {
3399    int iSequence = index[i];
3400    if (iSequence >= numberRows) {
3401      CoinBigIndex start = columnStart[iSequence];
3402      CoinBigIndex next = columnStart[iSequence + 1];
3403      double value = 0.0;
3404      for (CoinBigIndex j = start; j < next; j++) {
3405        int jRow = row[j];
3406        value -= other[jRow] * elementByColumn[j];
3407      }
3408      array[i] = value;
3409    } else {
3410      array[i] = -other[iSequence];
3411    }
3412  }
3413}
3414// Return <code>-x *A</code> in <code>z</code>
3415void AbcMatrix::transposeTimes(const CoinIndexedVector &x,
3416  CoinIndexedVector &z) const
3417{
3418  int maximumRows = model_->maximumAbcNumberRows();
3419  // get matrix data pointers
3420  const int *COIN_RESTRICT row = matrix_->getIndices();
3421  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
3422  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
3423  const double *COIN_RESTRICT other = x.denseVector();
3424  int numberNonZero = 0;
3425  int *COIN_RESTRICT index = z.getIndices();
3426  double *COIN_RESTRICT array = z.denseVector();
3427  int numberColumns = model_->numberColumns();
3428  double zeroTolerance = model_->zeroTolerance();
3429  for (int iSequence = maximumRows; iSequence < maximumRows + numberColumns; iSequence++) {
3430    CoinBigIndex start = columnStart[iSequence];
3431    CoinBigIndex next = columnStart[iSequence + 1];
3432    double value = 0.0;
3433    for (CoinBigIndex j = start; j < next; j++) {
3434      int jRow = row[j];
3435      value -= other[jRow] * elementByColumn[j];
3436    }
3437    if (fabs(value) > zeroTolerance) {
3438      // TEMP
3439      array[iSequence - maximumRows] = value;
3440      index[numberNonZero++] = iSequence - maximumRows;
3441    }
3442  }
3443  z.setNumElements(numberNonZero);
3444}
3445/* Return A * scalar(+-1) *x + y</code> in <code>y</code>.
3446   @pre <code>x</code> must be of size <code>numRows()</code>
3447   @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
3448void AbcMatrix::transposeTimesNonBasic(double scalar,
3449  const double *pi, double *y) const
3450{
3451  int numberTotal = model_->numberTotal();
3452  int maximumRows = model_->maximumAbcNumberRows();
3453  assert(scalar == -1.0);
3454  // get matrix data pointers
3455  const int *COIN_RESTRICT row = matrix_->getIndices();
3456  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
3457  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
3458  int numberRows = model_->numberRows();
3459  const unsigned char *COIN_RESTRICT status = model_->internalStatus();
3460  for (int i = 0; i < numberRows; i++) {
3461    y[i] += scalar * pi[i];
3462  }
3463  for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
3464    unsigned char type = status[iColumn] & 7;
3465    if (type < 6) {
3466      CoinBigIndex start = columnStart[iColumn];
3467      CoinBigIndex next = columnStart[iColumn + 1];
3468      double value = 0.0;
3469      for (CoinBigIndex j = start; j < next; j++) {
3470        int jRow = row[j];
3471        value += scalar * pi[jRow] * elementByColumn[j];
3472      }
3473      y[iColumn] += value;
3474    } else {
3475      y[iColumn] = 0.0; // may not be but .....
3476    }
3477  }
3478}
3479/* Return y - A * x</code> in <code>y</code>.
3480   @pre <code>x</code> must be of size <code>numRows()</code>
3481   @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
3482void AbcMatrix::transposeTimesAll(const double *pi, double *y) const
3483{
3484  int numberTotal = model_->numberTotal();
3485  int maximumRows = model_->maximumAbcNumberRows();
3486  // get matrix data pointers
3487  const int *COIN_RESTRICT row = matrix_->getIndices();
3488  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - maximumRows;
3489  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
3490  int numberRows = model_->numberRows();
3491  for (int i = 0; i < numberRows; i++) {
3492    y[i] -= pi[i];
3493  }
3494  for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
3495    CoinBigIndex start = columnStart[iColumn];
3496    CoinBigIndex next = columnStart[iColumn + 1];
3497    double value = 0.0;
3498    for (CoinBigIndex j = start; j < next; j++) {
3499      int jRow = row[j];
3500      value += pi[jRow] * elementByColumn[j];
3501    }
3502    y[iColumn] -= value;
3503  }
3504}
3505/* Return y  + A * scalar(+-1) *x</code> in <code>y</code>.
3506   @pre <code>x</code> must be of size <code>numRows()</code>
3507   @pre <code>y</code> must be of size <code>numRows()</code> */
3508void AbcMatrix::transposeTimesBasic(double scalar,
3509  const double *pi, double *y) const
3510{
3511  assert(scalar == -1.0);
3512  int numberRows = model_->numberRows();
3513  //AbcMemset0(y,numberRows);
3514  // get matrix data pointers
3515  const int *COIN_RESTRICT row = matrix_->getIndices();
3516  const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts() - model_->maximumAbcNumberRows();
3517  const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
3518  const int *COIN_RESTRICT pivotVariable = model_->pivotVariable();
3519  for (int i = 0; i < numberRows; i++) {
3520    int iSequence = pivotVariable[i];
3521    double value;
3522    if (iSequence >= numberRows) {
3523      CoinBigIndex start = columnStart[iSequence];
3524      CoinBigIndex next = columnStart[iSequence + 1];
3525      value = 0.0;
3526      for (CoinBigIndex j = start; j < next; j++) {
3527        int jRow = row[j];
3528        value += pi[jRow] * elementByColumn[j];
3529      }
3530    } else {
3531      value = pi[iSequence];
3532    }
3533    y[i] += scalar * value;
3534  }
3535}
3536// Adds multiple of a column (or slack) into an CoinIndexedvector
3537void AbcMatrix::add(CoinIndexedVector &rowArray, int iSequence, double multiplier) const
3538{
3539  int maximumRows = model_->maximumAbcNumberRows();
3540  if (iSequence >= maximumRows) {
3541    const int *COIN_RESTRICT row = matrix_->getIndices();
3542    iSequence -= maximumRows;
3543    const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts();
3544    const double *COIN_RESTRICT elementByColumn = matrix_->getElements();
3545    CoinBigIndex start = columnStart[iSequence];
3546    CoinBigIndex next = columnStart[iSequence + 1];
3547    for (CoinBigIndex j = start; j < next; j++) {
3548      int jRow = row[j];
3549      rowArray.quickAdd(jRow, elementByColumn[j] * multiplier);
3550    }
3551  } else {
3552    rowArray.quickAdd(iSequence, multiplier);
3553  }
3554}
3555/* Unpacks a column into an CoinIndexedVector
3556 */
3557void AbcMatrix::unpack(CoinIndexedVector &usefulArray,
3558  int sequence) const
3559{
3560  usefulArray.clear();
3561  int maximumRows = model_->maximumAbcNumberRows();
3562  if (sequence < maximumRows) {
3563    //slack
3564    usefulArray.insert(sequence, 1.0);
3565  } else {
3566    // column
3567    const CoinBigIndex *COIN_RESTRICT columnStart = matrix()->getVectorStarts() - maximumRows;
3568    CoinBigIndex start = columnStart[sequence];
3569    CoinBigIndex end = columnStart[sequence + 1];
3570    const int *COIN_RESTRICT row = matrix()->getIndices() + start;
3571    const double *COIN_RESTRICT elementByColumn = matrix()->getElements() + start;
3572    int *COIN_RESTRICT index = usefulArray.getIndices();
3573    double *COIN_RESTRICT array = usefulArray.denseVector();
3574    int numberNonZero = end - start;
3575    for (int j = 0; j < numberNonZero; j++) {
3576      int iRow = row[j];
3577      index[j] = iRow;
3578      array[iRow] = elementByColumn[j];
3579    }
3580    usefulArray.setNumElements(numberNonZero);
3581  }
3582}
3583// Row starts
3584CoinBigIndex *
3585AbcMatrix::rowStart() const
3586{
3587  return rowStart_;
3588}
3589// Row ends
3590CoinBigIndex *
3591AbcMatrix::rowEnd() const
3592{
3593  //if (numberRowBlocks_<2) {
3594  //return rowStart_+1;
3595  //} else {
3596  return rowStart_ + model_->numberRows() * (numberRowBlocks_ + 1);
3597  //}
3598}
3599// Row elements
3600double *
3601AbcMatrix::rowElements() const
3602{
3603  return element_;
3604}
3605// Row columnss
3606CoinSimplexInt *
3607AbcMatrix::rowColumns() const
3608{
3609  return column_;
3610}
3611
3612/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
3613*/
Note: See TracBrowser for help on using the repository browser.