source: trunk/Clp/src/ClpCholeskyDense.cpp @ 1371

Last change on this file since 1371 was 1371, checked in by forrest, 11 years ago

changes to try and make faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.5 KB
Line 
1/* $Id: ClpCholeskyDense.cpp 1371 2009-06-12 16:29:04Z forrest $ */
2/* Copyright (C) 2003, International Business Machines Corporation
3   and others.  All Rights Reserved. */
4#include "CoinPragma.hpp"
5#include "CoinHelperFunctions.hpp"
6#include "ClpHelperFunctions.hpp"
7
8#include "ClpInterior.hpp"
9#include "ClpCholeskyDense.hpp"
10#include "ClpMessage.hpp"
11#include "ClpQuadraticObjective.hpp"
12
13/*#############################################################################*/
14/* Constructors / Destructor / Assignment*/
15/*#############################################################################*/
16
17/*-------------------------------------------------------------------*/
18/* Default Constructor */
19/*-------------------------------------------------------------------*/
20ClpCholeskyDense::ClpCholeskyDense () 
21  : ClpCholeskyBase(),
22    borrowSpace_(false)
23{
24  type_=11;;
25}
26
27/*-------------------------------------------------------------------*/
28/* Copy constructor */
29/*-------------------------------------------------------------------*/
30ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs) 
31  : ClpCholeskyBase(rhs),
32    borrowSpace_(rhs.borrowSpace_)
33{
34  assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/
35}
36
37
38/*-------------------------------------------------------------------*/
39/* Destructor */
40/*-------------------------------------------------------------------*/
41ClpCholeskyDense::~ClpCholeskyDense ()
42{
43  if (borrowSpace_) {
44    /* set NULL*/
45    sparseFactor_=NULL;
46    workDouble_=NULL;
47    diagonal_=NULL;
48  }
49}
50
51/*----------------------------------------------------------------*/
52/* Assignment operator */
53/*-------------------------------------------------------------------*/
54ClpCholeskyDense &
55ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
56{
57  if (this != &rhs) {
58    assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/
59    ClpCholeskyBase::operator=(rhs);
60    borrowSpace_=rhs.borrowSpace_;
61  }
62  return *this;
63}
64/*-------------------------------------------------------------------*/
65/* Clone*/
66/*-------------------------------------------------------------------*/
67ClpCholeskyBase * ClpCholeskyDense::clone() const
68{
69  return new ClpCholeskyDense(*this);
70}
71/* If not power of 2 then need to redo a bit*/
72#define BLOCK 16
73#define BLOCKSHIFT 4
74/* Block unroll if power of 2 and at least 8*/
75#define BLOCKUNROLL
76
77#define BLOCKSQ ( BLOCK*BLOCK )
78#define BLOCKSQSHIFT ( BLOCKSHIFT+BLOCKSHIFT )
79#define number_blocks(x) (((x)+BLOCK-1)>>BLOCKSHIFT)
80#define number_rows(x) ((x)<<BLOCKSHIFT)
81#define number_entries(x) ((x)<<BLOCKSQSHIFT)
82/* Gets space */
83int 
84ClpCholeskyDense::reserveSpace(const ClpCholeskyBase * factor, int numberRows) 
85{
86  numberRows_ = numberRows;
87  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
88  /* allow one stripe extra*/
89  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
90  sizeFactor_=numberBlocks*BLOCKSQ;
91  /*#define CHOL_COMPARE*/
92#ifdef CHOL_COMPARE 
93  sizeFactor_ += 95000;
94#endif
95  if (!factor) {
96    sparseFactor_ = new longDouble [sizeFactor_];
97    rowsDropped_ = new char [numberRows_];
98    memset(rowsDropped_,0,numberRows_);
99    workDouble_ = new longDouble[numberRows_];
100    diagonal_ = new longDouble[numberRows_];
101  } else {
102    borrowSpace_=true;
103    int numberFull = factor->numberRows();
104    sparseFactor_ = factor->sparseFactor()+(factor->size()-sizeFactor_);
105    workDouble_ = factor->workDouble() + (numberFull-numberRows_);
106    diagonal_ = factor->diagonal() + (numberFull-numberRows_);
107  }
108  numberRowsDropped_=0;
109  return 0;
110}
111/* Returns space needed */
112CoinBigIndex
113ClpCholeskyDense::space( int numberRows) const
114{
115 int numberBlocks = (numberRows+BLOCK-1)>>BLOCKSHIFT;
116 /* allow one stripe extra*/
117  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
118  CoinBigIndex sizeFactor=numberBlocks*BLOCKSQ;
119#ifdef CHOL_COMPARE 
120  sizeFactor += 95000;
121#endif
122  return sizeFactor;
123 }
124/* Orders rows and saves pointer to matrix.and model */
125int 
126ClpCholeskyDense::order(ClpInterior * model) 
127{
128  model_=model;
129  int numberRows;
130  int numberRowsModel = model_->numberRows();
131  int numberColumns = model_->numberColumns();
132  if (!doKKT_) {
133    numberRows = numberRowsModel;
134  } else {
135    numberRows = 2*numberRowsModel+numberColumns;
136  }
137  reserveSpace(NULL,numberRows);
138  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
139  return 0;
140}
141/* Does Symbolic factorization given permutation.
142   This is called immediately after order.  If user provides this then
143   user must provide factorize and solve.  Otherwise the default factorization is used
144   returns non-zero if not enough memory */
145int 
146ClpCholeskyDense::symbolic()
147{
148  return 0;
149}
150/* Factorize - filling in rowsDropped and returning number dropped */
151int 
152ClpCholeskyDense::factorize(const CoinWorkDouble * diagonal, int * rowsDropped) 
153{
154  assert (!borrowSpace_);
155  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
156  const int * columnLength = model_->clpMatrix()->getVectorLengths();
157  const int * row = model_->clpMatrix()->getIndices();
158  const double * element = model_->clpMatrix()->getElements();
159  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
160  const int * rowLength = rowCopy_->getVectorLengths();
161  const int * column = rowCopy_->getIndices();
162  const double * elementByRow = rowCopy_->getElements();
163  int numberColumns=model_->clpMatrix()->getNumCols();
164  CoinZeroN(sparseFactor_,sizeFactor_);
165  /*perturbation*/
166  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
167  perturbation=perturbation*perturbation;
168  if (perturbation>1.0) {
169#ifdef COIN_DEVELOP
170    /*if (model_->model()->logLevel()&4) */
171      std::cout <<"large perturbation "<<perturbation<<std::endl;
172#endif
173    perturbation=CoinSqrt(perturbation);;
174    perturbation=1.0;
175  }
176  int iRow;
177  int newDropped=0;
178  CoinWorkDouble largest=1.0;
179  CoinWorkDouble smallest=COIN_DBL_MAX;
180  CoinWorkDouble delta2 = model_->delta(); /* add delta*delta to diagonal*/
181  delta2 *= delta2;
182  if (!doKKT_) {
183    longDouble * work = sparseFactor_;
184    work--; /* skip diagonal*/
185    int addOffset=numberRows_-1;
186    const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
187    /* largest in initial matrix*/
188    CoinWorkDouble largest2=1.0e-20;
189    for (iRow=0;iRow<numberRows_;iRow++) {
190      if (!rowsDropped_[iRow]) {
191        CoinBigIndex startRow=rowStart[iRow];
192        CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
193        CoinWorkDouble diagonalValue = diagonalSlack[iRow]+delta2;
194        for (CoinBigIndex k=startRow;k<endRow;k++) {
195          int iColumn=column[k];
196          CoinBigIndex start=columnStart[iColumn];
197          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
198          CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
199          for (CoinBigIndex j=start;j<end;j++) {
200            int jRow=row[j];
201            if (!rowsDropped_[jRow]) {
202              if (jRow>iRow) {
203                CoinWorkDouble value=element[j]*multiplier;
204                work[jRow] += value;
205              } else if (jRow==iRow) {
206                CoinWorkDouble value=element[j]*multiplier;
207                diagonalValue += value;
208              }
209            }
210          }
211        } 
212        for (int j=iRow+1;j<numberRows_;j++) 
213          largest2 = CoinMax(largest2,CoinAbs(work[j]));
214        diagonal_[iRow]=diagonalValue;
215        largest2 = CoinMax(largest2,CoinAbs(diagonalValue));
216      } else {
217        /* dropped*/
218        diagonal_[iRow]=1.0;
219      }
220      addOffset--;
221      work += addOffset;
222    }
223    /*check sizes*/
224    largest2*=1.0e-20;
225    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
226    int numberDroppedBefore=0;
227    for (iRow=0;iRow<numberRows_;iRow++) {
228      int dropped=rowsDropped_[iRow];
229      /* Move to int array*/
230      rowsDropped[iRow]=dropped;
231      if (!dropped) {
232        CoinWorkDouble diagonal = diagonal_[iRow];
233        if (diagonal>largest2) {
234          diagonal_[iRow]=diagonal+perturbation;
235        } else {
236          diagonal_[iRow]=diagonal+perturbation;
237          rowsDropped[iRow]=2;
238          numberDroppedBefore++;
239        } 
240      }
241    }
242    doubleParameters_[10]=CoinMax(1.0e-20,largest);
243    integerParameters_[20]=0;
244    doubleParameters_[3]=0.0;
245    doubleParameters_[4]=COIN_DBL_MAX;
246    integerParameters_[34]=0; /* say all must be positive*/
247#ifdef CHOL_COMPARE 
248    if (numberRows_<200)
249      factorizePart3(rowsDropped);
250#endif
251    factorizePart2(rowsDropped);
252    newDropped=integerParameters_[20]+numberDroppedBefore;
253    largest=doubleParameters_[3];
254    smallest=doubleParameters_[4];
255    if (model_->messageHandler()->logLevel()>1) 
256      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
257    choleskyCondition_=largest/smallest;
258    /*drop fresh makes some formADAT easier*/
259    if (newDropped||numberRowsDropped_) {
260      newDropped=0;
261      for (int i=0;i<numberRows_;i++) {
262        char dropped = static_cast<char>(rowsDropped[i]);
263        rowsDropped_[i]=dropped;
264        if (dropped==2) {
265          /*dropped this time*/
266          rowsDropped[newDropped++]=i;
267          rowsDropped_[i]=0;
268        } 
269      } 
270      numberRowsDropped_=newDropped;
271      newDropped=-(2+newDropped);
272    } 
273  } else {
274    /* KKT*/
275    CoinPackedMatrix * quadratic = NULL;
276    ClpQuadraticObjective * quadraticObj = 
277      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
278    if (quadraticObj) 
279      quadratic = quadraticObj->quadraticObjective();
280    /* matrix*/
281    int numberRowsModel = model_->numberRows();
282    int numberColumns = model_->numberColumns();
283    int numberTotal = numberColumns + numberRowsModel;
284    longDouble * work = sparseFactor_;
285    work--; /* skip diagonal*/
286    int addOffset=numberRows_-1;
287    int iColumn;
288    if (!quadratic) {
289      for (iColumn=0;iColumn<numberColumns;iColumn++) {
290        CoinWorkDouble value = diagonal[iColumn];
291        if (CoinAbs(value)>1.0e-100) {
292          value = 1.0/value;
293          largest = CoinMax(largest,CoinAbs(value));
294          diagonal_[iColumn] = -value;
295          CoinBigIndex start=columnStart[iColumn];
296          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
297          for (CoinBigIndex j=start;j<end;j++) {
298            /*choleskyRow_[numberElements]=row[j]+numberTotal;*/
299            /*sparseFactor_[numberElements++]=element[j];*/
300            work[row[j]+numberTotal]=element[j];
301            largest = CoinMax(largest,CoinAbs(element[j]));
302          }
303        } else {
304          diagonal_[iColumn] = -value;
305        }
306        addOffset--;
307        work += addOffset;
308      }
309    } else {
310      /* Quadratic*/
311      const int * columnQuadratic = quadratic->getIndices();
312      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
313      const int * columnQuadraticLength = quadratic->getVectorLengths();
314      const double * quadraticElement = quadratic->getElements();
315      for (iColumn=0;iColumn<numberColumns;iColumn++) {
316        CoinWorkDouble value = diagonal[iColumn];
317        CoinBigIndex j;
318        if (CoinAbs(value)>1.0e-100) {
319          value = 1.0/value;
320          for (j=columnQuadraticStart[iColumn];
321               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
322            int jColumn = columnQuadratic[j];
323            if (jColumn>iColumn) {
324              work[jColumn]=-quadraticElement[j];
325            } else if (iColumn==jColumn) {
326              value += quadraticElement[j];
327            }
328          }
329          largest = CoinMax(largest,CoinAbs(value));
330          diagonal_[iColumn] = -value;
331          CoinBigIndex start=columnStart[iColumn];
332          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
333          for (j=start;j<end;j++) {
334            work[row[j]+numberTotal]=element[j];
335            largest = CoinMax(largest,CoinAbs(element[j]));
336          }
337        } else {
338          value = 1.0e100;
339          diagonal_[iColumn] = -value;
340        }
341        addOffset--;
342        work += addOffset;
343      }
344    }
345    /* slacks*/
346    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
347      CoinWorkDouble value = diagonal[iColumn];
348      if (CoinAbs(value)>1.0e-100) {
349        value = 1.0/value;
350        largest = CoinMax(largest,CoinAbs(value));
351      } else {
352        value = 1.0e100;
353      }
354      diagonal_[iColumn] = -value;
355      work[iColumn-numberColumns+numberTotal]=-1.0;
356      addOffset--;
357      work += addOffset;
358    }
359    /* Finish diagonal*/
360    for (iRow=0;iRow<numberRowsModel;iRow++) {
361      diagonal_[iRow+numberTotal]=delta2;
362    }
363    /*check sizes*/
364    largest*=1.0e-20;
365    largest = CoinMin(largest,CHOL_SMALL_VALUE);
366    doubleParameters_[10]=CoinMax(1.0e-20,largest);
367    integerParameters_[20]=0;
368    doubleParameters_[3]=0.0;
369    doubleParameters_[4]=COIN_DBL_MAX;
370    /* Set up LDL cutoff*/
371    integerParameters_[34]=numberTotal;
372    /* KKT*/
373    int * rowsDropped2 = new int[numberRows_];
374    CoinZeroN(rowsDropped2,numberRows_);
375#ifdef CHOL_COMPARE 
376    if (numberRows_<200)
377      factorizePart3(rowsDropped2);
378#endif
379    factorizePart2(rowsDropped2);
380    newDropped=integerParameters_[20];
381    largest=doubleParameters_[3];
382    smallest=doubleParameters_[4];
383    if (model_->messageHandler()->logLevel()>1) 
384      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
385    choleskyCondition_=largest/smallest;
386    /* Should save adjustments in ..R_*/
387    int n1=0,n2=0;
388    CoinWorkDouble * primalR = model_->primalR();
389    CoinWorkDouble * dualR = model_->dualR();
390    for (iRow=0;iRow<numberTotal;iRow++) { 
391      if (rowsDropped2[iRow]) {
392        n1++;
393        /*printf("row region1 %d dropped\n",iRow);*/
394        /*rowsDropped_[iRow]=1;*/
395        rowsDropped_[iRow]=0;
396        primalR[iRow]=doubleParameters_[20];
397      } else {
398        rowsDropped_[iRow]=0;
399        primalR[iRow]=0.0;
400      }
401    }
402    for (;iRow<numberRows_;iRow++) {
403      if (rowsDropped2[iRow]) {
404        n2++;
405        /*printf("row region2 %d dropped\n",iRow);*/
406        /*rowsDropped_[iRow]=1;*/
407        rowsDropped_[iRow]=0;
408        dualR[iRow-numberTotal]=doubleParameters_[34];
409      } else {
410        rowsDropped_[iRow]=0;
411        dualR[iRow-numberTotal]=0.0;
412      }
413    }
414  }
415  return 0;
416}
417/* Factorize - filling in rowsDropped and returning number dropped */
418void 
419ClpCholeskyDense::factorizePart3( int * rowsDropped) 
420{
421  int iColumn;
422  longDouble * xx = sparseFactor_;
423  longDouble * yy = diagonal_;
424  diagonal_ = sparseFactor_ + 40000;
425  sparseFactor_=diagonal_ + numberRows_;
426  /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
427  CoinMemcpyN(xx,40000,sparseFactor_);
428  CoinMemcpyN(yy,numberRows_,diagonal_);
429  int numberDropped=0;
430  CoinWorkDouble largest=0.0;
431  CoinWorkDouble smallest=COIN_DBL_MAX;
432  double dropValue = doubleParameters_[10];
433  int firstPositive=integerParameters_[34];
434  longDouble * work = sparseFactor_;
435  /* Allow for triangular*/
436  int addOffset=numberRows_-1;
437  work--;
438  for (iColumn=0;iColumn<numberRows_;iColumn++) {
439    int iRow;
440    int addOffsetNow = numberRows_-1;;
441    longDouble * workNow=sparseFactor_-1+iColumn;
442    CoinWorkDouble diagonalValue = diagonal_[iColumn];
443    for (iRow=0;iRow<iColumn;iRow++) {
444      double aj = *workNow;
445      addOffsetNow--;
446      workNow += addOffsetNow;
447      diagonalValue -= aj*aj*workDouble_[iRow];
448    }
449    bool dropColumn=false;
450    if (iColumn<firstPositive) {
451      /* must be negative*/
452      if (diagonalValue<=-dropValue) {
453        smallest = CoinMin(smallest,-diagonalValue);
454        largest = CoinMax(largest,-diagonalValue);
455        workDouble_[iColumn]=diagonalValue;
456        diagonalValue = 1.0/diagonalValue;
457      } else {
458        dropColumn=true;
459        workDouble_[iColumn]=-1.0e100;
460        diagonalValue=0.0;
461        integerParameters_[20]++;
462      }
463    } else {
464      /* must be positive*/
465      if (diagonalValue>=dropValue) {
466        smallest = CoinMin(smallest,diagonalValue);
467        largest = CoinMax(largest,diagonalValue);
468        workDouble_[iColumn]=diagonalValue;
469        diagonalValue = 1.0/diagonalValue;
470      } else {
471        dropColumn=true;
472        workDouble_[iColumn]=1.0e100;
473        diagonalValue=0.0;
474        integerParameters_[20]++;
475      }
476    }
477    if (!dropColumn) {
478      diagonal_[iColumn]=diagonalValue;
479      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
480        double value = work[iRow];
481        workNow = sparseFactor_-1; 
482        int addOffsetNow = numberRows_-1;;
483        for (int jColumn=0;jColumn<iColumn;jColumn++) {
484          double aj = workNow[iColumn];
485          double multiplier = workDouble_[jColumn];
486          double ai = workNow[iRow];
487          addOffsetNow--;
488          workNow += addOffsetNow;
489          value -= aj*ai*multiplier;
490        }
491        work[iRow] = value*diagonalValue;
492      }
493    } else {
494      /* drop column*/
495      rowsDropped[iColumn]=2;
496      numberDropped++;
497      diagonal_[iColumn]=0.0;
498      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
499        work[iRow] = 0.0;
500      }
501    }
502    addOffset--;
503    work += addOffset;
504  }
505  doubleParameters_[3]=largest;
506  doubleParameters_[4]=smallest;
507  integerParameters_[20]=numberDropped;
508  sparseFactor_=xx;
509  diagonal_=yy;
510}
511/*#define POS_DEBUG*/
512#ifdef POS_DEBUG
513static int counter=0;
514int ClpCholeskyDense::bNumber(const longDouble * array, int &iRow, int &iCol)
515{
516  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
517  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
518  int k = array-a;
519  assert ((k%BLOCKSQ)==0);
520  iCol=0;
521  int kk=k>>BLOCKSQSHIFT;
522  /*printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);*/
523  k=kk;
524  while (k>=numberBlocks) {
525    iCol++;
526    k -= numberBlocks;
527    numberBlocks--;
528  }
529  iRow = iCol+k;
530  counter++;
531  if(counter>100000)
532    exit(77);
533  return kk;
534 }
535#endif
536/* Factorize - filling in rowsDropped and returning number dropped */
537void 
538ClpCholeskyDense::factorizePart2( int * rowsDropped) 
539{
540  int iColumn;
541  /*longDouble * xx = sparseFactor_;*/
542  /*longDouble * yy = diagonal_;*/
543  /*diagonal_ = sparseFactor_ + 40000;*/
544  /*sparseFactor_=diagonal_ + numberRows_;*/
545  /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
546  /*memcpy(sparseFactor_,xx,40000*sizeof(double));*/
547  /*memcpy(diagonal_,yy,numberRows_*sizeof(double));*/
548  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
549  /* later align on boundary*/
550  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
551  int n=numberRows_;
552  int nRound = numberRows_&(~(BLOCK-1));
553  /* adjust if exact*/
554  if (nRound==n)
555    nRound -= BLOCK;
556  int sizeLastBlock=n-nRound;
557  int get = n*(n-1)/2; /* ? as no diagonal*/
558  int block = numberBlocks*(numberBlocks+1)/2;
559  int ifOdd;
560  int rowLast;
561  if (sizeLastBlock!=BLOCK) {
562    longDouble * aa = &a[(block-1)*BLOCKSQ];
563    rowLast=nRound-1;
564    ifOdd=1;
565    int put = BLOCKSQ;
566    /* do last separately*/
567    put -= (BLOCK-sizeLastBlock)*(BLOCK+1);
568    for (iColumn=numberRows_-1;iColumn>=nRound;iColumn--) {
569      int put2 = put;
570      put -= BLOCK;
571      for (int iRow=numberRows_-1;iRow>iColumn;iRow--) {
572        aa[--put2] = sparseFactor_[--get];
573        assert (aa+put2>=sparseFactor_+get);
574      }
575      /* save diagonal as well*/
576      aa[--put2]=diagonal_[iColumn];
577    }
578    n=nRound;
579    block--;
580  } else {
581    /* exact fit*/
582    rowLast=numberRows_-1;
583    ifOdd=0;
584  }
585  /* Now main loop*/
586  int nBlock=0;
587  for (;n>0;n-=BLOCK) {
588    longDouble * aa = &a[(block-1)*BLOCKSQ];
589    longDouble * aaLast=NULL;
590    int put = BLOCKSQ;
591    int putLast=0;
592    /* see if we have small block*/
593    if (ifOdd) {
594      aaLast = &a[(block-1)*BLOCKSQ];
595      aa = aaLast-BLOCKSQ;
596      putLast = BLOCKSQ-BLOCK+sizeLastBlock;
597    }
598    for (iColumn=n-1;iColumn>=n-BLOCK;iColumn--) {
599      if (aaLast) {
600        /* last bit*/
601        for (int iRow=numberRows_-1;iRow>rowLast;iRow--) {
602          aaLast[--putLast] = sparseFactor_[--get];
603          assert (aaLast+putLast>=sparseFactor_+get);
604        }
605        putLast -= BLOCK-sizeLastBlock;
606      }
607      longDouble * aPut = aa;
608      int j=rowLast;
609      for (int jBlock=0;jBlock<=nBlock;jBlock++) {
610        int put2 = put;
611        int last = CoinMax(j-BLOCK,iColumn);
612        for (int iRow=j;iRow>last;iRow--) {
613          aPut[--put2] = sparseFactor_[--get];
614          assert (aPut+put2>=sparseFactor_+get);
615        }
616        if (j-BLOCK<iColumn) {
617          /* save diagonal as well*/
618          aPut[--put2]=diagonal_[iColumn];
619        }
620        j -= BLOCK;
621        aPut -=BLOCKSQ;
622      }
623      put -= BLOCK;
624    }
625    nBlock++;
626    block -= nBlock+ifOdd;
627  }
628  ClpCholeskyDenseC  info;
629  info.diagonal_=diagonal_;
630  info.doubleParameters_[0]=doubleParameters_[10];
631  info.integerParameters_[0]=integerParameters_[34];
632#ifndef CLP_CILK
633  ClpCholeskyCfactor(&info,a,numberRows_,numberBlocks,
634         diagonal_,workDouble_,rowsDropped);
635#else
636  info.a=a;
637  info.n=numberRows_;
638  info.numberBlocks=numberBlocks;
639  info.work=workDouble_;
640  info.rowsDropped=rowsDropped;
641  info.integerParameters_[1]=model_->numberThreads();
642  ClpCholeskySpawn(&info);
643#endif
644  double largest=0.0;
645  double smallest=COIN_DBL_MAX;
646  int numberDropped=0;
647  for (int i=0;i<numberRows_;i++) {
648    if (diagonal_[i]) {
649      largest = CoinMax(largest,CoinAbs(diagonal_[i]));
650      smallest = CoinMin(smallest,CoinAbs(diagonal_[i]));
651    } else {
652      numberDropped++;
653    }
654  } 
655  doubleParameters_[3]=CoinMax(doubleParameters_[3],1.0/smallest);
656  doubleParameters_[4]=CoinMin(doubleParameters_[4],1.0/largest);
657  integerParameters_[20]+= numberDropped;
658}
659  /* Non leaf recursive factor*/
660void 
661ClpCholeskyCfactor(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, int numberBlocks,
662            longDouble * diagonal, longDouble * work, int * rowsDropped)
663{
664  if (n<=BLOCK) {
665    ClpCholeskyCfactorLeaf(thisStruct, a,n,diagonal,work,rowsDropped);
666  } else {
667    int nb=number_blocks((n+1)>>1);
668    int nThis=number_rows(nb);
669    longDouble * aother;
670    int nLeft=n-nThis;
671    int nintri=(nb*(nb+1))>>1;
672    int nbelow=(numberBlocks-nb)*nb;
673    ClpCholeskyCfactor(thisStruct, a,nThis,numberBlocks,diagonal,work,rowsDropped);
674    ClpCholeskyCtriRec(thisStruct, a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
675    aother=a+number_entries(nintri+nbelow);
676    ClpCholeskyCrecTri(thisStruct, a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
677    ClpCholeskyCfactor(thisStruct, aother,nLeft,
678        numberBlocks-nb,diagonal+nThis,work+nThis,rowsDropped);
679  }
680}
681  /* Non leaf recursive triangle rectangle update*/
682void 
683ClpCholeskyCtriRec(ClpCholeskyDenseC * thisStruct, longDouble * aTri, int nThis, longDouble * aUnder,
684                         longDouble * diagonal, longDouble * work,
685                         int nLeft, int iBlock, int jBlock,
686                         int numberBlocks)
687{
688  if (nThis<=BLOCK&&nLeft<=BLOCK) {
689    ClpCholeskyCtriRecLeaf(thisStruct, aTri,aUnder,diagonal,work,nLeft);
690  } else if (nThis<nLeft) {
691    int nb=number_blocks((nLeft+1)>>1);
692    int nLeft2=number_rows(nb);
693    ClpCholeskyCtriRec(thisStruct, aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
694    ClpCholeskyCtriRec(thisStruct, aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
695          iBlock+nb,jBlock,numberBlocks);
696  } else {
697    int nb=number_blocks((nThis+1)>>1);
698    int nThis2=number_rows(nb);
699    longDouble * aother;
700    int kBlock=jBlock+nb;
701    int i;
702    int nintri=(nb*(nb+1))>>1;
703    int nbelow=(numberBlocks-nb)*nb;
704    ClpCholeskyCtriRec(thisStruct, aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
705    /* and rectangular update */
706    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
707      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
708    aother=aUnder+number_entries(i);
709    ClpCholeskyCrecRec(thisStruct, aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
710          work,kBlock,jBlock,numberBlocks);
711    ClpCholeskyCtriRec(thisStruct, aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
712           work+nThis2,nLeft,
713      iBlock-nb,kBlock-nb,numberBlocks-nb);
714  }
715}
716  /* Non leaf recursive rectangle triangle update*/
717void 
718ClpCholeskyCrecTri(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, int nTri, int nDo,
719                         int iBlock, int jBlock,longDouble * aTri,
720                         longDouble * diagonal, longDouble * work,
721                         int numberBlocks)
722{
723  if (nTri<=BLOCK&&nDo<=BLOCK) {
724    ClpCholeskyCrecTriLeaf(thisStruct, aUnder,aTri,diagonal,work,nTri);
725  } else if (nTri<nDo) {
726    int nb=number_blocks((nDo+1)>>1);
727    int nDo2=number_rows(nb);
728    longDouble * aother;
729    int i;
730    ClpCholeskyCrecTri(thisStruct, aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
731    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
732      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
733    aother=aUnder+number_entries(i);
734    ClpCholeskyCrecTri(thisStruct, aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
735           work+nDo2,numberBlocks-nb);
736  } else {
737    int nb=number_blocks((nTri+1)>>1);
738    int nTri2=number_rows(nb);
739    longDouble * aother;
740    int i;
741    ClpCholeskyCrecTri(thisStruct, aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
742    /* and rectangular update */
743    i=((numberBlocks-iBlock)*(numberBlocks-iBlock+1)-
744      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb+1))>>1;
745    aother=aTri+number_entries(nb);
746    ClpCholeskyCrecRec(thisStruct, aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
747           work,iBlock,jBlock,numberBlocks);
748    ClpCholeskyCrecTri(thisStruct, aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
749         aTri+number_entries(i),diagonal,work,numberBlocks);
750  }
751}
752/* Non leaf recursive rectangle rectangle update,
753   nUnder is number of rows in iBlock,
754   nUnderK is number of rows in kBlock
755*/
756void 
757ClpCholeskyCrecRec(ClpCholeskyDenseC * thisStruct, longDouble * above, int nUnder, int nUnderK,
758                         int nDo, longDouble * aUnder, longDouble *aOther,
759                         longDouble * work,
760                         int iBlock, int jBlock,
761                         int numberBlocks)
762{
763  if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) {
764    assert (nDo==BLOCK&&nUnder==BLOCK);
765    ClpCholeskyCrecRecLeaf(thisStruct, above , aUnder ,  aOther, work, nUnderK);
766  } else if (nDo<=nUnderK&&nUnder<=nUnderK) {
767    int nb=number_blocks((nUnderK+1)>>1);
768    int nUnder2=number_rows(nb);
769    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnder2,nDo,aUnder,aOther,work,
770               iBlock,jBlock,numberBlocks);
771    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
772        aOther+number_entries(nb),work,iBlock,jBlock,numberBlocks);
773  } else if (nUnderK<=nDo&&nUnder<=nDo) {
774    int nb=number_blocks((nDo+1)>>1);
775    int nDo2=number_rows(nb);
776    int i;
777    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnderK,nDo2,aUnder,aOther,work,
778         iBlock,jBlock,numberBlocks);
779    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
780      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
781    ClpCholeskyCrecRec(thisStruct, above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
782           aUnder+number_entries(i),
783           aOther,work+nDo2,
784           iBlock-nb,jBlock,numberBlocks-nb);
785  } else {
786    int nb=number_blocks((nUnder+1)>>1);
787    int nUnder2=number_rows(nb);
788    int i;
789    ClpCholeskyCrecRec(thisStruct, above,nUnder2,nUnderK,nDo,aUnder,aOther,work,
790               iBlock,jBlock,numberBlocks);
791    i=((numberBlocks-iBlock)*(numberBlocks-iBlock-1)-
792      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb-1))>>1;
793    ClpCholeskyCrecRec(thisStruct, above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
794        aOther+number_entries(i),work,iBlock+nb,jBlock,numberBlocks);
795  }
796}
797  /* Leaf recursive factor*/
798void 
799ClpCholeskyCfactorLeaf(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, 
800                longDouble * diagonal, longDouble * work, int * rowsDropped)
801{
802  double dropValue = thisStruct->doubleParameters_[0];
803  int firstPositive=thisStruct->integerParameters_[0];
804  int rowOffset=diagonal-thisStruct->diagonal_;
805  int i, j, k;
806  CoinWorkDouble t00,temp1;
807  longDouble * aa;
808  aa = a-BLOCK;
809  for (j = 0; j < n; j ++) {
810    bool dropColumn;
811    CoinWorkDouble useT00;
812    aa+=BLOCK;
813    t00 = aa[j];
814    for (k = 0; k < j; ++k) {
815      CoinWorkDouble multiplier = work[k];
816      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
817    }
818    dropColumn=false;
819    useT00=t00;
820    if (j+rowOffset<firstPositive) {
821      /* must be negative*/
822      if (t00<=-dropValue) {
823        /*aa[j]=t00;*/
824        t00 = 1.0/t00;
825      } else {
826        dropColumn=true;
827        /*aa[j]=-1.0e100;*/
828        useT00=-1.0e-100;
829        t00=0.0;
830      }
831    } else {
832      /* must be positive*/
833      if (t00>=dropValue) {
834        /*aa[j]=t00;*/
835        t00 = 1.0/t00;
836      } else {
837        dropColumn=true;
838        /*aa[j]=1.0e100;*/
839        useT00=1.0e-100;
840        t00=0.0;
841      }
842    }
843    if (!dropColumn) {
844      diagonal[j]=t00;
845      work[j]=useT00;
846      temp1 = t00;
847      for (i=j+1;i<n;i++) {
848        t00=aa[i];
849        for (k = 0; k < j; ++k) {
850          CoinWorkDouble multiplier = work[k];
851          t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
852        }
853        aa[i]=t00*temp1;
854      }
855    } else {
856      /* drop column*/
857      rowsDropped[j+rowOffset]=2;
858      diagonal[j]=0.0;
859      /*aa[j]=1.0e100;*/
860      work[j]=1.0e100;
861      for (i=j+1;i<n;i++) {
862        aa[i]=0.0;
863      }
864    }
865  }
866}
867  /* Leaf recursive triangle rectangle update*/
868void 
869ClpCholeskyCtriRecLeaf(ClpCholeskyDenseC * thisStruct, longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
870                int nUnder)
871{
872#ifdef POS_DEBUG
873  int iru,icu;
874  int iu=bNumber(aUnder,iru,icu);
875  int irt,ict;
876  int it=bNumber(aTri,irt,ict);
877  /*printf("%d %d\n",iu,it);*/
878  printf("trirecleaf  under (%d,%d), tri (%d,%d)\n",
879         iru,icu,irt,ict);
880  assert (diagonal==thisStruct->diagonal_+ict*BLOCK);
881#endif
882  int j;
883  longDouble * aa;
884#ifdef BLOCKUNROLL
885  if (nUnder==BLOCK) {
886    aa = aTri-2*BLOCK;
887    for (j = 0; j < BLOCK; j +=2) {
888      int i;
889      CoinWorkDouble temp0 = diagonal[j];
890      CoinWorkDouble temp1 = diagonal[j+1];
891      aa+=2*BLOCK;
892      for ( i=0;i<BLOCK;i+=2) {
893        CoinWorkDouble at1;
894        CoinWorkDouble t00=aUnder[i+j*BLOCK];
895        CoinWorkDouble t10=aUnder[i+ BLOCK + j*BLOCK];
896        CoinWorkDouble t01=aUnder[i+1+j*BLOCK];
897        CoinWorkDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
898        int k;
899        for (k = 0; k < j; ++k) {
900          CoinWorkDouble multiplier=work[k];
901          CoinWorkDouble au0 = aUnder[i + k * BLOCK]*multiplier;
902          CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
903          CoinWorkDouble at0 = aTri[j + k * BLOCK];
904          CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK];
905          t00 -= au0 * at0;
906          t10 -= au0 * at1;
907          t01 -= au1 * at0;
908          t11 -= au1 * at1;
909        }
910        t00 *= temp0;
911        at1 = aTri[j + 1 + j*BLOCK]*work[j];
912        t10 -= t00 * at1;
913        t01 *= temp0;
914        t11 -= t01 * at1;
915        aUnder[i+j*BLOCK]=t00;
916        aUnder[i+1+j*BLOCK]=t01;
917        aUnder[i+ BLOCK + j*BLOCK]=t10*temp1;
918        aUnder[i+1+ BLOCK + j*BLOCK]=t11*temp1;
919      }
920    }
921  } else {
922#endif
923    aa = aTri-BLOCK;
924    for (j = 0; j < BLOCK; j ++) {
925      int i;
926      CoinWorkDouble temp1 = diagonal[j];
927      aa+=BLOCK;
928      for (i=0;i<nUnder;i++) {
929        int k;
930        CoinWorkDouble t00=aUnder[i+j*BLOCK];
931        for ( k = 0; k < j; ++k) {
932          CoinWorkDouble multiplier=work[k];
933          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
934        }
935        aUnder[i+j*BLOCK]=t00*temp1;
936      }
937    }
938#ifdef BLOCKUNROLL
939  }
940#endif
941}
942  /* Leaf recursive rectangle triangle update*/
943void ClpCholeskyCrecTriLeaf(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, longDouble * aTri, 
944                                  longDouble * diagonal, longDouble * work, int nUnder)
945{
946#ifdef POS_DEBUG
947  int iru,icu;
948  int iu=bNumber(aUnder,iru,icu);
949  int irt,ict;
950  int it=bNumber(aTri,irt,ict);
951  /*printf("%d %d\n",iu,it);*/
952  printf("rectrileaf  under (%d,%d), tri (%d,%d)\n",
953         iru,icu,irt,ict);
954  assert (diagonal==thisStruct->diagonal_+icu*BLOCK);
955#endif
956  int i, j, k;
957  CoinWorkDouble t00;
958  longDouble * aa;
959#ifdef BLOCKUNROLL
960  if (nUnder==BLOCK) {
961    longDouble * aUnder2 = aUnder-2;
962    aa = aTri-2*BLOCK;
963    for (j = 0; j < BLOCK; j +=2) {
964      CoinWorkDouble t00,t01,t10,t11;
965      aa+=2*BLOCK;
966      aUnder2+=2;
967      t00=aa[j];
968      t01=aa[j+1];
969      t10=aa[j+1+BLOCK];
970      for (k = 0; k < BLOCK; ++k) {
971        CoinWorkDouble multiplier = work[k];
972        CoinWorkDouble a0 = aUnder2[k * BLOCK];
973        CoinWorkDouble a1 = aUnder2[1 + k * BLOCK];
974        CoinWorkDouble x0 = a0 * multiplier;
975        CoinWorkDouble x1 = a1 * multiplier;
976        t00 -= a0 * x0;
977        t01 -= a1 * x0;
978        t10 -= a1 * x1;
979      }
980      aa[j]=t00;
981      aa[j+1]=t01;
982      aa[j+1+BLOCK]=t10;
983      for (i=j+2;i< BLOCK;i+=2) {
984        t00=aa[i];
985        t01=aa[i+BLOCK];
986        t10=aa[i+1];
987        t11=aa[i+1+BLOCK];
988        for (k = 0; k < BLOCK; ++k) {
989          CoinWorkDouble multiplier = work[k];
990          CoinWorkDouble a0 = aUnder2[k * BLOCK]*multiplier;
991          CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
992          t00 -= aUnder[i + k * BLOCK] * a0;
993          t01 -= aUnder[i + k * BLOCK] * a1;
994          t10 -= aUnder[i + 1 + k * BLOCK] * a0;
995          t11 -= aUnder[i + 1 + k * BLOCK] * a1;
996        }
997        aa[i]=t00;
998        aa[i+BLOCK]=t01;
999        aa[i+1]=t10;
1000        aa[i+1+BLOCK]=t11;
1001      }
1002    }
1003  } else {
1004#endif
1005    aa = aTri-BLOCK;
1006    for (j = 0; j < nUnder; j ++) {
1007      aa+=BLOCK;
1008      for (i=j;i< nUnder;i++) {
1009        t00=aa[i];
1010        for (k = 0; k < BLOCK; ++k) {
1011          CoinWorkDouble multiplier = work[k];
1012          t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
1013        }
1014        aa[i]=t00;
1015      }
1016    }
1017#ifdef BLOCKUNROLL
1018  }
1019#endif
1020}
1021/* Leaf recursive rectangle rectangle update,
1022   nUnder is number of rows in iBlock,
1023   nUnderK is number of rows in kBlock
1024*/
1025void 
1026ClpCholeskyCrecRecLeaf(ClpCholeskyDenseC * thisStruct,
1027                                         const longDouble * COIN_RESTRICT above, 
1028                             const longDouble * COIN_RESTRICT aUnder, 
1029                             longDouble * COIN_RESTRICT aOther, 
1030                             const longDouble * COIN_RESTRICT work,
1031                             int nUnder)
1032{
1033#ifdef POS_DEBUG
1034  int ira,ica;
1035  int ia=bNumber(above,ira,ica);
1036  int iru,icu;
1037  int iu=bNumber(aUnder,iru,icu);
1038  int iro,ico;
1039  int io=bNumber(aOther,iro,ico);
1040  /*printf("%d %d %d\n",ia,iu,io);*/
1041  printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n",
1042         ira,ica,iru,icu,iro,ico);
1043#endif
1044  int i, j, k;
1045  longDouble * aa;
1046#ifdef BLOCKUNROLL
1047  aa = aOther-4*BLOCK;
1048  if (nUnder==BLOCK) {
1049    /*#define INTEL*/
1050#ifdef INTEL
1051    aa+=2*BLOCK;
1052    for (j = 0; j < BLOCK; j +=2) {
1053      aa+=2*BLOCK;
1054      for (i=0;i< BLOCK;i+=2) {
1055        CoinWorkDouble t00=aa[i+0*BLOCK];
1056        CoinWorkDouble t10=aa[i+1*BLOCK];
1057        CoinWorkDouble t01=aa[i+1+0*BLOCK];
1058        CoinWorkDouble t11=aa[i+1+1*BLOCK];
1059        for (k=0;k<BLOCK;k++) {
1060          CoinWorkDouble multiplier = work[k];
1061          CoinWorkDouble a00=aUnder[i+k*BLOCK]*multiplier;
1062          CoinWorkDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
1063          t00 -= a00 * above[j + 0 + k * BLOCK];
1064          t10 -= a00 * above[j + 1 + k * BLOCK];
1065          t01 -= a01 * above[j + 0 + k * BLOCK];
1066          t11 -= a01 * above[j + 1 + k * BLOCK];
1067        }
1068        aa[i+0*BLOCK]=t00;
1069        aa[i+1*BLOCK]=t10;
1070        aa[i+1+0*BLOCK]=t01;
1071        aa[i+1+1*BLOCK]=t11;
1072      }
1073    }
1074#else
1075    for (j = 0; j < BLOCK; j +=4) {
1076      aa+=4*BLOCK;
1077      for (i=0;i< BLOCK;i+=4) {
1078        CoinWorkDouble t00=aa[i+0+0*BLOCK];
1079        CoinWorkDouble t10=aa[i+0+1*BLOCK];
1080        CoinWorkDouble t20=aa[i+0+2*BLOCK];
1081        CoinWorkDouble t30=aa[i+0+3*BLOCK];
1082        CoinWorkDouble t01=aa[i+1+0*BLOCK];
1083        CoinWorkDouble t11=aa[i+1+1*BLOCK];
1084        CoinWorkDouble t21=aa[i+1+2*BLOCK];
1085        CoinWorkDouble t31=aa[i+1+3*BLOCK];
1086        CoinWorkDouble t02=aa[i+2+0*BLOCK];
1087        CoinWorkDouble t12=aa[i+2+1*BLOCK];
1088        CoinWorkDouble t22=aa[i+2+2*BLOCK];
1089        CoinWorkDouble t32=aa[i+2+3*BLOCK];
1090        CoinWorkDouble t03=aa[i+3+0*BLOCK];
1091        CoinWorkDouble t13=aa[i+3+1*BLOCK];
1092        CoinWorkDouble t23=aa[i+3+2*BLOCK];
1093        CoinWorkDouble t33=aa[i+3+3*BLOCK];
1094        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
1095        const longDouble * COIN_RESTRICT aboveNow = above+j;
1096        for (k=0;k<BLOCK;k++) {
1097          CoinWorkDouble multiplier = work[k];
1098          CoinWorkDouble a00=aUnderNow[0]*multiplier;
1099          CoinWorkDouble a01=aUnderNow[1]*multiplier;
1100          CoinWorkDouble a02=aUnderNow[2]*multiplier;
1101          CoinWorkDouble a03=aUnderNow[3]*multiplier;
1102          t00 -= a00 * aboveNow[0];
1103          t10 -= a00 * aboveNow[1];
1104          t20 -= a00 * aboveNow[2];
1105          t30 -= a00 * aboveNow[3];
1106          t01 -= a01 * aboveNow[0];
1107          t11 -= a01 * aboveNow[1];
1108          t21 -= a01 * aboveNow[2];
1109          t31 -= a01 * aboveNow[3];
1110          t02 -= a02 * aboveNow[0];
1111          t12 -= a02 * aboveNow[1];
1112          t22 -= a02 * aboveNow[2];
1113          t32 -= a02 * aboveNow[3];
1114          t03 -= a03 * aboveNow[0];
1115          t13 -= a03 * aboveNow[1];
1116          t23 -= a03 * aboveNow[2];
1117          t33 -= a03 * aboveNow[3];
1118          aUnderNow += BLOCK;
1119          aboveNow += BLOCK;
1120        }
1121        aa[i+0+0*BLOCK]=t00;
1122        aa[i+0+1*BLOCK]=t10;
1123        aa[i+0+2*BLOCK]=t20;
1124        aa[i+0+3*BLOCK]=t30;
1125        aa[i+1+0*BLOCK]=t01;
1126        aa[i+1+1*BLOCK]=t11;
1127        aa[i+1+2*BLOCK]=t21;
1128        aa[i+1+3*BLOCK]=t31;
1129        aa[i+2+0*BLOCK]=t02;
1130        aa[i+2+1*BLOCK]=t12;
1131        aa[i+2+2*BLOCK]=t22;
1132        aa[i+2+3*BLOCK]=t32;
1133        aa[i+3+0*BLOCK]=t03;
1134        aa[i+3+1*BLOCK]=t13;
1135        aa[i+3+2*BLOCK]=t23;
1136        aa[i+3+3*BLOCK]=t33;
1137      }
1138    }
1139#endif
1140  } else {
1141    int odd=nUnder&1;
1142    int n=nUnder-odd;
1143    for (j = 0; j < BLOCK; j +=4) {
1144      aa+=4*BLOCK;
1145      for (i=0;i< n;i+=2) {
1146        CoinWorkDouble t00=aa[i+0*BLOCK];
1147        CoinWorkDouble t10=aa[i+1*BLOCK];
1148        CoinWorkDouble t20=aa[i+2*BLOCK];
1149        CoinWorkDouble t30=aa[i+3*BLOCK];
1150        CoinWorkDouble t01=aa[i+1+0*BLOCK];
1151        CoinWorkDouble t11=aa[i+1+1*BLOCK];
1152        CoinWorkDouble t21=aa[i+1+2*BLOCK];
1153        CoinWorkDouble t31=aa[i+1+3*BLOCK];
1154        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
1155        const longDouble * COIN_RESTRICT aboveNow = above+j;
1156        for (k=0;k<BLOCK;k++) {
1157          CoinWorkDouble multiplier = work[k];
1158          CoinWorkDouble a00=aUnderNow[0]*multiplier;
1159          CoinWorkDouble a01=aUnderNow[1]*multiplier;
1160          t00 -= a00 * aboveNow[0];
1161          t10 -= a00 * aboveNow[1];
1162          t20 -= a00 * aboveNow[2];
1163          t30 -= a00 * aboveNow[3];
1164          t01 -= a01 * aboveNow[0];
1165          t11 -= a01 * aboveNow[1];
1166          t21 -= a01 * aboveNow[2];
1167          t31 -= a01 * aboveNow[3];
1168          aUnderNow += BLOCK;
1169          aboveNow += BLOCK;
1170        }
1171        aa[i+0*BLOCK]=t00;
1172        aa[i+1*BLOCK]=t10;
1173        aa[i+2*BLOCK]=t20;
1174        aa[i+3*BLOCK]=t30;
1175        aa[i+1+0*BLOCK]=t01;
1176        aa[i+1+1*BLOCK]=t11;
1177        aa[i+1+2*BLOCK]=t21;
1178        aa[i+1+3*BLOCK]=t31;
1179      }
1180      if (odd) {
1181        CoinWorkDouble t0=aa[n+0*BLOCK];
1182        CoinWorkDouble t1=aa[n+1*BLOCK];
1183        CoinWorkDouble t2=aa[n+2*BLOCK];
1184        CoinWorkDouble t3=aa[n+3*BLOCK];
1185        CoinWorkDouble a0;
1186        for (k=0;k<BLOCK;k++) {
1187          a0=aUnder[n+k*BLOCK]*work[k];
1188          t0 -= a0 * above[j + 0 + k * BLOCK];
1189          t1 -= a0 * above[j + 1 + k * BLOCK];
1190          t2 -= a0 * above[j + 2 + k * BLOCK];
1191          t3 -= a0 * above[j + 3 + k * BLOCK];
1192        }
1193        aa[n+0*BLOCK]=t0;
1194        aa[n+1*BLOCK]=t1;
1195        aa[n+2*BLOCK]=t2;
1196        aa[n+3*BLOCK]=t3;
1197      }
1198    }
1199  }
1200#else
1201  aa = aOther-BLOCK;
1202  for (j = 0; j < BLOCK; j ++) {
1203    aa+=BLOCK;
1204    for (i=0;i< nUnder;i++) {
1205      CoinWorkDouble t00=aa[i+0*BLOCK];
1206      for (k=0;k<BLOCK;k++) {
1207        CoinWorkDouble a00=aUnder[i+k*BLOCK]*work[k];
1208        t00 -= a00 * above[j + k * BLOCK];
1209      }
1210      aa[i]=t00;
1211    }
1212  }
1213#endif
1214}
1215/* Uses factorization to solve. */
1216void 
1217ClpCholeskyDense::solve (CoinWorkDouble * region) 
1218{
1219#ifdef CHOL_COMPARE
1220  double * region2=NULL;
1221  if (numberRows_<200) {
1222    longDouble * xx = sparseFactor_;
1223    longDouble * yy = diagonal_;
1224    diagonal_ = sparseFactor_ + 40000;
1225    sparseFactor_=diagonal_ + numberRows_;
1226    region2 = ClpCopyOfArray(region,numberRows_);
1227    int iRow,iColumn;
1228    int addOffset=numberRows_-1;
1229    longDouble * work = sparseFactor_-1;
1230    for (iColumn=0;iColumn<numberRows_;iColumn++) {
1231      double value = region2[iColumn];
1232      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
1233        region2[iRow] -= value*work[iRow];
1234      addOffset--;
1235      work += addOffset;
1236    }
1237    for (iColumn=numberRows_-1;iColumn>=0;iColumn--) {
1238      double value = region2[iColumn]*diagonal_[iColumn];
1239      work -= addOffset;
1240      addOffset++;
1241      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
1242        value -= region2[iRow]*work[iRow];
1243      region2[iColumn]=value;
1244    }
1245    sparseFactor_=xx;
1246    diagonal_=yy;
1247  }
1248#endif
1249  /*longDouble * xx = sparseFactor_;*/
1250  /*longDouble * yy = diagonal_;*/
1251  /*diagonal_ = sparseFactor_ + 40000;*/
1252  /*sparseFactor_=diagonal_ + numberRows_;*/
1253  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1254  /* later align on boundary*/
1255  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1256  int iBlock;
1257  longDouble * aa = a;
1258  int iColumn;
1259  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1260    int nChunk;
1261    int jBlock;
1262    int iDo = iBlock*BLOCK;
1263    int base=iDo;
1264    if (iDo+BLOCK>numberRows_) {
1265      nChunk=numberRows_-iDo;
1266    } else {
1267      nChunk=BLOCK;
1268    }
1269    solveF1(aa,nChunk,region+iDo);
1270    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1271      base+=BLOCK;
1272      aa+=BLOCKSQ;
1273      if (base+BLOCK>numberRows_) {
1274        nChunk=numberRows_-base;
1275      } else {
1276        nChunk=BLOCK;
1277      }
1278      solveF2(aa,nChunk,region+iDo,region+base);
1279    }
1280    aa+=BLOCKSQ;
1281  }
1282  /* do diagonal outside*/
1283  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1284    region[iColumn] *= diagonal_[iColumn];
1285  int offset=((numberBlocks*(numberBlocks+1))>>1);
1286  aa = a+number_entries(offset-1);
1287  int lBase=(numberBlocks-1)*BLOCK;
1288  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1289    int nChunk;
1290    int jBlock;
1291    int triBase=iBlock*BLOCK;
1292    int iBase=lBase;
1293    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1294      if (iBase+BLOCK>numberRows_) {
1295        nChunk=numberRows_-iBase;
1296      } else {
1297        nChunk=BLOCK;
1298      }
1299      solveB2(aa,nChunk,region+triBase,region+iBase);
1300      iBase-=BLOCK;
1301      aa-=BLOCKSQ;
1302    }
1303    if (triBase+BLOCK>numberRows_) {
1304      nChunk=numberRows_-triBase;
1305    } else {
1306      nChunk=BLOCK;
1307    }
1308    solveB1(aa,nChunk,region+triBase);
1309    aa-=BLOCKSQ;
1310  }
1311#ifdef CHOL_COMPARE
1312  if (numberRows_<200) {
1313    for (int i=0;i<numberRows_;i++) {
1314      assert(CoinAbs(region[i]-region2[i])<1.0e-3);
1315    }
1316    delete [] region2;
1317  }
1318#endif
1319}
1320/* Forward part of solve 1*/
1321void 
1322ClpCholeskyDense::solveF1(longDouble * a,int n,CoinWorkDouble * region)
1323{
1324  int j, k;
1325  CoinWorkDouble t00;
1326  for (j = 0; j < n; j ++) {
1327    t00 = region[j];
1328    for (k = 0; k < j; ++k) {
1329      t00 -= region[k] * a[j + k * BLOCK];
1330    }
1331    /*t00*=a[j + j * BLOCK];*/
1332    region[j] = t00;
1333  }
1334}
1335/* Forward part of solve 2*/
1336void 
1337ClpCholeskyDense::solveF2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1338{
1339  int j, k;
1340#ifdef BLOCKUNROLL
1341  if (n==BLOCK) {
1342    for (k = 0; k < BLOCK; k+=4) {
1343      CoinWorkDouble t0 = region2[0];
1344      CoinWorkDouble t1 = region2[1];
1345      CoinWorkDouble t2 = region2[2];
1346      CoinWorkDouble t3 = region2[3];
1347      t0 -= region[0] * a[0 + 0 * BLOCK];
1348      t1 -= region[0] * a[1 + 0 * BLOCK];
1349      t2 -= region[0] * a[2 + 0 * BLOCK];
1350      t3 -= region[0] * a[3 + 0 * BLOCK];
1351
1352      t0 -= region[1] * a[0 + 1 * BLOCK];
1353      t1 -= region[1] * a[1 + 1 * BLOCK];
1354      t2 -= region[1] * a[2 + 1 * BLOCK];
1355      t3 -= region[1] * a[3 + 1 * BLOCK];
1356
1357      t0 -= region[2] * a[0 + 2 * BLOCK];
1358      t1 -= region[2] * a[1 + 2 * BLOCK];
1359      t2 -= region[2] * a[2 + 2 * BLOCK];
1360      t3 -= region[2] * a[3 + 2 * BLOCK];
1361
1362      t0 -= region[3] * a[0 + 3 * BLOCK];
1363      t1 -= region[3] * a[1 + 3 * BLOCK];
1364      t2 -= region[3] * a[2 + 3 * BLOCK];
1365      t3 -= region[3] * a[3 + 3 * BLOCK];
1366
1367      t0 -= region[4] * a[0 + 4 * BLOCK];
1368      t1 -= region[4] * a[1 + 4 * BLOCK];
1369      t2 -= region[4] * a[2 + 4 * BLOCK];
1370      t3 -= region[4] * a[3 + 4 * BLOCK];
1371
1372      t0 -= region[5] * a[0 + 5 * BLOCK];
1373      t1 -= region[5] * a[1 + 5 * BLOCK];
1374      t2 -= region[5] * a[2 + 5 * BLOCK];
1375      t3 -= region[5] * a[3 + 5 * BLOCK];
1376
1377      t0 -= region[6] * a[0 + 6 * BLOCK];
1378      t1 -= region[6] * a[1 + 6 * BLOCK];
1379      t2 -= region[6] * a[2 + 6 * BLOCK];
1380      t3 -= region[6] * a[3 + 6 * BLOCK];
1381
1382      t0 -= region[7] * a[0 + 7 * BLOCK];
1383      t1 -= region[7] * a[1 + 7 * BLOCK];
1384      t2 -= region[7] * a[2 + 7 * BLOCK];
1385      t3 -= region[7] * a[3 + 7 * BLOCK];
1386#if BLOCK>8
1387      t0 -= region[8] * a[0 + 8 * BLOCK];
1388      t1 -= region[8] * a[1 + 8 * BLOCK];
1389      t2 -= region[8] * a[2 + 8 * BLOCK];
1390      t3 -= region[8] * a[3 + 8 * BLOCK];
1391
1392      t0 -= region[9] * a[0 + 9 * BLOCK];
1393      t1 -= region[9] * a[1 + 9 * BLOCK];
1394      t2 -= region[9] * a[2 + 9 * BLOCK];
1395      t3 -= region[9] * a[3 + 9 * BLOCK];
1396
1397      t0 -= region[10] * a[0 + 10 * BLOCK];
1398      t1 -= region[10] * a[1 + 10 * BLOCK];
1399      t2 -= region[10] * a[2 + 10 * BLOCK];
1400      t3 -= region[10] * a[3 + 10 * BLOCK];
1401
1402      t0 -= region[11] * a[0 + 11 * BLOCK];
1403      t1 -= region[11] * a[1 + 11 * BLOCK];
1404      t2 -= region[11] * a[2 + 11 * BLOCK];
1405      t3 -= region[11] * a[3 + 11 * BLOCK];
1406
1407      t0 -= region[12] * a[0 + 12 * BLOCK];
1408      t1 -= region[12] * a[1 + 12 * BLOCK];
1409      t2 -= region[12] * a[2 + 12 * BLOCK];
1410      t3 -= region[12] * a[3 + 12 * BLOCK];
1411
1412      t0 -= region[13] * a[0 + 13 * BLOCK];
1413      t1 -= region[13] * a[1 + 13 * BLOCK];
1414      t2 -= region[13] * a[2 + 13 * BLOCK];
1415      t3 -= region[13] * a[3 + 13 * BLOCK];
1416
1417      t0 -= region[14] * a[0 + 14 * BLOCK];
1418      t1 -= region[14] * a[1 + 14 * BLOCK];
1419      t2 -= region[14] * a[2 + 14 * BLOCK];
1420      t3 -= region[14] * a[3 + 14 * BLOCK];
1421
1422      t0 -= region[15] * a[0 + 15 * BLOCK];
1423      t1 -= region[15] * a[1 + 15 * BLOCK];
1424      t2 -= region[15] * a[2 + 15 * BLOCK];
1425      t3 -= region[15] * a[3 + 15 * BLOCK];
1426#endif
1427      region2[0] = t0;
1428      region2[1] = t1;
1429      region2[2] = t2;
1430      region2[3] = t3;
1431      region2+=4;
1432      a+=4;
1433    }
1434  } else {
1435#endif
1436    for (k = 0; k < n; ++k) {
1437      CoinWorkDouble t00 = region2[k];
1438      for (j = 0; j < BLOCK; j ++) {
1439        t00 -= region[j] * a[k + j * BLOCK];
1440      }
1441      region2[k] = t00;
1442    }
1443#ifdef BLOCKUNROLL
1444  }
1445#endif
1446}
1447/* Backward part of solve 1*/
1448void 
1449ClpCholeskyDense::solveB1(longDouble * a,int n,CoinWorkDouble * region)
1450{
1451  int j, k;
1452  CoinWorkDouble t00;
1453  for (j = n-1; j >=0; j --) {
1454    t00 = region[j];
1455    for (k = j+1; k < n; ++k) {
1456      t00 -= region[k] * a[k + j * BLOCK];
1457    }
1458    /*t00*=a[j + j * BLOCK];*/
1459    region[j] = t00;
1460  }
1461}
1462/* Backward part of solve 2*/
1463void 
1464ClpCholeskyDense::solveB2(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1465{
1466  int j, k;
1467#ifdef BLOCKUNROLL
1468  if (n==BLOCK) {
1469    for (j = 0; j < BLOCK; j +=4) {
1470      CoinWorkDouble t0 = region[0];
1471      CoinWorkDouble t1 = region[1];
1472      CoinWorkDouble t2 = region[2];
1473      CoinWorkDouble t3 = region[3];
1474      t0 -= region2[0] * a[0 + 0*BLOCK];
1475      t1 -= region2[0] * a[0 + 1*BLOCK];
1476      t2 -= region2[0] * a[0 + 2*BLOCK];
1477      t3 -= region2[0] * a[0 + 3*BLOCK];
1478
1479      t0 -= region2[1] * a[1 + 0*BLOCK];
1480      t1 -= region2[1] * a[1 + 1*BLOCK];
1481      t2 -= region2[1] * a[1 + 2*BLOCK];
1482      t3 -= region2[1] * a[1 + 3*BLOCK];
1483
1484      t0 -= region2[2] * a[2 + 0*BLOCK];
1485      t1 -= region2[2] * a[2 + 1*BLOCK];
1486      t2 -= region2[2] * a[2 + 2*BLOCK];
1487      t3 -= region2[2] * a[2 + 3*BLOCK];
1488
1489      t0 -= region2[3] * a[3 + 0*BLOCK];
1490      t1 -= region2[3] * a[3 + 1*BLOCK];
1491      t2 -= region2[3] * a[3 + 2*BLOCK];
1492      t3 -= region2[3] * a[3 + 3*BLOCK];
1493
1494      t0 -= region2[4] * a[4 + 0*BLOCK];
1495      t1 -= region2[4] * a[4 + 1*BLOCK];
1496      t2 -= region2[4] * a[4 + 2*BLOCK];
1497      t3 -= region2[4] * a[4 + 3*BLOCK];
1498
1499      t0 -= region2[5] * a[5 + 0*BLOCK];
1500      t1 -= region2[5] * a[5 + 1*BLOCK];
1501      t2 -= region2[5] * a[5 + 2*BLOCK];
1502      t3 -= region2[5] * a[5 + 3*BLOCK];
1503
1504      t0 -= region2[6] * a[6 + 0*BLOCK];
1505      t1 -= region2[6] * a[6 + 1*BLOCK];
1506      t2 -= region2[6] * a[6 + 2*BLOCK];
1507      t3 -= region2[6] * a[6 + 3*BLOCK];
1508
1509      t0 -= region2[7] * a[7 + 0*BLOCK];
1510      t1 -= region2[7] * a[7 + 1*BLOCK];
1511      t2 -= region2[7] * a[7 + 2*BLOCK];
1512      t3 -= region2[7] * a[7 + 3*BLOCK];
1513#if BLOCK>8
1514
1515      t0 -= region2[8] * a[8 + 0*BLOCK];
1516      t1 -= region2[8] * a[8 + 1*BLOCK];
1517      t2 -= region2[8] * a[8 + 2*BLOCK];
1518      t3 -= region2[8] * a[8 + 3*BLOCK];
1519
1520      t0 -= region2[9] * a[9 + 0*BLOCK];
1521      t1 -= region2[9] * a[9 + 1*BLOCK];
1522      t2 -= region2[9] * a[9 + 2*BLOCK];
1523      t3 -= region2[9] * a[9 + 3*BLOCK];
1524
1525      t0 -= region2[10] * a[10 + 0*BLOCK];
1526      t1 -= region2[10] * a[10 + 1*BLOCK];
1527      t2 -= region2[10] * a[10 + 2*BLOCK];
1528      t3 -= region2[10] * a[10 + 3*BLOCK];
1529
1530      t0 -= region2[11] * a[11 + 0*BLOCK];
1531      t1 -= region2[11] * a[11 + 1*BLOCK];
1532      t2 -= region2[11] * a[11 + 2*BLOCK];
1533      t3 -= region2[11] * a[11 + 3*BLOCK];
1534
1535      t0 -= region2[12] * a[12 + 0*BLOCK];
1536      t1 -= region2[12] * a[12 + 1*BLOCK];
1537      t2 -= region2[12] * a[12 + 2*BLOCK];
1538      t3 -= region2[12] * a[12 + 3*BLOCK];
1539
1540      t0 -= region2[13] * a[13 + 0*BLOCK];
1541      t1 -= region2[13] * a[13 + 1*BLOCK];
1542      t2 -= region2[13] * a[13 + 2*BLOCK];
1543      t3 -= region2[13] * a[13 + 3*BLOCK];
1544
1545      t0 -= region2[14] * a[14 + 0*BLOCK];
1546      t1 -= region2[14] * a[14 + 1*BLOCK];
1547      t2 -= region2[14] * a[14 + 2*BLOCK];
1548      t3 -= region2[14] * a[14 + 3*BLOCK];
1549
1550      t0 -= region2[15] * a[15 + 0*BLOCK];
1551      t1 -= region2[15] * a[15 + 1*BLOCK];
1552      t2 -= region2[15] * a[15 + 2*BLOCK];
1553      t3 -= region2[15] * a[15 + 3*BLOCK];
1554#endif
1555      region[0] = t0;
1556      region[1] = t1;
1557      region[2] = t2;
1558      region[3] = t3;
1559      a+=4*BLOCK;
1560      region+=4;
1561    }
1562  } else {
1563#endif
1564    for (j = 0; j < BLOCK; j ++) {
1565      CoinWorkDouble t00 = region[j];
1566      for (k = 0; k < n; ++k) {
1567        t00 -= region2[k] * a[k + j * BLOCK];
1568      }
1569      region[j] = t00;
1570    }
1571#ifdef BLOCKUNROLL
1572  }
1573#endif
1574}
Note: See TracBrowser for help on using the repository browser.