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

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