source: releases/1.8.0/Clp/src/ClpCholeskyDense.cpp @ 2257

Last change on this file since 2257 was 1197, checked in by forrest, 12 years ago

many changes to try and improve performance

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