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

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

add ids

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 65.6 KB
Line 
1/* $Id: ClpCholeskyDense.cpp 1370 2009-06-04 09:37:13Z 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 double * 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 double * 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 (double * 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,double * 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,double * region, double * 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,double * 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,double * region, double * 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}
1575/* Uses factorization to solve. */
1576void 
1577ClpCholeskyDense::solveLong (CoinWorkDouble * region) 
1578{
1579  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1580  /* later align on boundary*/
1581  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1582  int iBlock;
1583  longDouble * aa = a;
1584  int iColumn;
1585  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1586    int nChunk;
1587    int jBlock;
1588    int iDo = iBlock*BLOCK;
1589    int base=iDo;
1590    if (iDo+BLOCK>numberRows_) {
1591      nChunk=numberRows_-iDo;
1592    } else {
1593      nChunk=BLOCK;
1594    }
1595    solveF1Long(aa,nChunk,region+iDo);
1596    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1597      base+=BLOCK;
1598      aa+=BLOCKSQ;
1599      if (base+BLOCK>numberRows_) {
1600        nChunk=numberRows_-base;
1601      } else {
1602        nChunk=BLOCK;
1603      }
1604      solveF2Long(aa,nChunk,region+iDo,region+base);
1605    }
1606    aa+=BLOCKSQ;
1607  }
1608  /* do diagonal outside*/
1609  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1610    region[iColumn] *= diagonal_[iColumn];
1611  int offset=((numberBlocks*(numberBlocks+1))>>1);
1612  aa = a+number_entries(offset-1);
1613  int lBase=(numberBlocks-1)*BLOCK;
1614  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1615    int nChunk;
1616    int jBlock;
1617    int triBase=iBlock*BLOCK;
1618    int iBase=lBase;
1619    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1620      if (iBase+BLOCK>numberRows_) {
1621        nChunk=numberRows_-iBase;
1622      } else {
1623        nChunk=BLOCK;
1624      }
1625      solveB2Long(aa,nChunk,region+triBase,region+iBase);
1626      iBase-=BLOCK;
1627      aa-=BLOCKSQ;
1628    }
1629    if (triBase+BLOCK>numberRows_) {
1630      nChunk=numberRows_-triBase;
1631    } else {
1632      nChunk=BLOCK;
1633    }
1634    solveB1Long(aa,nChunk,region+triBase);
1635    aa-=BLOCKSQ;
1636  }
1637}
1638/* Forward part of solve 1*/
1639void 
1640ClpCholeskyDense::solveF1Long(longDouble * a,int n,CoinWorkDouble * region)
1641{
1642  int j, k;
1643  CoinWorkDouble t00;
1644  for (j = 0; j < n; j ++) {
1645    t00 = region[j];
1646    for (k = 0; k < j; ++k) {
1647      t00 -= region[k] * a[j + k * BLOCK];
1648    }
1649    /*t00*=a[j + j * BLOCK];*/
1650    region[j] = t00;
1651  }
1652}
1653/* Forward part of solve 2*/
1654void 
1655ClpCholeskyDense::solveF2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1656{
1657  int j, k;
1658#ifdef BLOCKUNROLL
1659  if (n==BLOCK) {
1660    for (k = 0; k < BLOCK; k+=4) {
1661      CoinWorkDouble t0 = region2[0];
1662      CoinWorkDouble t1 = region2[1];
1663      CoinWorkDouble t2 = region2[2];
1664      CoinWorkDouble t3 = region2[3];
1665      t0 -= region[0] * a[0 + 0 * BLOCK];
1666      t1 -= region[0] * a[1 + 0 * BLOCK];
1667      t2 -= region[0] * a[2 + 0 * BLOCK];
1668      t3 -= region[0] * a[3 + 0 * BLOCK];
1669
1670      t0 -= region[1] * a[0 + 1 * BLOCK];
1671      t1 -= region[1] * a[1 + 1 * BLOCK];
1672      t2 -= region[1] * a[2 + 1 * BLOCK];
1673      t3 -= region[1] * a[3 + 1 * BLOCK];
1674
1675      t0 -= region[2] * a[0 + 2 * BLOCK];
1676      t1 -= region[2] * a[1 + 2 * BLOCK];
1677      t2 -= region[2] * a[2 + 2 * BLOCK];
1678      t3 -= region[2] * a[3 + 2 * BLOCK];
1679
1680      t0 -= region[3] * a[0 + 3 * BLOCK];
1681      t1 -= region[3] * a[1 + 3 * BLOCK];
1682      t2 -= region[3] * a[2 + 3 * BLOCK];
1683      t3 -= region[3] * a[3 + 3 * BLOCK];
1684
1685      t0 -= region[4] * a[0 + 4 * BLOCK];
1686      t1 -= region[4] * a[1 + 4 * BLOCK];
1687      t2 -= region[4] * a[2 + 4 * BLOCK];
1688      t3 -= region[4] * a[3 + 4 * BLOCK];
1689
1690      t0 -= region[5] * a[0 + 5 * BLOCK];
1691      t1 -= region[5] * a[1 + 5 * BLOCK];
1692      t2 -= region[5] * a[2 + 5 * BLOCK];
1693      t3 -= region[5] * a[3 + 5 * BLOCK];
1694
1695      t0 -= region[6] * a[0 + 6 * BLOCK];
1696      t1 -= region[6] * a[1 + 6 * BLOCK];
1697      t2 -= region[6] * a[2 + 6 * BLOCK];
1698      t3 -= region[6] * a[3 + 6 * BLOCK];
1699
1700      t0 -= region[7] * a[0 + 7 * BLOCK];
1701      t1 -= region[7] * a[1 + 7 * BLOCK];
1702      t2 -= region[7] * a[2 + 7 * BLOCK];
1703      t3 -= region[7] * a[3 + 7 * BLOCK];
1704#if BLOCK>8
1705      t0 -= region[8] * a[0 + 8 * BLOCK];
1706      t1 -= region[8] * a[1 + 8 * BLOCK];
1707      t2 -= region[8] * a[2 + 8 * BLOCK];
1708      t3 -= region[8] * a[3 + 8 * BLOCK];
1709
1710      t0 -= region[9] * a[0 + 9 * BLOCK];
1711      t1 -= region[9] * a[1 + 9 * BLOCK];
1712      t2 -= region[9] * a[2 + 9 * BLOCK];
1713      t3 -= region[9] * a[3 + 9 * BLOCK];
1714
1715      t0 -= region[10] * a[0 + 10 * BLOCK];
1716      t1 -= region[10] * a[1 + 10 * BLOCK];
1717      t2 -= region[10] * a[2 + 10 * BLOCK];
1718      t3 -= region[10] * a[3 + 10 * BLOCK];
1719
1720      t0 -= region[11] * a[0 + 11 * BLOCK];
1721      t1 -= region[11] * a[1 + 11 * BLOCK];
1722      t2 -= region[11] * a[2 + 11 * BLOCK];
1723      t3 -= region[11] * a[3 + 11 * BLOCK];
1724
1725      t0 -= region[12] * a[0 + 12 * BLOCK];
1726      t1 -= region[12] * a[1 + 12 * BLOCK];
1727      t2 -= region[12] * a[2 + 12 * BLOCK];
1728      t3 -= region[12] * a[3 + 12 * BLOCK];
1729
1730      t0 -= region[13] * a[0 + 13 * BLOCK];
1731      t1 -= region[13] * a[1 + 13 * BLOCK];
1732      t2 -= region[13] * a[2 + 13 * BLOCK];
1733      t3 -= region[13] * a[3 + 13 * BLOCK];
1734
1735      t0 -= region[14] * a[0 + 14 * BLOCK];
1736      t1 -= region[14] * a[1 + 14 * BLOCK];
1737      t2 -= region[14] * a[2 + 14 * BLOCK];
1738      t3 -= region[14] * a[3 + 14 * BLOCK];
1739
1740      t0 -= region[15] * a[0 + 15 * BLOCK];
1741      t1 -= region[15] * a[1 + 15 * BLOCK];
1742      t2 -= region[15] * a[2 + 15 * BLOCK];
1743      t3 -= region[15] * a[3 + 15 * BLOCK];
1744#endif
1745      region2[0] = t0;
1746      region2[1] = t1;
1747      region2[2] = t2;
1748      region2[3] = t3;
1749      region2+=4;
1750      a+=4;
1751    }
1752  } else {
1753#endif
1754    for (k = 0; k < n; ++k) {
1755      CoinWorkDouble t00 = region2[k];
1756      for (j = 0; j < BLOCK; j ++) {
1757        t00 -= region[j] * a[k + j * BLOCK];
1758      }
1759      region2[k] = t00;
1760    }
1761#ifdef BLOCKUNROLL
1762  }
1763#endif
1764}
1765/* Backward part of solve 1*/
1766void 
1767ClpCholeskyDense::solveB1Long(longDouble * a,int n,CoinWorkDouble * region)
1768{
1769  int j, k;
1770  CoinWorkDouble t00;
1771  for (j = n-1; j >=0; j --) {
1772    t00 = region[j];
1773    for (k = j+1; k < n; ++k) {
1774      t00 -= region[k] * a[k + j * BLOCK];
1775    }
1776    /*t00*=a[j + j * BLOCK];*/
1777    region[j] = t00;
1778  }
1779}
1780/* Backward part of solve 2*/
1781void 
1782ClpCholeskyDense::solveB2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1783{
1784  int j, k;
1785#ifdef BLOCKUNROLL
1786  if (n==BLOCK) {
1787    for (j = 0; j < BLOCK; j +=4) {
1788      CoinWorkDouble t0 = region[0];
1789      CoinWorkDouble t1 = region[1];
1790      CoinWorkDouble t2 = region[2];
1791      CoinWorkDouble t3 = region[3];
1792      t0 -= region2[0] * a[0 + 0*BLOCK];
1793      t1 -= region2[0] * a[0 + 1*BLOCK];
1794      t2 -= region2[0] * a[0 + 2*BLOCK];
1795      t3 -= region2[0] * a[0 + 3*BLOCK];
1796
1797      t0 -= region2[1] * a[1 + 0*BLOCK];
1798      t1 -= region2[1] * a[1 + 1*BLOCK];
1799      t2 -= region2[1] * a[1 + 2*BLOCK];
1800      t3 -= region2[1] * a[1 + 3*BLOCK];
1801
1802      t0 -= region2[2] * a[2 + 0*BLOCK];
1803      t1 -= region2[2] * a[2 + 1*BLOCK];
1804      t2 -= region2[2] * a[2 + 2*BLOCK];
1805      t3 -= region2[2] * a[2 + 3*BLOCK];
1806
1807      t0 -= region2[3] * a[3 + 0*BLOCK];
1808      t1 -= region2[3] * a[3 + 1*BLOCK];
1809      t2 -= region2[3] * a[3 + 2*BLOCK];
1810      t3 -= region2[3] * a[3 + 3*BLOCK];
1811
1812      t0 -= region2[4] * a[4 + 0*BLOCK];
1813      t1 -= region2[4] * a[4 + 1*BLOCK];
1814      t2 -= region2[4] * a[4 + 2*BLOCK];
1815      t3 -= region2[4] * a[4 + 3*BLOCK];
1816
1817      t0 -= region2[5] * a[5 + 0*BLOCK];
1818      t1 -= region2[5] * a[5 + 1*BLOCK];
1819      t2 -= region2[5] * a[5 + 2*BLOCK];
1820      t3 -= region2[5] * a[5 + 3*BLOCK];
1821
1822      t0 -= region2[6] * a[6 + 0*BLOCK];
1823      t1 -= region2[6] * a[6 + 1*BLOCK];
1824      t2 -= region2[6] * a[6 + 2*BLOCK];
1825      t3 -= region2[6] * a[6 + 3*BLOCK];
1826
1827      t0 -= region2[7] * a[7 + 0*BLOCK];
1828      t1 -= region2[7] * a[7 + 1*BLOCK];
1829      t2 -= region2[7] * a[7 + 2*BLOCK];
1830      t3 -= region2[7] * a[7 + 3*BLOCK];
1831#if BLOCK>8
1832
1833      t0 -= region2[8] * a[8 + 0*BLOCK];
1834      t1 -= region2[8] * a[8 + 1*BLOCK];
1835      t2 -= region2[8] * a[8 + 2*BLOCK];
1836      t3 -= region2[8] * a[8 + 3*BLOCK];
1837
1838      t0 -= region2[9] * a[9 + 0*BLOCK];
1839      t1 -= region2[9] * a[9 + 1*BLOCK];
1840      t2 -= region2[9] * a[9 + 2*BLOCK];
1841      t3 -= region2[9] * a[9 + 3*BLOCK];
1842
1843      t0 -= region2[10] * a[10 + 0*BLOCK];
1844      t1 -= region2[10] * a[10 + 1*BLOCK];
1845      t2 -= region2[10] * a[10 + 2*BLOCK];
1846      t3 -= region2[10] * a[10 + 3*BLOCK];
1847
1848      t0 -= region2[11] * a[11 + 0*BLOCK];
1849      t1 -= region2[11] * a[11 + 1*BLOCK];
1850      t2 -= region2[11] * a[11 + 2*BLOCK];
1851      t3 -= region2[11] * a[11 + 3*BLOCK];
1852
1853      t0 -= region2[12] * a[12 + 0*BLOCK];
1854      t1 -= region2[12] * a[12 + 1*BLOCK];
1855      t2 -= region2[12] * a[12 + 2*BLOCK];
1856      t3 -= region2[12] * a[12 + 3*BLOCK];
1857
1858      t0 -= region2[13] * a[13 + 0*BLOCK];
1859      t1 -= region2[13] * a[13 + 1*BLOCK];
1860      t2 -= region2[13] * a[13 + 2*BLOCK];
1861      t3 -= region2[13] * a[13 + 3*BLOCK];
1862
1863      t0 -= region2[14] * a[14 + 0*BLOCK];
1864      t1 -= region2[14] * a[14 + 1*BLOCK];
1865      t2 -= region2[14] * a[14 + 2*BLOCK];
1866      t3 -= region2[14] * a[14 + 3*BLOCK];
1867
1868      t0 -= region2[15] * a[15 + 0*BLOCK];
1869      t1 -= region2[15] * a[15 + 1*BLOCK];
1870      t2 -= region2[15] * a[15 + 2*BLOCK];
1871      t3 -= region2[15] * a[15 + 3*BLOCK];
1872#endif
1873      region[0] = t0;
1874      region[1] = t1;
1875      region[2] = t2;
1876      region[3] = t3;
1877      a+=4*BLOCK;
1878      region+=4;
1879    }
1880  } else {
1881#endif
1882    for (j = 0; j < BLOCK; j ++) {
1883      CoinWorkDouble t00 = region[j];
1884      for (k = 0; k < n; ++k) {
1885        t00 -= region2[k] * a[k + j * BLOCK];
1886      }
1887      region[j] = t00;
1888    }
1889#ifdef BLOCKUNROLL
1890  }
1891#endif
1892}
1893/* Uses factorization to solve. */
1894void 
1895ClpCholeskyDense::solveLongWork (CoinWorkDouble * region) 
1896{
1897  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1898  /* later align on boundary*/
1899  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1900  int iBlock;
1901  longDouble * aa = a;
1902  int iColumn;
1903  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1904    int nChunk;
1905    int jBlock;
1906    int iDo = iBlock*BLOCK;
1907    int base=iDo;
1908    if (iDo+BLOCK>numberRows_) {
1909      nChunk=numberRows_-iDo;
1910    } else {
1911      nChunk=BLOCK;
1912    }
1913    solveF1LongWork(aa,nChunk,region+iDo);
1914    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1915      base+=BLOCK;
1916      aa+=BLOCKSQ;
1917      if (base+BLOCK>numberRows_) {
1918        nChunk=numberRows_-base;
1919      } else {
1920        nChunk=BLOCK;
1921      }
1922      solveF2LongWork(aa,nChunk,region+iDo,region+base);
1923    }
1924    aa+=BLOCKSQ;
1925  }
1926  /* do diagonal outside*/
1927  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1928    region[iColumn] *= diagonal_[iColumn];
1929  int offset=((numberBlocks*(numberBlocks+1))>>1);
1930  aa = a+number_entries(offset-1);
1931  int lBase=(numberBlocks-1)*BLOCK;
1932  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1933    int nChunk;
1934    int jBlock;
1935    int triBase=iBlock*BLOCK;
1936    int iBase=lBase;
1937    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1938      if (iBase+BLOCK>numberRows_) {
1939        nChunk=numberRows_-iBase;
1940      } else {
1941        nChunk=BLOCK;
1942      }
1943      solveB2LongWork(aa,nChunk,region+triBase,region+iBase);
1944      iBase-=BLOCK;
1945      aa-=BLOCKSQ;
1946    }
1947    if (triBase+BLOCK>numberRows_) {
1948      nChunk=numberRows_-triBase;
1949    } else {
1950      nChunk=BLOCK;
1951    }
1952    solveB1LongWork(aa,nChunk,region+triBase);
1953    aa-=BLOCKSQ;
1954  }
1955}
1956/* Forward part of solve 1*/
1957void 
1958ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region)
1959{
1960  int j, k;
1961  CoinWorkDouble t00;
1962  for (j = 0; j < n; j ++) {
1963    t00 = region[j];
1964    for (k = 0; k < j; ++k) {
1965      t00 -= region[k] * a[j + k * BLOCK];
1966    }
1967    /*t00*=a[j + j * BLOCK];*/
1968    region[j] = t00;
1969  }
1970}
1971/* Forward part of solve 2*/
1972void 
1973ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1974{
1975  int j, k;
1976#ifdef BLOCKUNROLL
1977  if (n==BLOCK) {
1978    for (k = 0; k < BLOCK; k+=4) {
1979      CoinWorkDouble t0 = region2[0];
1980      CoinWorkDouble t1 = region2[1];
1981      CoinWorkDouble t2 = region2[2];
1982      CoinWorkDouble t3 = region2[3];
1983      t0 -= region[0] * a[0 + 0 * BLOCK];
1984      t1 -= region[0] * a[1 + 0 * BLOCK];
1985      t2 -= region[0] * a[2 + 0 * BLOCK];
1986      t3 -= region[0] * a[3 + 0 * BLOCK];
1987
1988      t0 -= region[1] * a[0 + 1 * BLOCK];
1989      t1 -= region[1] * a[1 + 1 * BLOCK];
1990      t2 -= region[1] * a[2 + 1 * BLOCK];
1991      t3 -= region[1] * a[3 + 1 * BLOCK];
1992
1993      t0 -= region[2] * a[0 + 2 * BLOCK];
1994      t1 -= region[2] * a[1 + 2 * BLOCK];
1995      t2 -= region[2] * a[2 + 2 * BLOCK];
1996      t3 -= region[2] * a[3 + 2 * BLOCK];
1997
1998      t0 -= region[3] * a[0 + 3 * BLOCK];
1999      t1 -= region[3] * a[1 + 3 * BLOCK];
2000      t2 -= region[3] * a[2 + 3 * BLOCK];
2001      t3 -= region[3] * a[3 + 3 * BLOCK];
2002
2003      t0 -= region[4] * a[0 + 4 * BLOCK];
2004      t1 -= region[4] * a[1 + 4 * BLOCK];
2005      t2 -= region[4] * a[2 + 4 * BLOCK];
2006      t3 -= region[4] * a[3 + 4 * BLOCK];
2007
2008      t0 -= region[5] * a[0 + 5 * BLOCK];
2009      t1 -= region[5] * a[1 + 5 * BLOCK];
2010      t2 -= region[5] * a[2 + 5 * BLOCK];
2011      t3 -= region[5] * a[3 + 5 * BLOCK];
2012
2013      t0 -= region[6] * a[0 + 6 * BLOCK];
2014      t1 -= region[6] * a[1 + 6 * BLOCK];
2015      t2 -= region[6] * a[2 + 6 * BLOCK];
2016      t3 -= region[6] * a[3 + 6 * BLOCK];
2017
2018      t0 -= region[7] * a[0 + 7 * BLOCK];
2019      t1 -= region[7] * a[1 + 7 * BLOCK];
2020      t2 -= region[7] * a[2 + 7 * BLOCK];
2021      t3 -= region[7] * a[3 + 7 * BLOCK];
2022#if BLOCK>8
2023      t0 -= region[8] * a[0 + 8 * BLOCK];
2024      t1 -= region[8] * a[1 + 8 * BLOCK];
2025      t2 -= region[8] * a[2 + 8 * BLOCK];
2026      t3 -= region[8] * a[3 + 8 * BLOCK];
2027
2028      t0 -= region[9] * a[0 + 9 * BLOCK];
2029      t1 -= region[9] * a[1 + 9 * BLOCK];
2030      t2 -= region[9] * a[2 + 9 * BLOCK];
2031      t3 -= region[9] * a[3 + 9 * BLOCK];
2032
2033      t0 -= region[10] * a[0 + 10 * BLOCK];
2034      t1 -= region[10] * a[1 + 10 * BLOCK];
2035      t2 -= region[10] * a[2 + 10 * BLOCK];
2036      t3 -= region[10] * a[3 + 10 * BLOCK];
2037
2038      t0 -= region[11] * a[0 + 11 * BLOCK];
2039      t1 -= region[11] * a[1 + 11 * BLOCK];
2040      t2 -= region[11] * a[2 + 11 * BLOCK];
2041      t3 -= region[11] * a[3 + 11 * BLOCK];
2042
2043      t0 -= region[12] * a[0 + 12 * BLOCK];
2044      t1 -= region[12] * a[1 + 12 * BLOCK];
2045      t2 -= region[12] * a[2 + 12 * BLOCK];
2046      t3 -= region[12] * a[3 + 12 * BLOCK];
2047
2048      t0 -= region[13] * a[0 + 13 * BLOCK];
2049      t1 -= region[13] * a[1 + 13 * BLOCK];
2050      t2 -= region[13] * a[2 + 13 * BLOCK];
2051      t3 -= region[13] * a[3 + 13 * BLOCK];
2052
2053      t0 -= region[14] * a[0 + 14 * BLOCK];
2054      t1 -= region[14] * a[1 + 14 * BLOCK];
2055      t2 -= region[14] * a[2 + 14 * BLOCK];
2056      t3 -= region[14] * a[3 + 14 * BLOCK];
2057
2058      t0 -= region[15] * a[0 + 15 * BLOCK];
2059      t1 -= region[15] * a[1 + 15 * BLOCK];
2060      t2 -= region[15] * a[2 + 15 * BLOCK];
2061      t3 -= region[15] * a[3 + 15 * BLOCK];
2062#endif
2063      region2[0] = t0;
2064      region2[1] = t1;
2065      region2[2] = t2;
2066      region2[3] = t3;
2067      region2+=4;
2068      a+=4;
2069    }
2070  } else {
2071#endif
2072    for (k = 0; k < n; ++k) {
2073      CoinWorkDouble t00 = region2[k];
2074      for (j = 0; j < BLOCK; j ++) {
2075        t00 -= region[j] * a[k + j * BLOCK];
2076      }
2077      region2[k] = t00;
2078    }
2079#ifdef BLOCKUNROLL
2080  }
2081#endif
2082}
2083/* Backward part of solve 1*/
2084void 
2085ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region)
2086{
2087  int j, k;
2088  CoinWorkDouble t00;
2089  for (j = n-1; j >=0; j --) {
2090    t00 = region[j];
2091    for (k = j+1; k < n; ++k) {
2092      t00 -= region[k] * a[k + j * BLOCK];
2093    }
2094    /*t00*=a[j + j * BLOCK];*/
2095    region[j] = t00;
2096  }
2097}
2098/* Backward part of solve 2*/
2099void 
2100ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
2101{
2102  int j, k;
2103#ifdef BLOCKUNROLL
2104  if (n==BLOCK) {
2105    for (j = 0; j < BLOCK; j +=4) {
2106      CoinWorkDouble t0 = region[0];
2107      CoinWorkDouble t1 = region[1];
2108      CoinWorkDouble t2 = region[2];
2109      CoinWorkDouble t3 = region[3];
2110      t0 -= region2[0] * a[0 + 0*BLOCK];
2111      t1 -= region2[0] * a[0 + 1*BLOCK];
2112      t2 -= region2[0] * a[0 + 2*BLOCK];
2113      t3 -= region2[0] * a[0 + 3*BLOCK];
2114
2115      t0 -= region2[1] * a[1 + 0*BLOCK];
2116      t1 -= region2[1] * a[1 + 1*BLOCK];
2117      t2 -= region2[1] * a[1 + 2*BLOCK];
2118      t3 -= region2[1] * a[1 + 3*BLOCK];
2119
2120      t0 -= region2[2] * a[2 + 0*BLOCK];
2121      t1 -= region2[2] * a[2 + 1*BLOCK];
2122      t2 -= region2[2] * a[2 + 2*BLOCK];
2123      t3 -= region2[2] * a[2 + 3*BLOCK];
2124
2125      t0 -= region2[3] * a[3 + 0*BLOCK];
2126      t1 -= region2[3] * a[3 + 1*BLOCK];
2127      t2 -= region2[3] * a[3 + 2*BLOCK];
2128      t3 -= region2[3] * a[3 + 3*BLOCK];
2129
2130      t0 -= region2[4] * a[4 + 0*BLOCK];
2131      t1 -= region2[4] * a[4 + 1*BLOCK];
2132      t2 -= region2[4] * a[4 + 2*BLOCK];
2133      t3 -= region2[4] * a[4 + 3*BLOCK];
2134
2135      t0 -= region2[5] * a[5 + 0*BLOCK];
2136      t1 -= region2[5] * a[5 + 1*BLOCK];
2137      t2 -= region2[5] * a[5 + 2*BLOCK];
2138      t3 -= region2[5] * a[5 + 3*BLOCK];
2139
2140      t0 -= region2[6] * a[6 + 0*BLOCK];
2141      t1 -= region2[6] * a[6 + 1*BLOCK];
2142      t2 -= region2[6] * a[6 + 2*BLOCK];
2143      t3 -= region2[6] * a[6 + 3*BLOCK];
2144
2145      t0 -= region2[7] * a[7 + 0*BLOCK];
2146      t1 -= region2[7] * a[7 + 1*BLOCK];
2147      t2 -= region2[7] * a[7 + 2*BLOCK];
2148      t3 -= region2[7] * a[7 + 3*BLOCK];
2149#if BLOCK>8
2150
2151      t0 -= region2[8] * a[8 + 0*BLOCK];
2152      t1 -= region2[8] * a[8 + 1*BLOCK];
2153      t2 -= region2[8] * a[8 + 2*BLOCK];
2154      t3 -= region2[8] * a[8 + 3*BLOCK];
2155
2156      t0 -= region2[9] * a[9 + 0*BLOCK];
2157      t1 -= region2[9] * a[9 + 1*BLOCK];
2158      t2 -= region2[9] * a[9 + 2*BLOCK];
2159      t3 -= region2[9] * a[9 + 3*BLOCK];
2160
2161      t0 -= region2[10] * a[10 + 0*BLOCK];
2162      t1 -= region2[10] * a[10 + 1*BLOCK];
2163      t2 -= region2[10] * a[10 + 2*BLOCK];
2164      t3 -= region2[10] * a[10 + 3*BLOCK];
2165
2166      t0 -= region2[11] * a[11 + 0*BLOCK];
2167      t1 -= region2[11] * a[11 + 1*BLOCK];
2168      t2 -= region2[11] * a[11 + 2*BLOCK];
2169      t3 -= region2[11] * a[11 + 3*BLOCK];
2170
2171      t0 -= region2[12] * a[12 + 0*BLOCK];
2172      t1 -= region2[12] * a[12 + 1*BLOCK];
2173      t2 -= region2[12] * a[12 + 2*BLOCK];
2174      t3 -= region2[12] * a[12 + 3*BLOCK];
2175
2176      t0 -= region2[13] * a[13 + 0*BLOCK];
2177      t1 -= region2[13] * a[13 + 1*BLOCK];
2178      t2 -= region2[13] * a[13 + 2*BLOCK];
2179      t3 -= region2[13] * a[13 + 3*BLOCK];
2180
2181      t0 -= region2[14] * a[14 + 0*BLOCK];
2182      t1 -= region2[14] * a[14 + 1*BLOCK];
2183      t2 -= region2[14] * a[14 + 2*BLOCK];
2184      t3 -= region2[14] * a[14 + 3*BLOCK];
2185
2186      t0 -= region2[15] * a[15 + 0*BLOCK];
2187      t1 -= region2[15] * a[15 + 1*BLOCK];
2188      t2 -= region2[15] * a[15 + 2*BLOCK];
2189      t3 -= region2[15] * a[15 + 3*BLOCK];
2190#endif
2191      region[0] = t0;
2192      region[1] = t1;
2193      region[2] = t2;
2194      region[3] = t3;
2195      a+=4*BLOCK;
2196      region+=4;
2197    }
2198  } else {
2199#endif
2200    for (j = 0; j < BLOCK; j ++) {
2201      CoinWorkDouble t00 = region[j];
2202      for (k = 0; k < n; ++k) {
2203        t00 -= region2[k] * a[k + j * BLOCK];
2204      }
2205      region[j] = t00;
2206    }
2207#ifdef BLOCKUNROLL
2208  }
2209#endif
2210}
Note: See TracBrowser for help on using the repository browser.