source: stable/1.10/Clp/src/ClpCholeskyDense.cpp @ 1404

Last change on this file since 1404 was 1404, checked in by tkr, 11 years ago

Merging r1375:1393 from stable/BSP to stable/1.10
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 63.4 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6#include "CoinPragma.hpp"
7#include "CoinHelperFunctions.hpp"
8
9#include "ClpInterior.hpp"
10#include "ClpCholeskyDense.hpp"
11#include "ClpMessage.hpp"
12#include "ClpQuadraticObjective.hpp"
13
14//#############################################################################
15// Constructors / Destructor / Assignment
16//#############################################################################
17
18//-------------------------------------------------------------------
19// Default Constructor
20//-------------------------------------------------------------------
21ClpCholeskyDense::ClpCholeskyDense () 
22  : ClpCholeskyBase(),
23    borrowSpace_(false)
24{
25  type_=11;;
26}
27
28//-------------------------------------------------------------------
29// Copy constructor
30//-------------------------------------------------------------------
31ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs) 
32  : ClpCholeskyBase(rhs),
33    borrowSpace_(rhs.borrowSpace_)
34{
35  assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
36}
37
38
39//-------------------------------------------------------------------
40// Destructor
41//-------------------------------------------------------------------
42ClpCholeskyDense::~ClpCholeskyDense ()
43{
44  if (borrowSpace_) {
45    // set NULL
46    sparseFactor_=NULL;
47    workDouble_=NULL;
48    diagonal_=NULL;
49  }
50}
51
52//----------------------------------------------------------------
53// Assignment operator
54//-------------------------------------------------------------------
55ClpCholeskyDense &
56ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
57{
58  if (this != &rhs) {
59    assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
60    ClpCholeskyBase::operator=(rhs);
61    borrowSpace_=rhs.borrowSpace_;
62  }
63  return *this;
64}
65//-------------------------------------------------------------------
66// Clone
67//-------------------------------------------------------------------
68ClpCholeskyBase * ClpCholeskyDense::clone() const
69{
70  return new ClpCholeskyDense(*this);
71}
72// If not power of 2 then need to redo a bit
73#define BLOCK 16
74#define BLOCKSHIFT 4
75// Block unroll if power of 2 and at least 8
76#define BLOCKUNROLL
77
78#define BLOCKSQ ( BLOCK*BLOCK )
79#define BLOCKSQSHIFT ( BLOCKSHIFT+BLOCKSHIFT )
80#define number_blocks(x) (((x)+BLOCK-1)>>BLOCKSHIFT)
81#define number_rows(x) ((x)<<BLOCKSHIFT)
82#define number_entries(x) ((x)<<BLOCKSQSHIFT)
83/* Gets space */
84int 
85ClpCholeskyDense::reserveSpace(const ClpCholeskyBase * factor, int numberRows) 
86{
87  numberRows_ = numberRows;
88  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
89  // allow one stripe extra
90  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
91  sizeFactor_=numberBlocks*BLOCKSQ;
92  //#define CHOL_COMPARE
93#ifdef CHOL_COMPARE 
94  sizeFactor_ += 95000;
95#endif
96  if (!factor) {
97    sparseFactor_ = new longDouble [sizeFactor_];
98    rowsDropped_ = new char [numberRows_];
99    memset(rowsDropped_,0,numberRows_);
100    workDouble_ = new longDouble[numberRows_];
101    diagonal_ = new longDouble[numberRows_];
102  } else {
103    borrowSpace_=true;
104    int numberFull = factor->numberRows();
105    sparseFactor_ = factor->sparseFactor()+(factor->size()-sizeFactor_);
106    workDouble_ = factor->workDouble() + (numberFull-numberRows_);
107    diagonal_ = factor->diagonal() + (numberFull-numberRows_);
108  }
109  numberRowsDropped_=0;
110  return 0;
111}
112/* Returns space needed */
113CoinBigIndex
114ClpCholeskyDense::space( int numberRows) const
115{
116 int numberBlocks = (numberRows+BLOCK-1)>>BLOCKSHIFT;
117  // allow one stripe extra
118  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
119  CoinBigIndex sizeFactor=numberBlocks*BLOCKSQ;
120#ifdef CHOL_COMPARE 
121  sizeFactor += 95000;
122#endif
123  return sizeFactor;
124 }
125/* Orders rows and saves pointer to matrix.and model */
126int 
127ClpCholeskyDense::order(ClpInterior * model) 
128{
129  model_=model;
130  int numberRows;
131  int numberRowsModel = model_->numberRows();
132  int numberColumns = model_->numberColumns();
133  if (!doKKT_) {
134    numberRows = numberRowsModel;
135  } else {
136    numberRows = 2*numberRowsModel+numberColumns;
137  }
138  reserveSpace(NULL,numberRows);
139  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
140  return 0;
141}
142/* Does Symbolic factorization given permutation.
143   This is called immediately after order.  If user provides this then
144   user must provide factorize and solve.  Otherwise the default factorization is used
145   returns non-zero if not enough memory */
146int 
147ClpCholeskyDense::symbolic()
148{
149  return 0;
150}
151/* Factorize - filling in rowsDropped and returning number dropped */
152int 
153ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped) 
154{
155  assert (!borrowSpace_);
156  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
157  const int * columnLength = model_->clpMatrix()->getVectorLengths();
158  const int * row = model_->clpMatrix()->getIndices();
159  const double * element = model_->clpMatrix()->getElements();
160  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
161  const int * rowLength = rowCopy_->getVectorLengths();
162  const int * column = rowCopy_->getIndices();
163  const double * elementByRow = rowCopy_->getElements();
164  int numberColumns=model_->clpMatrix()->getNumCols();
165  CoinZeroN(sparseFactor_,sizeFactor_);
166  //perturbation
167  longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
168  perturbation=perturbation*perturbation;
169  if (perturbation>1.0) {
170#ifdef COIN_DEVELOP
171    //if (model_->model()->logLevel()&4)
172      std::cout <<"large perturbation "<<perturbation<<std::endl;
173#endif
174    perturbation=sqrt(perturbation);;
175    perturbation=1.0;
176  }
177  int iRow;
178  int newDropped=0;
179  double largest=1.0;
180  double smallest=COIN_DBL_MAX;
181  longDouble delta2 = model_->delta(); // add delta*delta to diagonal
182  delta2 *= delta2;
183  if (!doKKT_) {
184    longDouble * work = sparseFactor_;
185    work--; // skip diagonal
186    int addOffset=numberRows_-1;
187    const double * diagonalSlack = diagonal + numberColumns;
188    // largest in initial matrix
189    double largest2=1.0e-20;
190    for (iRow=0;iRow<numberRows_;iRow++) {
191      if (!rowsDropped_[iRow]) {
192        CoinBigIndex startRow=rowStart[iRow];
193        CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
194        longDouble diagonalValue = diagonalSlack[iRow]+delta2;
195        for (CoinBigIndex k=startRow;k<endRow;k++) {
196          int iColumn=column[k];
197          CoinBigIndex start=columnStart[iColumn];
198          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
199          longDouble multiplier = diagonal[iColumn]*elementByRow[k];
200          for (CoinBigIndex j=start;j<end;j++) {
201            int jRow=row[j];
202            if (!rowsDropped_[jRow]) {
203              if (jRow>iRow) {
204                longDouble value=element[j]*multiplier;
205                work[jRow] += value;
206              } else if (jRow==iRow) {
207                longDouble value=element[j]*multiplier;
208                diagonalValue += value;
209              }
210            }
211          }
212        } 
213        for (int j=iRow+1;j<numberRows_;j++) 
214          largest2 = CoinMax(largest2,fabs(work[j]));
215        diagonal_[iRow]=diagonalValue;
216        largest2 = CoinMax(largest2,fabs(diagonalValue));
217      } else {
218        // dropped
219        diagonal_[iRow]=1.0;
220      }
221      addOffset--;
222      work += addOffset;
223    }
224    //check sizes
225    largest2*=1.0e-20;
226    largest = CoinMin(largest2,1.0e-11);
227    int numberDroppedBefore=0;
228    for (iRow=0;iRow<numberRows_;iRow++) {
229      int dropped=rowsDropped_[iRow];
230      // Move to int array
231      rowsDropped[iRow]=dropped;
232      if (!dropped) {
233        longDouble diagonal = diagonal_[iRow];
234        if (diagonal>largest2) {
235          diagonal_[iRow]=diagonal+perturbation;
236        } else {
237          diagonal_[iRow]=diagonal+perturbation;
238          rowsDropped[iRow]=2;
239          numberDroppedBefore++;
240        } 
241      }
242    }
243    doubleParameters_[10]=CoinMax(1.0e-20,largest);
244    integerParameters_[20]=0;
245    doubleParameters_[3]=0.0;
246    doubleParameters_[4]=COIN_DBL_MAX;
247    integerParameters_[34]=0; // say all must be positive
248#ifdef CHOL_COMPARE 
249    if (numberRows_<200)
250      factorizePart3(rowsDropped);
251#endif
252    factorizePart2(rowsDropped);
253    newDropped=integerParameters_[20]+numberDroppedBefore;
254    largest=doubleParameters_[3];
255    smallest=doubleParameters_[4];
256    if (model_->messageHandler()->logLevel()>1) 
257      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
258    choleskyCondition_=largest/smallest;
259    //drop fresh makes some formADAT easier
260    if (newDropped||numberRowsDropped_) {
261      newDropped=0;
262      for (int i=0;i<numberRows_;i++) {
263        char dropped = static_cast<char>(rowsDropped[i]);
264        rowsDropped_[i]=dropped;
265        if (dropped==2) {
266          //dropped this time
267          rowsDropped[newDropped++]=i;
268          rowsDropped_[i]=0;
269        } 
270      } 
271      numberRowsDropped_=newDropped;
272      newDropped=-(2+newDropped);
273    } 
274  } else {
275    // KKT
276    CoinPackedMatrix * quadratic = NULL;
277    ClpQuadraticObjective * quadraticObj = 
278      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
279    if (quadraticObj) 
280      quadratic = quadraticObj->quadraticObjective();
281    // matrix
282    int numberRowsModel = model_->numberRows();
283    int numberColumns = model_->numberColumns();
284    int numberTotal = numberColumns + numberRowsModel;
285    longDouble * work = sparseFactor_;
286    work--; // skip diagonal
287    int addOffset=numberRows_-1;
288    int iColumn;
289    if (!quadratic) {
290      for (iColumn=0;iColumn<numberColumns;iColumn++) {
291        longDouble value = diagonal[iColumn];
292        if (fabs(value)>1.0e-100) {
293          value = 1.0/value;
294          largest = CoinMax(largest,fabs(value));
295          diagonal_[iColumn] = -value;
296          CoinBigIndex start=columnStart[iColumn];
297          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
298          for (CoinBigIndex j=start;j<end;j++) {
299            //choleskyRow_[numberElements]=row[j]+numberTotal;
300            //sparseFactor_[numberElements++]=element[j];
301            work[row[j]+numberTotal]=element[j];
302            largest = CoinMax(largest,fabs(element[j]));
303          }
304        } else {
305          diagonal_[iColumn] = -value;
306        }
307        addOffset--;
308        work += addOffset;
309      }
310    } else {
311      // Quadratic
312      const int * columnQuadratic = quadratic->getIndices();
313      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
314      const int * columnQuadraticLength = quadratic->getVectorLengths();
315      const double * quadraticElement = quadratic->getElements();
316      for (iColumn=0;iColumn<numberColumns;iColumn++) {
317        longDouble value = diagonal[iColumn];
318        CoinBigIndex j;
319        if (fabs(value)>1.0e-100) {
320          value = 1.0/value;
321          for (j=columnQuadraticStart[iColumn];
322               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
323            int jColumn = columnQuadratic[j];
324            if (jColumn>iColumn) {
325              work[jColumn]=-quadraticElement[j];
326            } else if (iColumn==jColumn) {
327              value += quadraticElement[j];
328            }
329          }
330          largest = CoinMax(largest,fabs(value));
331          diagonal_[iColumn] = -value;
332          CoinBigIndex start=columnStart[iColumn];
333          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
334          for (j=start;j<end;j++) {
335            work[row[j]+numberTotal]=element[j];
336            largest = CoinMax(largest,fabs(element[j]));
337          }
338        } else {
339          value = 1.0e100;
340          diagonal_[iColumn] = -value;
341        }
342        addOffset--;
343        work += addOffset;
344      }
345    }
346    // slacks
347    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
348      longDouble value = diagonal[iColumn];
349      if (fabs(value)>1.0e-100) {
350        value = 1.0/value;
351        largest = CoinMax(largest,fabs(value));
352      } else {
353        value = 1.0e100;
354      }
355      diagonal_[iColumn] = -value;
356      work[iColumn-numberColumns+numberTotal]=-1.0;
357      addOffset--;
358      work += addOffset;
359    }
360    // Finish diagonal
361    for (iRow=0;iRow<numberRowsModel;iRow++) {
362      diagonal_[iRow+numberTotal]=delta2;
363    }
364    //check sizes
365    largest*=1.0e-20;
366    largest = CoinMin(largest,1.0e-11);
367    doubleParameters_[10]=CoinMax(1.0e-20,largest);
368    integerParameters_[20]=0;
369    doubleParameters_[3]=0.0;
370    doubleParameters_[4]=COIN_DBL_MAX;
371    // Set up LDL cutoff
372    integerParameters_[34]=numberTotal;
373    // KKT
374    int * rowsDropped2 = new int[numberRows_];
375    CoinZeroN(rowsDropped2,numberRows_);
376#ifdef CHOL_COMPARE 
377    if (numberRows_<200)
378      factorizePart3(rowsDropped2);
379#endif
380    factorizePart2(rowsDropped2);
381    newDropped=integerParameters_[20];
382    largest=doubleParameters_[3];
383    smallest=doubleParameters_[4];
384    if (model_->messageHandler()->logLevel()>1) 
385      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
386    choleskyCondition_=largest/smallest;
387    // Should save adjustments in ..R_
388    int n1=0,n2=0;
389    double * primalR = model_->primalR();
390    double * dualR = model_->dualR();
391    for (iRow=0;iRow<numberTotal;iRow++) { 
392      if (rowsDropped2[iRow]) {
393        n1++;
394        //printf("row region1 %d dropped\n",iRow);
395        //rowsDropped_[iRow]=1;
396        rowsDropped_[iRow]=0;
397        primalR[iRow]=doubleParameters_[20];
398      } else {
399        rowsDropped_[iRow]=0;
400        primalR[iRow]=0.0;
401      }
402    }
403    for (;iRow<numberRows_;iRow++) {
404      if (rowsDropped2[iRow]) {
405        n2++;
406        //printf("row region2 %d dropped\n",iRow);
407        //rowsDropped_[iRow]=1;
408        rowsDropped_[iRow]=0;
409        dualR[iRow-numberTotal]=doubleParameters_[34];
410      } else {
411        rowsDropped_[iRow]=0;
412        dualR[iRow-numberTotal]=0.0;
413      }
414    }
415  }
416  return 0;
417}
418/* Factorize - filling in rowsDropped and returning number dropped */
419void 
420ClpCholeskyDense::factorizePart3( int * rowsDropped) 
421{
422  int iColumn;
423  longDouble * xx = sparseFactor_;
424  longDouble * yy = diagonal_;
425  diagonal_ = sparseFactor_ + 40000;
426  sparseFactor_=diagonal_ + numberRows_;
427  //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
428  CoinMemcpyN(xx,40000,sparseFactor_);
429  CoinMemcpyN(yy,numberRows_,diagonal_);
430  int numberDropped=0;
431  double largest=0.0;
432  double smallest=COIN_DBL_MAX;
433  double dropValue = doubleParameters_[10];
434  int firstPositive=integerParameters_[34];
435  longDouble * work = sparseFactor_;
436  // Allow for triangular
437  int addOffset=numberRows_-1;
438  work--;
439  for (iColumn=0;iColumn<numberRows_;iColumn++) {
440    int iRow;
441    int addOffsetNow = numberRows_-1;;
442    longDouble * workNow=sparseFactor_-1+iColumn;
443    double diagonalValue = diagonal_[iColumn];
444    for (iRow=0;iRow<iColumn;iRow++) {
445      double aj = *workNow;
446      addOffsetNow--;
447      workNow += addOffsetNow;
448      double multiplier = workDouble_[iRow];
449      diagonalValue -= aj*aj*multiplier;
450    }
451    bool dropColumn=false;
452    if (iColumn<firstPositive) {
453      // must be negative
454      if (diagonalValue<=-dropValue) {
455        smallest = CoinMin(smallest,-diagonalValue);
456        largest = CoinMax(largest,-diagonalValue);
457        workDouble_[iColumn]=diagonalValue;
458        diagonalValue = 1.0/diagonalValue;
459      } else {
460        dropColumn=true;
461        workDouble_[iColumn]=-1.0e100;
462        diagonalValue=0.0;
463        integerParameters_[20]++;
464      }
465    } else {
466      // must be positive
467      if (diagonalValue>=dropValue) {
468        smallest = CoinMin(smallest,diagonalValue);
469        largest = CoinMax(largest,diagonalValue);
470        workDouble_[iColumn]=diagonalValue;
471        diagonalValue = 1.0/diagonalValue;
472      } else {
473        dropColumn=true;
474        workDouble_[iColumn]=1.0e100;
475        diagonalValue=0.0;
476        integerParameters_[20]++;
477      }
478    }
479    if (!dropColumn) {
480      diagonal_[iColumn]=diagonalValue;
481      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
482        double value = work[iRow];
483        workNow = sparseFactor_-1; 
484        int addOffsetNow = numberRows_-1;;
485        for (int jColumn=0;jColumn<iColumn;jColumn++) {
486          double aj = workNow[iColumn];
487          double multiplier = workDouble_[jColumn];
488          double ai = workNow[iRow];
489          addOffsetNow--;
490          workNow += addOffsetNow;
491          value -= aj*ai*multiplier;
492        }
493        work[iRow] = value*diagonalValue;
494      }
495    } else {
496      // drop column
497      rowsDropped[iColumn]=2;
498      numberDropped++;
499      diagonal_[iColumn]=0.0;
500      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
501        work[iRow] = 0.0;
502      }
503    }
504    addOffset--;
505    work += addOffset;
506  }
507  doubleParameters_[3]=largest;
508  doubleParameters_[4]=smallest;
509  integerParameters_[20]=numberDropped;
510  sparseFactor_=xx;
511  diagonal_=yy;
512}
513//#define POS_DEBUG
514#ifdef POS_DEBUG
515static int counter=0;
516int ClpCholeskyDense::bNumber(const longDouble * array, int &iRow, int &iCol)
517{
518  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
519  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
520  int k = array-a;
521  assert ((k%BLOCKSQ)==0);
522  iCol=0;
523  int kk=k>>BLOCKSQSHIFT;
524  //printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);
525  k=kk;
526  while (k>=numberBlocks) {
527    iCol++;
528    k -= numberBlocks;
529    numberBlocks--;
530  }
531  iRow = iCol+k;
532  counter++;
533  if(counter>100000)
534    exit(77);
535  return kk;
536 }
537#endif
538/* Factorize - filling in rowsDropped and returning number dropped */
539void 
540ClpCholeskyDense::factorizePart2( int * rowsDropped) 
541{
542  int iColumn;
543  //longDouble * xx = sparseFactor_;
544  //longDouble * yy = diagonal_;
545  //diagonal_ = sparseFactor_ + 40000;
546  //sparseFactor_=diagonal_ + numberRows_;
547  //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
548  //memcpy(sparseFactor_,xx,40000*sizeof(double));
549  //memcpy(diagonal_,yy,numberRows_*sizeof(double));
550  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
551  // later align on boundary
552  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
553  int n=numberRows_;
554  int nRound = numberRows_&(~(BLOCK-1));
555  // adjust if exact
556  if (nRound==n)
557    nRound -= BLOCK;
558  int sizeLastBlock=n-nRound;
559  int get = n*(n-1)/2; // ? as no diagonal
560  int block = numberBlocks*(numberBlocks+1)/2;
561  int ifOdd;
562  int rowLast;
563  if (sizeLastBlock!=BLOCK) {
564    longDouble * aa = &a[(block-1)*BLOCKSQ];
565    rowLast=nRound-1;
566    ifOdd=1;
567    int put = BLOCKSQ;
568    // do last separately
569    put -= (BLOCK-sizeLastBlock)*(BLOCK+1);
570    for (iColumn=numberRows_-1;iColumn>=nRound;iColumn--) {
571      int put2 = put;
572      put -= BLOCK;
573      for (int iRow=numberRows_-1;iRow>iColumn;iRow--) {
574        aa[--put2] = sparseFactor_[--get];
575        assert (aa+put2>=sparseFactor_+get);
576      }
577      // save diagonal as well
578      aa[--put2]=diagonal_[iColumn];
579    }
580    n=nRound;
581    block--;
582  } else {
583    // exact fit
584    rowLast=numberRows_-1;
585    ifOdd=0;
586  }
587  // Now main loop
588  int nBlock=0;
589  for (;n>0;n-=BLOCK) {
590    longDouble * aa = &a[(block-1)*BLOCKSQ];
591    longDouble * aaLast=NULL;
592    int put = BLOCKSQ;
593    int putLast=0;
594    // see if we have small block
595    if (ifOdd) {
596      aaLast = &a[(block-1)*BLOCKSQ];
597      aa = aaLast-BLOCKSQ;
598      putLast = BLOCKSQ-BLOCK+sizeLastBlock;
599    }
600    for (iColumn=n-1;iColumn>=n-BLOCK;iColumn--) {
601      if (aaLast) {
602        // last bit
603        for (int iRow=numberRows_-1;iRow>rowLast;iRow--) {
604          aaLast[--putLast] = sparseFactor_[--get];
605          assert (aaLast+putLast>=sparseFactor_+get);
606        }
607        putLast -= BLOCK-sizeLastBlock;
608      }
609      longDouble * aPut = aa;
610      int j=rowLast;
611      for (int jBlock=0;jBlock<=nBlock;jBlock++) {
612        int put2 = put;
613        int last = CoinMax(j-BLOCK,iColumn);
614        for (int iRow=j;iRow>last;iRow--) {
615          aPut[--put2] = sparseFactor_[--get];
616          assert (aPut+put2>=sparseFactor_+get);
617        }
618        if (j-BLOCK<iColumn) {
619          // save diagonal as well
620          aPut[--put2]=diagonal_[iColumn];
621        }
622        j -= BLOCK;
623        aPut -=BLOCKSQ;
624      }
625      put -= BLOCK;
626    }
627    nBlock++;
628    block -= nBlock+ifOdd;
629  }
630  factor(a,numberRows_,numberBlocks,
631         diagonal_,workDouble_,rowsDropped);
632  //sparseFactor_=xx;
633  //diagonal_=yy;
634}
635// Non leaf recursive factor
636void 
637ClpCholeskyDense::factor(longDouble * a, int n, int numberBlocks,
638            longDouble * diagonal, longDouble * work, int * rowsDropped)
639{
640  if (n<=BLOCK) {
641    factorLeaf(a,n,diagonal,work,rowsDropped);
642  } else {
643    int nb=number_blocks((n+1)>>1);
644    int nThis=number_rows(nb);
645    longDouble * aother;
646    int nLeft=n-nThis;
647    int nintri=(nb*(nb+1))>>1;
648    int nbelow=(numberBlocks-nb)*nb;
649    factor(a,nThis,numberBlocks,diagonal,work,rowsDropped);
650    triRec(a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
651    aother=a+number_entries(nintri+nbelow);
652    recTri(a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
653    factor(aother,nLeft,
654        numberBlocks-nb,diagonal+nThis,work+nThis,rowsDropped);
655  }
656}
657// Non leaf recursive triangle rectangle update
658void 
659ClpCholeskyDense::triRec(longDouble * aTri, int nThis, longDouble * aUnder,
660                         longDouble * diagonal, longDouble * work,
661                         int nLeft, int iBlock, int jBlock,
662                         int numberBlocks)
663{
664  if (nThis<=BLOCK&&nLeft<=BLOCK) {
665    triRecLeaf(aTri,aUnder,diagonal,work,nLeft);
666  } else if (nThis<nLeft) {
667    int nb=number_blocks((nLeft+1)>>1);
668    int nLeft2=number_rows(nb);
669    triRec(aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
670    triRec(aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
671          iBlock+nb,jBlock,numberBlocks);
672  } else {
673    int nb=number_blocks((nThis+1)>>1);
674    int nThis2=number_rows(nb);
675    longDouble * aother;
676    int kBlock=jBlock+nb;
677    int i;
678    int nintri=(nb*(nb+1))>>1;
679    int nbelow=(numberBlocks-nb)*nb;
680    triRec(aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
681    /* and rectangular update */
682    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
683      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
684    aother=aUnder+number_entries(i);
685    recRec(aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
686          work,kBlock,jBlock,numberBlocks);
687    triRec(aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
688           work+nThis2,nLeft,
689      iBlock-nb,kBlock-nb,numberBlocks-nb);
690  }
691}
692// Non leaf recursive rectangle triangle update
693void 
694ClpCholeskyDense::recTri(longDouble * aUnder, int nTri, int nDo,
695                         int iBlock, int jBlock,longDouble * aTri,
696                         longDouble * diagonal, longDouble * work,
697                         int numberBlocks)
698{
699  if (nTri<=BLOCK&&nDo<=BLOCK) {
700    recTriLeaf(aUnder,aTri,diagonal,work,nTri);
701  } else if (nTri<nDo) {
702    int nb=number_blocks((nDo+1)>>1);
703    int nDo2=number_rows(nb);
704    longDouble * aother;
705    int i;
706    recTri(aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
707    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
708      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
709    aother=aUnder+number_entries(i);
710    recTri(aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
711           work+nDo2,numberBlocks-nb);
712  } else {
713    int nb=number_blocks((nTri+1)>>1);
714    int nTri2=number_rows(nb);
715    longDouble * aother;
716    int i;
717    recTri(aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
718    /* and rectangular update */
719    i=((numberBlocks-iBlock)*(numberBlocks-iBlock+1)-
720      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb+1))>>1;
721    aother=aTri+number_entries(nb);
722    recRec(aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
723           work,iBlock,jBlock,numberBlocks);
724    recTri(aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
725         aTri+number_entries(i),diagonal,work,numberBlocks);
726  }
727}
728/* Non leaf recursive rectangle rectangle update,
729   nUnder is number of rows in iBlock,
730   nUnderK is number of rows in kBlock
731*/
732void 
733ClpCholeskyDense::recRec(longDouble * above, int nUnder, int nUnderK,
734                         int nDo, longDouble * aUnder, longDouble *aOther,
735                         longDouble * work,
736                         int iBlock, int jBlock,
737                         int numberBlocks)
738{
739  if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) {
740    assert (nDo==BLOCK&&nUnder==BLOCK);
741    recRecLeaf(above , aUnder ,  aOther, work, nUnderK);
742  } else if (nDo<=nUnderK&&nUnder<=nUnderK) {
743    int nb=number_blocks((nUnderK+1)>>1);
744    int nUnder2=number_rows(nb);
745    recRec(above,nUnder,nUnder2,nDo,aUnder,aOther,work,
746               iBlock,jBlock,numberBlocks);
747    recRec(above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
748        aOther+number_entries(nb),work,iBlock,jBlock,numberBlocks);
749  } else if (nUnderK<=nDo&&nUnder<=nDo) {
750    int nb=number_blocks((nDo+1)>>1);
751    int nDo2=number_rows(nb);
752    int i;
753    recRec(above,nUnder,nUnderK,nDo2,aUnder,aOther,work,
754         iBlock,jBlock,numberBlocks);
755    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
756      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
757    recRec(above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
758           aUnder+number_entries(i),
759           aOther,work+nDo2,
760           iBlock-nb,jBlock,numberBlocks-nb);
761  } else {
762    int nb=number_blocks((nUnder+1)>>1);
763    int nUnder2=number_rows(nb);
764    int i;
765    recRec(above,nUnder2,nUnderK,nDo,aUnder,aOther,work,
766               iBlock,jBlock,numberBlocks);
767    i=((numberBlocks-iBlock)*(numberBlocks-iBlock-1)-
768      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb-1))>>1;
769    recRec(above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
770        aOther+number_entries(i),work,iBlock+nb,jBlock,numberBlocks);
771  }
772}
773// Leaf recursive factor
774void 
775ClpCholeskyDense::factorLeaf(longDouble * a, int n, 
776                longDouble * diagonal, longDouble * work, int * rowsDropped)
777{
778  longDouble largest=doubleParameters_[3];
779  longDouble smallest=doubleParameters_[4];
780  double dropValue = doubleParameters_[10];
781  int firstPositive=integerParameters_[34];
782  size_t rowOffset=diagonal-diagonal_;
783  int numberDropped=0;
784  int i, j, k;
785  longDouble t00,temp1;
786  longDouble * aa;
787  aa = a-BLOCK;
788  for (j = 0; j < n; j ++) {
789    aa+=BLOCK;
790    t00 = aa[j];
791    for (k = 0; k < j; ++k) {
792      longDouble multiplier = work[k];
793      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
794    }
795    bool dropColumn=false;
796    longDouble useT00=t00;
797    if ((int)(j+rowOffset)<firstPositive) {
798      // must be negative
799      if (t00<=-dropValue) {
800        smallest = CoinMin(smallest,-t00);
801        largest = CoinMax(largest,-t00);
802        //aa[j]=t00;
803        t00 = 1.0/t00;
804      } else {
805        dropColumn=true;
806        //aa[j]=-1.0e100;
807        useT00=-1.0e-100;
808        t00=0.0;
809        integerParameters_[20]++;
810      }
811    } else {
812      // must be positive
813      if (t00>=dropValue) {
814        smallest = CoinMin(smallest,t00);
815        largest = CoinMax(largest,t00);
816        //aa[j]=t00;
817        t00 = 1.0/t00;
818      } else {
819        dropColumn=true;
820        //aa[j]=1.0e100;
821        useT00=1.0e-100;
822        t00=0.0;
823        integerParameters_[20]++;
824      }
825    }
826    if (!dropColumn) {
827      diagonal[j]=t00;
828      work[j]=useT00;
829      temp1 = t00;
830      for (i=j+1;i<n;i++) {
831        t00=aa[i];
832        for (k = 0; k < j; ++k) {
833          longDouble multiplier = work[k];
834          t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
835        }
836        aa[i]=t00*temp1;
837      }
838    } else {
839      // drop column
840      rowsDropped[j+rowOffset]=2;
841      numberDropped++;
842      diagonal[j]=0.0;
843      //aa[j]=1.0e100;
844      work[j]=1.0e100;
845      for (i=j+1;i<n;i++) {
846        aa[i]=0.0;
847      }
848    }
849  }
850  doubleParameters_[3]=largest;
851  doubleParameters_[4]=smallest;
852  integerParameters_[20]+=numberDropped;
853}
854// Leaf recursive triangle rectangle update
855void 
856ClpCholeskyDense::triRecLeaf(longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
857                int nUnder)
858{
859#ifdef POS_DEBUG
860  int iru,icu;
861  int iu=bNumber(aUnder,iru,icu);
862  int irt,ict;
863  int it=bNumber(aTri,irt,ict);
864  //printf("%d %d\n",iu,it);
865  printf("trirecleaf  under (%d,%d), tri (%d,%d)\n",
866         iru,icu,irt,ict);
867  assert (diagonal==diagonal_+ict*BLOCK);
868#endif
869  int j;
870  longDouble * aa;
871#ifdef BLOCKUNROLL
872  if (nUnder==BLOCK) {
873    aa = aTri-2*BLOCK;
874    for (j = 0; j < BLOCK; j +=2) {
875      aa+=2*BLOCK;
876      longDouble temp0 = diagonal[j];
877      longDouble temp1 = diagonal[j+1];
878      for (int i=0;i<BLOCK;i+=2) {
879        longDouble t00=aUnder[i+j*BLOCK];
880        longDouble t10=aUnder[i+ BLOCK + j*BLOCK];
881        longDouble t01=aUnder[i+1+j*BLOCK];
882        longDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
883        for (int k = 0; k < j; ++k) {
884          longDouble multiplier=work[k];
885          longDouble au0 = aUnder[i + k * BLOCK]*multiplier;
886          longDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
887          longDouble at0 = aTri[j + k * BLOCK];
888          longDouble at1 = aTri[j + 1 + k * BLOCK];
889          t00 -= au0 * at0;
890          t10 -= au0 * at1;
891          t01 -= au1 * at0;
892          t11 -= au1 * at1;
893        }
894        t00 *= temp0;
895        longDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
896        t10 -= t00 * at1;
897        t01 *= temp0;
898        t11 -= t01 * at1;
899        aUnder[i+j*BLOCK]=t00;
900        aUnder[i+1+j*BLOCK]=t01;
901        aUnder[i+ BLOCK + j*BLOCK]=t10*temp1;
902        aUnder[i+1+ BLOCK + j*BLOCK]=t11*temp1;
903      }
904    }
905  } else {
906#endif
907    aa = aTri-BLOCK;
908    for (j = 0; j < BLOCK; j ++) {
909      aa+=BLOCK;
910      longDouble temp1 = diagonal[j];
911      for (int i=0;i<nUnder;i++) {
912        longDouble t00=aUnder[i+j*BLOCK];
913        for (int k = 0; k < j; ++k) {
914          longDouble multiplier=work[k];
915          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
916        }
917        aUnder[i+j*BLOCK]=t00*temp1;
918      }
919    }
920#ifdef BLOCKUNROLL
921  }
922#endif
923}
924// Leaf recursive rectangle triangle update
925void ClpCholeskyDense::recTriLeaf(longDouble * aUnder, longDouble * aTri, 
926                                  longDouble * diagonal, longDouble * work, int nUnder)
927{
928#ifdef POS_DEBUG
929  int iru,icu;
930  int iu=bNumber(aUnder,iru,icu);
931  int irt,ict;
932  int it=bNumber(aTri,irt,ict);
933  //printf("%d %d\n",iu,it);
934  printf("rectrileaf  under (%d,%d), tri (%d,%d)\n",
935         iru,icu,irt,ict);
936  assert (diagonal==diagonal_+icu*BLOCK);
937#endif
938  int i, j, k;
939  longDouble t00;
940  longDouble * aa;
941#ifdef BLOCKUNROLL
942  if (nUnder==BLOCK) {
943    longDouble * aUnder2 = aUnder-2;
944    aa = aTri-2*BLOCK;
945    for (j = 0; j < BLOCK; j +=2) {
946      longDouble t00,t01,t10,t11;
947      aa+=2*BLOCK;
948      aUnder2+=2;
949      t00=aa[j];
950      t01=aa[j+1];
951      t10=aa[j+1+BLOCK];
952      for (k = 0; k < BLOCK; ++k) {
953        longDouble multiplier = work[k];
954        longDouble a0 = aUnder2[k * BLOCK];
955        longDouble a1 = aUnder2[1 + k * BLOCK];
956        longDouble x0 = a0 * multiplier;
957        longDouble x1 = a1 * multiplier;
958        t00 -= a0 * x0;
959        t01 -= a1 * x0;
960        t10 -= a1 * x1;
961      }
962      aa[j]=t00;
963      aa[j+1]=t01;
964      aa[j+1+BLOCK]=t10;
965      for (i=j+2;i< BLOCK;i+=2) {
966        t00=aa[i];
967        t01=aa[i+BLOCK];
968        t10=aa[i+1];
969        t11=aa[i+1+BLOCK];
970        for (k = 0; k < BLOCK; ++k) {
971          longDouble multiplier = work[k];
972          longDouble a0 = aUnder2[k * BLOCK]*multiplier;
973          longDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
974          t00 -= aUnder[i + k * BLOCK] * a0;
975          t01 -= aUnder[i + k * BLOCK] * a1;
976          t10 -= aUnder[i + 1 + k * BLOCK] * a0;
977          t11 -= aUnder[i + 1 + k * BLOCK] * a1;
978        }
979        aa[i]=t00;
980        aa[i+BLOCK]=t01;
981        aa[i+1]=t10;
982        aa[i+1+BLOCK]=t11;
983      }
984    }
985  } else {
986#endif
987    aa = aTri-BLOCK;
988    for (j = 0; j < nUnder; j ++) {
989      aa+=BLOCK;
990      for (i=j;i< nUnder;i++) {
991        t00=aa[i];
992        for (k = 0; k < BLOCK; ++k) {
993          longDouble multiplier = work[k];
994          t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
995        }
996        aa[i]=t00;
997      }
998    }
999#ifdef BLOCKUNROLL
1000  }
1001#endif
1002}
1003/* Leaf recursive rectangle rectangle update,
1004   nUnder is number of rows in iBlock,
1005   nUnderK is number of rows in kBlock
1006*/
1007void 
1008ClpCholeskyDense::recRecLeaf(const longDouble * COIN_RESTRICT above, 
1009                             const longDouble * COIN_RESTRICT aUnder, 
1010                             longDouble * COIN_RESTRICT aOther, 
1011                             const longDouble * COIN_RESTRICT work,
1012                             int nUnder)
1013{
1014#ifdef POS_DEBUG
1015  int ira,ica;
1016  int ia=bNumber(above,ira,ica);
1017  int iru,icu;
1018  int iu=bNumber(aUnder,iru,icu);
1019  int iro,ico;
1020  int io=bNumber(aOther,iro,ico);
1021  //printf("%d %d %d\n",ia,iu,io);
1022  printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n",
1023         ira,ica,iru,icu,iro,ico);
1024#endif
1025  int i, j, k;
1026  longDouble * aa;
1027#ifdef BLOCKUNROLL
1028  aa = aOther-4*BLOCK;
1029  if (nUnder==BLOCK) {
1030    //#define INTEL
1031#ifdef INTEL
1032    aa+=2*BLOCK;
1033    for (j = 0; j < BLOCK; j +=2) {
1034      aa+=2*BLOCK;
1035      for (i=0;i< BLOCK;i+=2) {
1036        longDouble t00=aa[i+0*BLOCK];
1037        longDouble t10=aa[i+1*BLOCK];
1038        longDouble t01=aa[i+1+0*BLOCK];
1039        longDouble t11=aa[i+1+1*BLOCK];
1040        for (k=0;k<BLOCK;k++) {
1041          longDouble multiplier = work[k];
1042          longDouble a00=aUnder[i+k*BLOCK]*multiplier;
1043          longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
1044          t00 -= a00 * above[j + 0 + k * BLOCK];
1045          t10 -= a00 * above[j + 1 + k * BLOCK];
1046          t01 -= a01 * above[j + 0 + k * BLOCK];
1047          t11 -= a01 * above[j + 1 + k * BLOCK];
1048        }
1049        aa[i+0*BLOCK]=t00;
1050        aa[i+1*BLOCK]=t10;
1051        aa[i+1+0*BLOCK]=t01;
1052        aa[i+1+1*BLOCK]=t11;
1053      }
1054    }
1055#else
1056    for (j = 0; j < BLOCK; j +=4) {
1057      aa+=4*BLOCK;
1058      for (i=0;i< BLOCK;i+=4) {
1059        longDouble t00=aa[i+0+0*BLOCK];
1060        longDouble t10=aa[i+0+1*BLOCK];
1061        longDouble t20=aa[i+0+2*BLOCK];
1062        longDouble t30=aa[i+0+3*BLOCK];
1063        longDouble t01=aa[i+1+0*BLOCK];
1064        longDouble t11=aa[i+1+1*BLOCK];
1065        longDouble t21=aa[i+1+2*BLOCK];
1066        longDouble t31=aa[i+1+3*BLOCK];
1067        longDouble t02=aa[i+2+0*BLOCK];
1068        longDouble t12=aa[i+2+1*BLOCK];
1069        longDouble t22=aa[i+2+2*BLOCK];
1070        longDouble t32=aa[i+2+3*BLOCK];
1071        longDouble t03=aa[i+3+0*BLOCK];
1072        longDouble t13=aa[i+3+1*BLOCK];
1073        longDouble t23=aa[i+3+2*BLOCK];
1074        longDouble t33=aa[i+3+3*BLOCK];
1075        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
1076        const longDouble * COIN_RESTRICT aboveNow = above+j;
1077        for (k=0;k<BLOCK;k++) {
1078          longDouble multiplier = work[k];
1079          longDouble a00=aUnderNow[0]*multiplier;
1080          longDouble a01=aUnderNow[1]*multiplier;
1081          longDouble a02=aUnderNow[2]*multiplier;
1082          longDouble a03=aUnderNow[3]*multiplier;
1083          t00 -= a00 * aboveNow[0];
1084          t10 -= a00 * aboveNow[1];
1085          t20 -= a00 * aboveNow[2];
1086          t30 -= a00 * aboveNow[3];
1087          t01 -= a01 * aboveNow[0];
1088          t11 -= a01 * aboveNow[1];
1089          t21 -= a01 * aboveNow[2];
1090          t31 -= a01 * aboveNow[3];
1091          t02 -= a02 * aboveNow[0];
1092          t12 -= a02 * aboveNow[1];
1093          t22 -= a02 * aboveNow[2];
1094          t32 -= a02 * aboveNow[3];
1095          t03 -= a03 * aboveNow[0];
1096          t13 -= a03 * aboveNow[1];
1097          t23 -= a03 * aboveNow[2];
1098          t33 -= a03 * aboveNow[3];
1099          aUnderNow += BLOCK;
1100          aboveNow += BLOCK;
1101        }
1102        aa[i+0+0*BLOCK]=t00;
1103        aa[i+0+1*BLOCK]=t10;
1104        aa[i+0+2*BLOCK]=t20;
1105        aa[i+0+3*BLOCK]=t30;
1106        aa[i+1+0*BLOCK]=t01;
1107        aa[i+1+1*BLOCK]=t11;
1108        aa[i+1+2*BLOCK]=t21;
1109        aa[i+1+3*BLOCK]=t31;
1110        aa[i+2+0*BLOCK]=t02;
1111        aa[i+2+1*BLOCK]=t12;
1112        aa[i+2+2*BLOCK]=t22;
1113        aa[i+2+3*BLOCK]=t32;
1114        aa[i+3+0*BLOCK]=t03;
1115        aa[i+3+1*BLOCK]=t13;
1116        aa[i+3+2*BLOCK]=t23;
1117        aa[i+3+3*BLOCK]=t33;
1118      }
1119    }
1120#endif
1121  } else {
1122    int odd=nUnder&1;
1123    int n=nUnder-odd;
1124    for (j = 0; j < BLOCK; j +=4) {
1125      aa+=4*BLOCK;
1126      for (i=0;i< n;i+=2) {
1127        longDouble t00=aa[i+0*BLOCK];
1128        longDouble t10=aa[i+1*BLOCK];
1129        longDouble t20=aa[i+2*BLOCK];
1130        longDouble t30=aa[i+3*BLOCK];
1131        longDouble t01=aa[i+1+0*BLOCK];
1132        longDouble t11=aa[i+1+1*BLOCK];
1133        longDouble t21=aa[i+1+2*BLOCK];
1134        longDouble t31=aa[i+1+3*BLOCK];
1135        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
1136        const longDouble * COIN_RESTRICT aboveNow = above+j;
1137        for (k=0;k<BLOCK;k++) {
1138          longDouble multiplier = work[k];
1139          longDouble a00=aUnderNow[0]*multiplier;
1140          longDouble a01=aUnderNow[1]*multiplier;
1141          t00 -= a00 * aboveNow[0];
1142          t10 -= a00 * aboveNow[1];
1143          t20 -= a00 * aboveNow[2];
1144          t30 -= a00 * aboveNow[3];
1145          t01 -= a01 * aboveNow[0];
1146          t11 -= a01 * aboveNow[1];
1147          t21 -= a01 * aboveNow[2];
1148          t31 -= a01 * aboveNow[3];
1149          aUnderNow += BLOCK;
1150          aboveNow += BLOCK;
1151        }
1152        aa[i+0*BLOCK]=t00;
1153        aa[i+1*BLOCK]=t10;
1154        aa[i+2*BLOCK]=t20;
1155        aa[i+3*BLOCK]=t30;
1156        aa[i+1+0*BLOCK]=t01;
1157        aa[i+1+1*BLOCK]=t11;
1158        aa[i+1+2*BLOCK]=t21;
1159        aa[i+1+3*BLOCK]=t31;
1160      }
1161      if (odd) {
1162        longDouble t0=aa[n+0*BLOCK];
1163        longDouble t1=aa[n+1*BLOCK];
1164        longDouble t2=aa[n+2*BLOCK];
1165        longDouble t3=aa[n+3*BLOCK];
1166        longDouble a0;
1167        for (k=0;k<BLOCK;k++) {
1168          a0=aUnder[n+k*BLOCK]*work[k];
1169          t0 -= a0 * above[j + 0 + k * BLOCK];
1170          t1 -= a0 * above[j + 1 + k * BLOCK];
1171          t2 -= a0 * above[j + 2 + k * BLOCK];
1172          t3 -= a0 * above[j + 3 + k * BLOCK];
1173        }
1174        aa[n+0*BLOCK]=t0;
1175        aa[n+1*BLOCK]=t1;
1176        aa[n+2*BLOCK]=t2;
1177        aa[n+3*BLOCK]=t3;
1178      }
1179    }
1180  }
1181#else
1182  aa = aOther-BLOCK;
1183  for (j = 0; j < BLOCK; j ++) {
1184    aa+=BLOCK;
1185    for (i=0;i< nUnder;i++) {
1186      longDouble t00=aa[i+0*BLOCK];
1187      for (k=0;k<BLOCK;k++) {
1188        longDouble a00=aUnder[i+k*BLOCK]*work[k];
1189        t00 -= a00 * above[j + k * BLOCK];
1190      }
1191      aa[i]=t00;
1192    }
1193  }
1194#endif
1195}
1196/* Uses factorization to solve. */
1197void 
1198ClpCholeskyDense::solve (double * region) 
1199{
1200#ifdef CHOL_COMPARE
1201  double * region2=NULL;
1202  if (numberRows_<200) {
1203    longDouble * xx = sparseFactor_;
1204    longDouble * yy = diagonal_;
1205    diagonal_ = sparseFactor_ + 40000;
1206    sparseFactor_=diagonal_ + numberRows_;
1207    region2 = ClpCopyOfArray(region,numberRows_);
1208    int iRow,iColumn;
1209    int addOffset=numberRows_-1;
1210    longDouble * work = sparseFactor_-1;
1211    for (iColumn=0;iColumn<numberRows_;iColumn++) {
1212      double value = region2[iColumn];
1213      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
1214        region2[iRow] -= value*work[iRow];
1215      addOffset--;
1216      work += addOffset;
1217    }
1218    for (iColumn=numberRows_-1;iColumn>=0;iColumn--) {
1219      double value = region2[iColumn]*diagonal_[iColumn];
1220      work -= addOffset;
1221      addOffset++;
1222      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
1223        value -= region2[iRow]*work[iRow];
1224      region2[iColumn]=value;
1225    }
1226    sparseFactor_=xx;
1227    diagonal_=yy;
1228  }
1229#endif
1230  //longDouble * xx = sparseFactor_;
1231  //longDouble * yy = diagonal_;
1232  //diagonal_ = sparseFactor_ + 40000;
1233  //sparseFactor_=diagonal_ + numberRows_;
1234  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1235  // later align on boundary
1236  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1237  int iBlock;
1238  longDouble * aa = a;
1239  int iColumn;
1240  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1241    int nChunk;
1242    int jBlock;
1243    int iDo = iBlock*BLOCK;
1244    int base=iDo;
1245    if (iDo+BLOCK>numberRows_) {
1246      nChunk=numberRows_-iDo;
1247    } else {
1248      nChunk=BLOCK;
1249    }
1250    solveF1(aa,nChunk,region+iDo);
1251    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1252      base+=BLOCK;
1253      aa+=BLOCKSQ;
1254      if (base+BLOCK>numberRows_) {
1255        nChunk=numberRows_-base;
1256      } else {
1257        nChunk=BLOCK;
1258      }
1259      solveF2(aa,nChunk,region+iDo,region+base);
1260    }
1261    aa+=BLOCKSQ;
1262  }
1263  // do diagonal outside
1264  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1265    region[iColumn] *= diagonal_[iColumn];
1266  int offset=((numberBlocks*(numberBlocks+1))>>1);
1267  aa = a+number_entries(offset-1);
1268  int lBase=(numberBlocks-1)*BLOCK;
1269  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1270    int nChunk;
1271    int jBlock;
1272    int triBase=iBlock*BLOCK;
1273    int iBase=lBase;
1274    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1275      if (iBase+BLOCK>numberRows_) {
1276        nChunk=numberRows_-iBase;
1277      } else {
1278        nChunk=BLOCK;
1279      }
1280      solveB2(aa,nChunk,region+triBase,region+iBase);
1281      iBase-=BLOCK;
1282      aa-=BLOCKSQ;
1283    }
1284    if (triBase+BLOCK>numberRows_) {
1285      nChunk=numberRows_-triBase;
1286    } else {
1287      nChunk=BLOCK;
1288    }
1289    solveB1(aa,nChunk,region+triBase);
1290    aa-=BLOCKSQ;
1291  }
1292#ifdef CHOL_COMPARE
1293  if (numberRows_<200) {
1294    for (int i=0;i<numberRows_;i++) {
1295      assert(fabs(region[i]-region2[i])<1.0e-3);
1296    }
1297    delete [] region2;
1298  }
1299#endif
1300}
1301// Forward part of solve 1
1302void 
1303ClpCholeskyDense::solveF1(longDouble * a,int n,double * region)
1304{
1305  int j, k;
1306  longDouble t00;
1307  for (j = 0; j < n; j ++) {
1308    t00 = region[j];
1309    for (k = 0; k < j; ++k) {
1310      t00 -= region[k] * a[j + k * BLOCK];
1311    }
1312    //t00*=a[j + j * BLOCK];
1313    region[j] = t00;
1314  }
1315}
1316// Forward part of solve 2
1317void 
1318ClpCholeskyDense::solveF2(longDouble * a,int n,double * region, double * region2)
1319{
1320  int j, k;
1321#ifdef BLOCKUNROLL
1322  if (n==BLOCK) {
1323    for (k = 0; k < BLOCK; k+=4) {
1324      longDouble t0 = region2[0];
1325      longDouble t1 = region2[1];
1326      longDouble t2 = region2[2];
1327      longDouble t3 = region2[3];
1328      t0 -= region[0] * a[0 + 0 * BLOCK];
1329      t1 -= region[0] * a[1 + 0 * BLOCK];
1330      t2 -= region[0] * a[2 + 0 * BLOCK];
1331      t3 -= region[0] * a[3 + 0 * BLOCK];
1332
1333      t0 -= region[1] * a[0 + 1 * BLOCK];
1334      t1 -= region[1] * a[1 + 1 * BLOCK];
1335      t2 -= region[1] * a[2 + 1 * BLOCK];
1336      t3 -= region[1] * a[3 + 1 * BLOCK];
1337
1338      t0 -= region[2] * a[0 + 2 * BLOCK];
1339      t1 -= region[2] * a[1 + 2 * BLOCK];
1340      t2 -= region[2] * a[2 + 2 * BLOCK];
1341      t3 -= region[2] * a[3 + 2 * BLOCK];
1342
1343      t0 -= region[3] * a[0 + 3 * BLOCK];
1344      t1 -= region[3] * a[1 + 3 * BLOCK];
1345      t2 -= region[3] * a[2 + 3 * BLOCK];
1346      t3 -= region[3] * a[3 + 3 * BLOCK];
1347
1348      t0 -= region[4] * a[0 + 4 * BLOCK];
1349      t1 -= region[4] * a[1 + 4 * BLOCK];
1350      t2 -= region[4] * a[2 + 4 * BLOCK];
1351      t3 -= region[4] * a[3 + 4 * BLOCK];
1352
1353      t0 -= region[5] * a[0 + 5 * BLOCK];
1354      t1 -= region[5] * a[1 + 5 * BLOCK];
1355      t2 -= region[5] * a[2 + 5 * BLOCK];
1356      t3 -= region[5] * a[3 + 5 * BLOCK];
1357
1358      t0 -= region[6] * a[0 + 6 * BLOCK];
1359      t1 -= region[6] * a[1 + 6 * BLOCK];
1360      t2 -= region[6] * a[2 + 6 * BLOCK];
1361      t3 -= region[6] * a[3 + 6 * BLOCK];
1362
1363      t0 -= region[7] * a[0 + 7 * BLOCK];
1364      t1 -= region[7] * a[1 + 7 * BLOCK];
1365      t2 -= region[7] * a[2 + 7 * BLOCK];
1366      t3 -= region[7] * a[3 + 7 * BLOCK];
1367#if BLOCK>8
1368      t0 -= region[8] * a[0 + 8 * BLOCK];
1369      t1 -= region[8] * a[1 + 8 * BLOCK];
1370      t2 -= region[8] * a[2 + 8 * BLOCK];
1371      t3 -= region[8] * a[3 + 8 * BLOCK];
1372
1373      t0 -= region[9] * a[0 + 9 * BLOCK];
1374      t1 -= region[9] * a[1 + 9 * BLOCK];
1375      t2 -= region[9] * a[2 + 9 * BLOCK];
1376      t3 -= region[9] * a[3 + 9 * BLOCK];
1377
1378      t0 -= region[10] * a[0 + 10 * BLOCK];
1379      t1 -= region[10] * a[1 + 10 * BLOCK];
1380      t2 -= region[10] * a[2 + 10 * BLOCK];
1381      t3 -= region[10] * a[3 + 10 * BLOCK];
1382
1383      t0 -= region[11] * a[0 + 11 * BLOCK];
1384      t1 -= region[11] * a[1 + 11 * BLOCK];
1385      t2 -= region[11] * a[2 + 11 * BLOCK];
1386      t3 -= region[11] * a[3 + 11 * BLOCK];
1387
1388      t0 -= region[12] * a[0 + 12 * BLOCK];
1389      t1 -= region[12] * a[1 + 12 * BLOCK];
1390      t2 -= region[12] * a[2 + 12 * BLOCK];
1391      t3 -= region[12] * a[3 + 12 * BLOCK];
1392
1393      t0 -= region[13] * a[0 + 13 * BLOCK];
1394      t1 -= region[13] * a[1 + 13 * BLOCK];
1395      t2 -= region[13] * a[2 + 13 * BLOCK];
1396      t3 -= region[13] * a[3 + 13 * BLOCK];
1397
1398      t0 -= region[14] * a[0 + 14 * BLOCK];
1399      t1 -= region[14] * a[1 + 14 * BLOCK];
1400      t2 -= region[14] * a[2 + 14 * BLOCK];
1401      t3 -= region[14] * a[3 + 14 * BLOCK];
1402
1403      t0 -= region[15] * a[0 + 15 * BLOCK];
1404      t1 -= region[15] * a[1 + 15 * BLOCK];
1405      t2 -= region[15] * a[2 + 15 * BLOCK];
1406      t3 -= region[15] * a[3 + 15 * BLOCK];
1407#endif
1408      region2[0] = t0;
1409      region2[1] = t1;
1410      region2[2] = t2;
1411      region2[3] = t3;
1412      region2+=4;
1413      a+=4;
1414    }
1415  } else {
1416#endif
1417    for (k = 0; k < n; ++k) {
1418      longDouble t00 = region2[k];
1419      for (j = 0; j < BLOCK; j ++) {
1420        t00 -= region[j] * a[k + j * BLOCK];
1421      }
1422      region2[k] = t00;
1423    }
1424#ifdef BLOCKUNROLL
1425  }
1426#endif
1427}
1428// Backward part of solve 1
1429void 
1430ClpCholeskyDense::solveB1(longDouble * a,int n,double * region)
1431{
1432  int j, k;
1433  longDouble t00;
1434  for (j = n-1; j >=0; j --) {
1435    t00 = region[j];
1436    for (k = j+1; k < n; ++k) {
1437      t00 -= region[k] * a[k + j * BLOCK];
1438    }
1439    //t00*=a[j + j * BLOCK];
1440    region[j] = t00;
1441  }
1442}
1443// Backward part of solve 2
1444void 
1445ClpCholeskyDense::solveB2(longDouble * a,int n,double * region, double * region2)
1446{
1447  int j, k;
1448#ifdef BLOCKUNROLL
1449  if (n==BLOCK) {
1450    for (j = 0; j < BLOCK; j +=4) {
1451      longDouble t0 = region[0];
1452      longDouble t1 = region[1];
1453      longDouble t2 = region[2];
1454      longDouble t3 = region[3];
1455      t0 -= region2[0] * a[0 + 0*BLOCK];
1456      t1 -= region2[0] * a[0 + 1*BLOCK];
1457      t2 -= region2[0] * a[0 + 2*BLOCK];
1458      t3 -= region2[0] * a[0 + 3*BLOCK];
1459
1460      t0 -= region2[1] * a[1 + 0*BLOCK];
1461      t1 -= region2[1] * a[1 + 1*BLOCK];
1462      t2 -= region2[1] * a[1 + 2*BLOCK];
1463      t3 -= region2[1] * a[1 + 3*BLOCK];
1464
1465      t0 -= region2[2] * a[2 + 0*BLOCK];
1466      t1 -= region2[2] * a[2 + 1*BLOCK];
1467      t2 -= region2[2] * a[2 + 2*BLOCK];
1468      t3 -= region2[2] * a[2 + 3*BLOCK];
1469
1470      t0 -= region2[3] * a[3 + 0*BLOCK];
1471      t1 -= region2[3] * a[3 + 1*BLOCK];
1472      t2 -= region2[3] * a[3 + 2*BLOCK];
1473      t3 -= region2[3] * a[3 + 3*BLOCK];
1474
1475      t0 -= region2[4] * a[4 + 0*BLOCK];
1476      t1 -= region2[4] * a[4 + 1*BLOCK];
1477      t2 -= region2[4] * a[4 + 2*BLOCK];
1478      t3 -= region2[4] * a[4 + 3*BLOCK];
1479
1480      t0 -= region2[5] * a[5 + 0*BLOCK];
1481      t1 -= region2[5] * a[5 + 1*BLOCK];
1482      t2 -= region2[5] * a[5 + 2*BLOCK];
1483      t3 -= region2[5] * a[5 + 3*BLOCK];
1484
1485      t0 -= region2[6] * a[6 + 0*BLOCK];
1486      t1 -= region2[6] * a[6 + 1*BLOCK];
1487      t2 -= region2[6] * a[6 + 2*BLOCK];
1488      t3 -= region2[6] * a[6 + 3*BLOCK];
1489
1490      t0 -= region2[7] * a[7 + 0*BLOCK];
1491      t1 -= region2[7] * a[7 + 1*BLOCK];
1492      t2 -= region2[7] * a[7 + 2*BLOCK];
1493      t3 -= region2[7] * a[7 + 3*BLOCK];
1494#if BLOCK>8
1495
1496      t0 -= region2[8] * a[8 + 0*BLOCK];
1497      t1 -= region2[8] * a[8 + 1*BLOCK];
1498      t2 -= region2[8] * a[8 + 2*BLOCK];
1499      t3 -= region2[8] * a[8 + 3*BLOCK];
1500
1501      t0 -= region2[9] * a[9 + 0*BLOCK];
1502      t1 -= region2[9] * a[9 + 1*BLOCK];
1503      t2 -= region2[9] * a[9 + 2*BLOCK];
1504      t3 -= region2[9] * a[9 + 3*BLOCK];
1505
1506      t0 -= region2[10] * a[10 + 0*BLOCK];
1507      t1 -= region2[10] * a[10 + 1*BLOCK];
1508      t2 -= region2[10] * a[10 + 2*BLOCK];
1509      t3 -= region2[10] * a[10 + 3*BLOCK];
1510
1511      t0 -= region2[11] * a[11 + 0*BLOCK];
1512      t1 -= region2[11] * a[11 + 1*BLOCK];
1513      t2 -= region2[11] * a[11 + 2*BLOCK];
1514      t3 -= region2[11] * a[11 + 3*BLOCK];
1515
1516      t0 -= region2[12] * a[12 + 0*BLOCK];
1517      t1 -= region2[12] * a[12 + 1*BLOCK];
1518      t2 -= region2[12] * a[12 + 2*BLOCK];
1519      t3 -= region2[12] * a[12 + 3*BLOCK];
1520
1521      t0 -= region2[13] * a[13 + 0*BLOCK];
1522      t1 -= region2[13] * a[13 + 1*BLOCK];
1523      t2 -= region2[13] * a[13 + 2*BLOCK];
1524      t3 -= region2[13] * a[13 + 3*BLOCK];
1525
1526      t0 -= region2[14] * a[14 + 0*BLOCK];
1527      t1 -= region2[14] * a[14 + 1*BLOCK];
1528      t2 -= region2[14] * a[14 + 2*BLOCK];
1529      t3 -= region2[14] * a[14 + 3*BLOCK];
1530
1531      t0 -= region2[15] * a[15 + 0*BLOCK];
1532      t1 -= region2[15] * a[15 + 1*BLOCK];
1533      t2 -= region2[15] * a[15 + 2*BLOCK];
1534      t3 -= region2[15] * a[15 + 3*BLOCK];
1535#endif
1536      region[0] = t0;
1537      region[1] = t1;
1538      region[2] = t2;
1539      region[3] = t3;
1540      a+=4*BLOCK;
1541      region+=4;
1542    }
1543  } else {
1544#endif
1545    for (j = 0; j < BLOCK; j ++) {
1546      longDouble t00 = region[j];
1547      for (k = 0; k < n; ++k) {
1548        t00 -= region2[k] * a[k + j * BLOCK];
1549      }
1550      region[j] = t00;
1551    }
1552#ifdef BLOCKUNROLL
1553  }
1554#endif
1555}
1556/* Uses factorization to solve. */
1557void 
1558ClpCholeskyDense::solveLong (longDouble * region) 
1559{
1560  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1561  // later align on boundary
1562  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1563  int iBlock;
1564  longDouble * aa = a;
1565  int iColumn;
1566  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1567    int nChunk;
1568    int jBlock;
1569    int iDo = iBlock*BLOCK;
1570    int base=iDo;
1571    if (iDo+BLOCK>numberRows_) {
1572      nChunk=numberRows_-iDo;
1573    } else {
1574      nChunk=BLOCK;
1575    }
1576    solveF1Long(aa,nChunk,region+iDo);
1577    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1578      base+=BLOCK;
1579      aa+=BLOCKSQ;
1580      if (base+BLOCK>numberRows_) {
1581        nChunk=numberRows_-base;
1582      } else {
1583        nChunk=BLOCK;
1584      }
1585      solveF2Long(aa,nChunk,region+iDo,region+base);
1586    }
1587    aa+=BLOCKSQ;
1588  }
1589  // do diagonal outside
1590  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1591    region[iColumn] *= diagonal_[iColumn];
1592  int offset=((numberBlocks*(numberBlocks+1))>>1);
1593  aa = a+number_entries(offset-1);
1594  int lBase=(numberBlocks-1)*BLOCK;
1595  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1596    int nChunk;
1597    int jBlock;
1598    int triBase=iBlock*BLOCK;
1599    int iBase=lBase;
1600    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1601      if (iBase+BLOCK>numberRows_) {
1602        nChunk=numberRows_-iBase;
1603      } else {
1604        nChunk=BLOCK;
1605      }
1606      solveB2Long(aa,nChunk,region+triBase,region+iBase);
1607      iBase-=BLOCK;
1608      aa-=BLOCKSQ;
1609    }
1610    if (triBase+BLOCK>numberRows_) {
1611      nChunk=numberRows_-triBase;
1612    } else {
1613      nChunk=BLOCK;
1614    }
1615    solveB1Long(aa,nChunk,region+triBase);
1616    aa-=BLOCKSQ;
1617  }
1618}
1619// Forward part of solve 1
1620void 
1621ClpCholeskyDense::solveF1Long(longDouble * a,int n,longDouble * region)
1622{
1623  int j, k;
1624  longDouble t00;
1625  for (j = 0; j < n; j ++) {
1626    t00 = region[j];
1627    for (k = 0; k < j; ++k) {
1628      t00 -= region[k] * a[j + k * BLOCK];
1629    }
1630    //t00*=a[j + j * BLOCK];
1631    region[j] = t00;
1632  }
1633}
1634// Forward part of solve 2
1635void 
1636ClpCholeskyDense::solveF2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
1637{
1638  int j, k;
1639#ifdef BLOCKUNROLL
1640  if (n==BLOCK) {
1641    for (k = 0; k < BLOCK; k+=4) {
1642      longDouble t0 = region2[0];
1643      longDouble t1 = region2[1];
1644      longDouble t2 = region2[2];
1645      longDouble t3 = region2[3];
1646      t0 -= region[0] * a[0 + 0 * BLOCK];
1647      t1 -= region[0] * a[1 + 0 * BLOCK];
1648      t2 -= region[0] * a[2 + 0 * BLOCK];
1649      t3 -= region[0] * a[3 + 0 * BLOCK];
1650
1651      t0 -= region[1] * a[0 + 1 * BLOCK];
1652      t1 -= region[1] * a[1 + 1 * BLOCK];
1653      t2 -= region[1] * a[2 + 1 * BLOCK];
1654      t3 -= region[1] * a[3 + 1 * BLOCK];
1655
1656      t0 -= region[2] * a[0 + 2 * BLOCK];
1657      t1 -= region[2] * a[1 + 2 * BLOCK];
1658      t2 -= region[2] * a[2 + 2 * BLOCK];
1659      t3 -= region[2] * a[3 + 2 * BLOCK];
1660
1661      t0 -= region[3] * a[0 + 3 * BLOCK];
1662      t1 -= region[3] * a[1 + 3 * BLOCK];
1663      t2 -= region[3] * a[2 + 3 * BLOCK];
1664      t3 -= region[3] * a[3 + 3 * BLOCK];
1665
1666      t0 -= region[4] * a[0 + 4 * BLOCK];
1667      t1 -= region[4] * a[1 + 4 * BLOCK];
1668      t2 -= region[4] * a[2 + 4 * BLOCK];
1669      t3 -= region[4] * a[3 + 4 * BLOCK];
1670
1671      t0 -= region[5] * a[0 + 5 * BLOCK];
1672      t1 -= region[5] * a[1 + 5 * BLOCK];
1673      t2 -= region[5] * a[2 + 5 * BLOCK];
1674      t3 -= region[5] * a[3 + 5 * BLOCK];
1675
1676      t0 -= region[6] * a[0 + 6 * BLOCK];
1677      t1 -= region[6] * a[1 + 6 * BLOCK];
1678      t2 -= region[6] * a[2 + 6 * BLOCK];
1679      t3 -= region[6] * a[3 + 6 * BLOCK];
1680
1681      t0 -= region[7] * a[0 + 7 * BLOCK];
1682      t1 -= region[7] * a[1 + 7 * BLOCK];
1683      t2 -= region[7] * a[2 + 7 * BLOCK];
1684      t3 -= region[7] * a[3 + 7 * BLOCK];
1685#if BLOCK>8
1686      t0 -= region[8] * a[0 + 8 * BLOCK];
1687      t1 -= region[8] * a[1 + 8 * BLOCK];
1688      t2 -= region[8] * a[2 + 8 * BLOCK];
1689      t3 -= region[8] * a[3 + 8 * BLOCK];
1690
1691      t0 -= region[9] * a[0 + 9 * BLOCK];
1692      t1 -= region[9] * a[1 + 9 * BLOCK];
1693      t2 -= region[9] * a[2 + 9 * BLOCK];
1694      t3 -= region[9] * a[3 + 9 * BLOCK];
1695
1696      t0 -= region[10] * a[0 + 10 * BLOCK];
1697      t1 -= region[10] * a[1 + 10 * BLOCK];
1698      t2 -= region[10] * a[2 + 10 * BLOCK];
1699      t3 -= region[10] * a[3 + 10 * BLOCK];
1700
1701      t0 -= region[11] * a[0 + 11 * BLOCK];
1702      t1 -= region[11] * a[1 + 11 * BLOCK];
1703      t2 -= region[11] * a[2 + 11 * BLOCK];
1704      t3 -= region[11] * a[3 + 11 * BLOCK];
1705
1706      t0 -= region[12] * a[0 + 12 * BLOCK];
1707      t1 -= region[12] * a[1 + 12 * BLOCK];
1708      t2 -= region[12] * a[2 + 12 * BLOCK];
1709      t3 -= region[12] * a[3 + 12 * BLOCK];
1710
1711      t0 -= region[13] * a[0 + 13 * BLOCK];
1712      t1 -= region[13] * a[1 + 13 * BLOCK];
1713      t2 -= region[13] * a[2 + 13 * BLOCK];
1714      t3 -= region[13] * a[3 + 13 * BLOCK];
1715
1716      t0 -= region[14] * a[0 + 14 * BLOCK];
1717      t1 -= region[14] * a[1 + 14 * BLOCK];
1718      t2 -= region[14] * a[2 + 14 * BLOCK];
1719      t3 -= region[14] * a[3 + 14 * BLOCK];
1720
1721      t0 -= region[15] * a[0 + 15 * BLOCK];
1722      t1 -= region[15] * a[1 + 15 * BLOCK];
1723      t2 -= region[15] * a[2 + 15 * BLOCK];
1724      t3 -= region[15] * a[3 + 15 * BLOCK];
1725#endif
1726      region2[0] = t0;
1727      region2[1] = t1;
1728      region2[2] = t2;
1729      region2[3] = t3;
1730      region2+=4;
1731      a+=4;
1732    }
1733  } else {
1734#endif
1735    for (k = 0; k < n; ++k) {
1736      longDouble t00 = region2[k];
1737      for (j = 0; j < BLOCK; j ++) {
1738        t00 -= region[j] * a[k + j * BLOCK];
1739      }
1740      region2[k] = t00;
1741    }
1742#ifdef BLOCKUNROLL
1743  }
1744#endif
1745}
1746// Backward part of solve 1
1747void 
1748ClpCholeskyDense::solveB1Long(longDouble * a,int n,longDouble * region)
1749{
1750  int j, k;
1751  longDouble t00;
1752  for (j = n-1; j >=0; j --) {
1753    t00 = region[j];
1754    for (k = j+1; k < n; ++k) {
1755      t00 -= region[k] * a[k + j * BLOCK];
1756    }
1757    //t00*=a[j + j * BLOCK];
1758    region[j] = t00;
1759  }
1760}
1761// Backward part of solve 2
1762void 
1763ClpCholeskyDense::solveB2Long(longDouble * a,int n,longDouble * region, longDouble * region2)
1764{
1765  int j, k;
1766#ifdef BLOCKUNROLL
1767  if (n==BLOCK) {
1768    for (j = 0; j < BLOCK; j +=4) {
1769      longDouble t0 = region[0];
1770      longDouble t1 = region[1];
1771      longDouble t2 = region[2];
1772      longDouble t3 = region[3];
1773      t0 -= region2[0] * a[0 + 0*BLOCK];
1774      t1 -= region2[0] * a[0 + 1*BLOCK];
1775      t2 -= region2[0] * a[0 + 2*BLOCK];
1776      t3 -= region2[0] * a[0 + 3*BLOCK];
1777
1778      t0 -= region2[1] * a[1 + 0*BLOCK];
1779      t1 -= region2[1] * a[1 + 1*BLOCK];
1780      t2 -= region2[1] * a[1 + 2*BLOCK];
1781      t3 -= region2[1] * a[1 + 3*BLOCK];
1782
1783      t0 -= region2[2] * a[2 + 0*BLOCK];
1784      t1 -= region2[2] * a[2 + 1*BLOCK];
1785      t2 -= region2[2] * a[2 + 2*BLOCK];
1786      t3 -= region2[2] * a[2 + 3*BLOCK];
1787
1788      t0 -= region2[3] * a[3 + 0*BLOCK];
1789      t1 -= region2[3] * a[3 + 1*BLOCK];
1790      t2 -= region2[3] * a[3 + 2*BLOCK];
1791      t3 -= region2[3] * a[3 + 3*BLOCK];
1792
1793      t0 -= region2[4] * a[4 + 0*BLOCK];
1794      t1 -= region2[4] * a[4 + 1*BLOCK];
1795      t2 -= region2[4] * a[4 + 2*BLOCK];
1796      t3 -= region2[4] * a[4 + 3*BLOCK];
1797
1798      t0 -= region2[5] * a[5 + 0*BLOCK];
1799      t1 -= region2[5] * a[5 + 1*BLOCK];
1800      t2 -= region2[5] * a[5 + 2*BLOCK];
1801      t3 -= region2[5] * a[5 + 3*BLOCK];
1802
1803      t0 -= region2[6] * a[6 + 0*BLOCK];
1804      t1 -= region2[6] * a[6 + 1*BLOCK];
1805      t2 -= region2[6] * a[6 + 2*BLOCK];
1806      t3 -= region2[6] * a[6 + 3*BLOCK];
1807
1808      t0 -= region2[7] * a[7 + 0*BLOCK];
1809      t1 -= region2[7] * a[7 + 1*BLOCK];
1810      t2 -= region2[7] * a[7 + 2*BLOCK];
1811      t3 -= region2[7] * a[7 + 3*BLOCK];
1812#if BLOCK>8
1813
1814      t0 -= region2[8] * a[8 + 0*BLOCK];
1815      t1 -= region2[8] * a[8 + 1*BLOCK];
1816      t2 -= region2[8] * a[8 + 2*BLOCK];
1817      t3 -= region2[8] * a[8 + 3*BLOCK];
1818
1819      t0 -= region2[9] * a[9 + 0*BLOCK];
1820      t1 -= region2[9] * a[9 + 1*BLOCK];
1821      t2 -= region2[9] * a[9 + 2*BLOCK];
1822      t3 -= region2[9] * a[9 + 3*BLOCK];
1823
1824      t0 -= region2[10] * a[10 + 0*BLOCK];
1825      t1 -= region2[10] * a[10 + 1*BLOCK];
1826      t2 -= region2[10] * a[10 + 2*BLOCK];
1827      t3 -= region2[10] * a[10 + 3*BLOCK];
1828
1829      t0 -= region2[11] * a[11 + 0*BLOCK];
1830      t1 -= region2[11] * a[11 + 1*BLOCK];
1831      t2 -= region2[11] * a[11 + 2*BLOCK];
1832      t3 -= region2[11] * a[11 + 3*BLOCK];
1833
1834      t0 -= region2[12] * a[12 + 0*BLOCK];
1835      t1 -= region2[12] * a[12 + 1*BLOCK];
1836      t2 -= region2[12] * a[12 + 2*BLOCK];
1837      t3 -= region2[12] * a[12 + 3*BLOCK];
1838
1839      t0 -= region2[13] * a[13 + 0*BLOCK];
1840      t1 -= region2[13] * a[13 + 1*BLOCK];
1841      t2 -= region2[13] * a[13 + 2*BLOCK];
1842      t3 -= region2[13] * a[13 + 3*BLOCK];
1843
1844      t0 -= region2[14] * a[14 + 0*BLOCK];
1845      t1 -= region2[14] * a[14 + 1*BLOCK];
1846      t2 -= region2[14] * a[14 + 2*BLOCK];
1847      t3 -= region2[14] * a[14 + 3*BLOCK];
1848
1849      t0 -= region2[15] * a[15 + 0*BLOCK];
1850      t1 -= region2[15] * a[15 + 1*BLOCK];
1851      t2 -= region2[15] * a[15 + 2*BLOCK];
1852      t3 -= region2[15] * a[15 + 3*BLOCK];
1853#endif
1854      region[0] = t0;
1855      region[1] = t1;
1856      region[2] = t2;
1857      region[3] = t3;
1858      a+=4*BLOCK;
1859      region+=4;
1860    }
1861  } else {
1862#endif
1863    for (j = 0; j < BLOCK; j ++) {
1864      longDouble t00 = region[j];
1865      for (k = 0; k < n; ++k) {
1866        t00 -= region2[k] * a[k + j * BLOCK];
1867      }
1868      region[j] = t00;
1869    }
1870#ifdef BLOCKUNROLL
1871  }
1872#endif
1873}
1874/* Uses factorization to solve. */
1875void 
1876ClpCholeskyDense::solveLongWork (longWork * region) 
1877{
1878  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1879  // later align on boundary
1880  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1881  int iBlock;
1882  longDouble * aa = a;
1883  int iColumn;
1884  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1885    int nChunk;
1886    int jBlock;
1887    int iDo = iBlock*BLOCK;
1888    int base=iDo;
1889    if (iDo+BLOCK>numberRows_) {
1890      nChunk=numberRows_-iDo;
1891    } else {
1892      nChunk=BLOCK;
1893    }
1894    solveF1LongWork(aa,nChunk,region+iDo);
1895    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1896      base+=BLOCK;
1897      aa+=BLOCKSQ;
1898      if (base+BLOCK>numberRows_) {
1899        nChunk=numberRows_-base;
1900      } else {
1901        nChunk=BLOCK;
1902      }
1903      solveF2LongWork(aa,nChunk,region+iDo,region+base);
1904    }
1905    aa+=BLOCKSQ;
1906  }
1907  // do diagonal outside
1908  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1909    region[iColumn] *= diagonal_[iColumn];
1910  int offset=((numberBlocks*(numberBlocks+1))>>1);
1911  aa = a+number_entries(offset-1);
1912  int lBase=(numberBlocks-1)*BLOCK;
1913  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1914    int nChunk;
1915    int jBlock;
1916    int triBase=iBlock*BLOCK;
1917    int iBase=lBase;
1918    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1919      if (iBase+BLOCK>numberRows_) {
1920        nChunk=numberRows_-iBase;
1921      } else {
1922        nChunk=BLOCK;
1923      }
1924      solveB2LongWork(aa,nChunk,region+triBase,region+iBase);
1925      iBase-=BLOCK;
1926      aa-=BLOCKSQ;
1927    }
1928    if (triBase+BLOCK>numberRows_) {
1929      nChunk=numberRows_-triBase;
1930    } else {
1931      nChunk=BLOCK;
1932    }
1933    solveB1LongWork(aa,nChunk,region+triBase);
1934    aa-=BLOCKSQ;
1935  }
1936}
1937// Forward part of solve 1
1938void 
1939ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,longWork * region)
1940{
1941  int j, k;
1942  longWork t00;
1943  for (j = 0; j < n; j ++) {
1944    t00 = region[j];
1945    for (k = 0; k < j; ++k) {
1946      t00 -= region[k] * a[j + k * BLOCK];
1947    }
1948    //t00*=a[j + j * BLOCK];
1949    region[j] = t00;
1950  }
1951}
1952// Forward part of solve 2
1953void 
1954ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
1955{
1956  int j, k;
1957#ifdef BLOCKUNROLL
1958  if (n==BLOCK) {
1959    for (k = 0; k < BLOCK; k+=4) {
1960      longWork t0 = region2[0];
1961      longWork t1 = region2[1];
1962      longWork t2 = region2[2];
1963      longWork t3 = region2[3];
1964      t0 -= region[0] * a[0 + 0 * BLOCK];
1965      t1 -= region[0] * a[1 + 0 * BLOCK];
1966      t2 -= region[0] * a[2 + 0 * BLOCK];
1967      t3 -= region[0] * a[3 + 0 * BLOCK];
1968
1969      t0 -= region[1] * a[0 + 1 * BLOCK];
1970      t1 -= region[1] * a[1 + 1 * BLOCK];
1971      t2 -= region[1] * a[2 + 1 * BLOCK];
1972      t3 -= region[1] * a[3 + 1 * BLOCK];
1973
1974      t0 -= region[2] * a[0 + 2 * BLOCK];
1975      t1 -= region[2] * a[1 + 2 * BLOCK];
1976      t2 -= region[2] * a[2 + 2 * BLOCK];
1977      t3 -= region[2] * a[3 + 2 * BLOCK];
1978
1979      t0 -= region[3] * a[0 + 3 * BLOCK];
1980      t1 -= region[3] * a[1 + 3 * BLOCK];
1981      t2 -= region[3] * a[2 + 3 * BLOCK];
1982      t3 -= region[3] * a[3 + 3 * BLOCK];
1983
1984      t0 -= region[4] * a[0 + 4 * BLOCK];
1985      t1 -= region[4] * a[1 + 4 * BLOCK];
1986      t2 -= region[4] * a[2 + 4 * BLOCK];
1987      t3 -= region[4] * a[3 + 4 * BLOCK];
1988
1989      t0 -= region[5] * a[0 + 5 * BLOCK];
1990      t1 -= region[5] * a[1 + 5 * BLOCK];
1991      t2 -= region[5] * a[2 + 5 * BLOCK];
1992      t3 -= region[5] * a[3 + 5 * BLOCK];
1993
1994      t0 -= region[6] * a[0 + 6 * BLOCK];
1995      t1 -= region[6] * a[1 + 6 * BLOCK];
1996      t2 -= region[6] * a[2 + 6 * BLOCK];
1997      t3 -= region[6] * a[3 + 6 * BLOCK];
1998
1999      t0 -= region[7] * a[0 + 7 * BLOCK];
2000      t1 -= region[7] * a[1 + 7 * BLOCK];
2001      t2 -= region[7] * a[2 + 7 * BLOCK];
2002      t3 -= region[7] * a[3 + 7 * BLOCK];
2003#if BLOCK>8
2004      t0 -= region[8] * a[0 + 8 * BLOCK];
2005      t1 -= region[8] * a[1 + 8 * BLOCK];
2006      t2 -= region[8] * a[2 + 8 * BLOCK];
2007      t3 -= region[8] * a[3 + 8 * BLOCK];
2008
2009      t0 -= region[9] * a[0 + 9 * BLOCK];
2010      t1 -= region[9] * a[1 + 9 * BLOCK];
2011      t2 -= region[9] * a[2 + 9 * BLOCK];
2012      t3 -= region[9] * a[3 + 9 * BLOCK];
2013
2014      t0 -= region[10] * a[0 + 10 * BLOCK];
2015      t1 -= region[10] * a[1 + 10 * BLOCK];
2016      t2 -= region[10] * a[2 + 10 * BLOCK];
2017      t3 -= region[10] * a[3 + 10 * BLOCK];
2018
2019      t0 -= region[11] * a[0 + 11 * BLOCK];
2020      t1 -= region[11] * a[1 + 11 * BLOCK];
2021      t2 -= region[11] * a[2 + 11 * BLOCK];
2022      t3 -= region[11] * a[3 + 11 * BLOCK];
2023
2024      t0 -= region[12] * a[0 + 12 * BLOCK];
2025      t1 -= region[12] * a[1 + 12 * BLOCK];
2026      t2 -= region[12] * a[2 + 12 * BLOCK];
2027      t3 -= region[12] * a[3 + 12 * BLOCK];
2028
2029      t0 -= region[13] * a[0 + 13 * BLOCK];
2030      t1 -= region[13] * a[1 + 13 * BLOCK];
2031      t2 -= region[13] * a[2 + 13 * BLOCK];
2032      t3 -= region[13] * a[3 + 13 * BLOCK];
2033
2034      t0 -= region[14] * a[0 + 14 * BLOCK];
2035      t1 -= region[14] * a[1 + 14 * BLOCK];
2036      t2 -= region[14] * a[2 + 14 * BLOCK];
2037      t3 -= region[14] * a[3 + 14 * BLOCK];
2038
2039      t0 -= region[15] * a[0 + 15 * BLOCK];
2040      t1 -= region[15] * a[1 + 15 * BLOCK];
2041      t2 -= region[15] * a[2 + 15 * BLOCK];
2042      t3 -= region[15] * a[3 + 15 * BLOCK];
2043#endif
2044      region2[0] = t0;
2045      region2[1] = t1;
2046      region2[2] = t2;
2047      region2[3] = t3;
2048      region2+=4;
2049      a+=4;
2050    }
2051  } else {
2052#endif
2053    for (k = 0; k < n; ++k) {
2054      longWork t00 = region2[k];
2055      for (j = 0; j < BLOCK; j ++) {
2056        t00 -= region[j] * a[k + j * BLOCK];
2057      }
2058      region2[k] = t00;
2059    }
2060#ifdef BLOCKUNROLL
2061  }
2062#endif
2063}
2064// Backward part of solve 1
2065void 
2066ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,longWork * region)
2067{
2068  int j, k;
2069  longWork t00;
2070  for (j = n-1; j >=0; j --) {
2071    t00 = region[j];
2072    for (k = j+1; k < n; ++k) {
2073      t00 -= region[k] * a[k + j * BLOCK];
2074    }
2075    //t00*=a[j + j * BLOCK];
2076    region[j] = t00;
2077  }
2078}
2079// Backward part of solve 2
2080void 
2081ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,longWork * region, longWork * region2)
2082{
2083  int j, k;
2084#ifdef BLOCKUNROLL
2085  if (n==BLOCK) {
2086    for (j = 0; j < BLOCK; j +=4) {
2087      longWork t0 = region[0];
2088      longWork t1 = region[1];
2089      longWork t2 = region[2];
2090      longWork t3 = region[3];
2091      t0 -= region2[0] * a[0 + 0*BLOCK];
2092      t1 -= region2[0] * a[0 + 1*BLOCK];
2093      t2 -= region2[0] * a[0 + 2*BLOCK];
2094      t3 -= region2[0] * a[0 + 3*BLOCK];
2095
2096      t0 -= region2[1] * a[1 + 0*BLOCK];
2097      t1 -= region2[1] * a[1 + 1*BLOCK];
2098      t2 -= region2[1] * a[1 + 2*BLOCK];
2099      t3 -= region2[1] * a[1 + 3*BLOCK];
2100
2101      t0 -= region2[2] * a[2 + 0*BLOCK];
2102      t1 -= region2[2] * a[2 + 1*BLOCK];
2103      t2 -= region2[2] * a[2 + 2*BLOCK];
2104      t3 -= region2[2] * a[2 + 3*BLOCK];
2105
2106      t0 -= region2[3] * a[3 + 0*BLOCK];
2107      t1 -= region2[3] * a[3 + 1*BLOCK];
2108      t2 -= region2[3] * a[3 + 2*BLOCK];
2109      t3 -= region2[3] * a[3 + 3*BLOCK];
2110
2111      t0 -= region2[4] * a[4 + 0*BLOCK];
2112      t1 -= region2[4] * a[4 + 1*BLOCK];
2113      t2 -= region2[4] * a[4 + 2*BLOCK];
2114      t3 -= region2[4] * a[4 + 3*BLOCK];
2115
2116      t0 -= region2[5] * a[5 + 0*BLOCK];
2117      t1 -= region2[5] * a[5 + 1*BLOCK];
2118      t2 -= region2[5] * a[5 + 2*BLOCK];
2119      t3 -= region2[5] * a[5 + 3*BLOCK];
2120
2121      t0 -= region2[6] * a[6 + 0*BLOCK];
2122      t1 -= region2[6] * a[6 + 1*BLOCK];
2123      t2 -= region2[6] * a[6 + 2*BLOCK];
2124      t3 -= region2[6] * a[6 + 3*BLOCK];
2125
2126      t0 -= region2[7] * a[7 + 0*BLOCK];
2127      t1 -= region2[7] * a[7 + 1*BLOCK];
2128      t2 -= region2[7] * a[7 + 2*BLOCK];
2129      t3 -= region2[7] * a[7 + 3*BLOCK];
2130#if BLOCK>8
2131
2132      t0 -= region2[8] * a[8 + 0*BLOCK];
2133      t1 -= region2[8] * a[8 + 1*BLOCK];
2134      t2 -= region2[8] * a[8 + 2*BLOCK];
2135      t3 -= region2[8] * a[8 + 3*BLOCK];
2136
2137      t0 -= region2[9] * a[9 + 0*BLOCK];
2138      t1 -= region2[9] * a[9 + 1*BLOCK];
2139      t2 -= region2[9] * a[9 + 2*BLOCK];
2140      t3 -= region2[9] * a[9 + 3*BLOCK];
2141
2142      t0 -= region2[10] * a[10 + 0*BLOCK];
2143      t1 -= region2[10] * a[10 + 1*BLOCK];
2144      t2 -= region2[10] * a[10 + 2*BLOCK];
2145      t3 -= region2[10] * a[10 + 3*BLOCK];
2146
2147      t0 -= region2[11] * a[11 + 0*BLOCK];
2148      t1 -= region2[11] * a[11 + 1*BLOCK];
2149      t2 -= region2[11] * a[11 + 2*BLOCK];
2150      t3 -= region2[11] * a[11 + 3*BLOCK];
2151
2152      t0 -= region2[12] * a[12 + 0*BLOCK];
2153      t1 -= region2[12] * a[12 + 1*BLOCK];
2154      t2 -= region2[12] * a[12 + 2*BLOCK];
2155      t3 -= region2[12] * a[12 + 3*BLOCK];
2156
2157      t0 -= region2[13] * a[13 + 0*BLOCK];
2158      t1 -= region2[13] * a[13 + 1*BLOCK];
2159      t2 -= region2[13] * a[13 + 2*BLOCK];
2160      t3 -= region2[13] * a[13 + 3*BLOCK];
2161
2162      t0 -= region2[14] * a[14 + 0*BLOCK];
2163      t1 -= region2[14] * a[14 + 1*BLOCK];
2164      t2 -= region2[14] * a[14 + 2*BLOCK];
2165      t3 -= region2[14] * a[14 + 3*BLOCK];
2166
2167      t0 -= region2[15] * a[15 + 0*BLOCK];
2168      t1 -= region2[15] * a[15 + 1*BLOCK];
2169      t2 -= region2[15] * a[15 + 2*BLOCK];
2170      t3 -= region2[15] * a[15 + 3*BLOCK];
2171#endif
2172      region[0] = t0;
2173      region[1] = t1;
2174      region[2] = t2;
2175      region[3] = t3;
2176      a+=4*BLOCK;
2177      region+=4;
2178    }
2179  } else {
2180#endif
2181    for (j = 0; j < BLOCK; j ++) {
2182      longWork t00 = region[j];
2183      for (k = 0; k < n; ++k) {
2184        t00 -= region2[k] * a[k + j * BLOCK];
2185      }
2186      region[j] = t00;
2187    }
2188#ifdef BLOCKUNROLL
2189  }
2190#endif
2191}
Note: See TracBrowser for help on using the repository browser.