source: trunk/ClpCholeskyWssmp.cpp @ 320

Last change on this file since 320 was 320, checked in by forrest, 16 years ago

Cholesky stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.8 KB
Line 
1#ifdef REAL_BARRIER
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5
6
7#include "CoinPragma.hpp"
8#include "CoinHelperFunctions.hpp"
9#include "ClpHelperFunctions.hpp"
10
11#include "ClpInterior.hpp"
12#include "ClpCholeskyWssmp.hpp"
13#include "ClpCholeskyDense.hpp"
14#include "ClpMessage.hpp"
15
16//#############################################################################
17// Constructors / Destructor / Assignment
18//#############################################################################
19
20//-------------------------------------------------------------------
21// Default Constructor
22//-------------------------------------------------------------------
23ClpCholeskyWssmp::ClpCholeskyWssmp (int denseThreshold) 
24  : ClpCholeskyBase(),
25    sparseFactor_(NULL),
26    choleskyStart_(NULL),
27    choleskyRow_(NULL),
28    sizeFactor_(0),
29    rowCopy_(NULL),
30    whichDense_(NULL),
31    denseColumn_(NULL),
32    dense_(NULL),
33    denseThreshold_(denseThreshold)
34{
35  type_=12;
36  memset(integerParameters_,0,64*sizeof(int));
37  memset(doubleParameters_,0,64*sizeof(double));
38}
39
40//-------------------------------------------------------------------
41// Copy constructor
42//-------------------------------------------------------------------
43ClpCholeskyWssmp::ClpCholeskyWssmp (const ClpCholeskyWssmp & rhs) 
44: ClpCholeskyBase(rhs)
45{
46  type_=rhs.type_;
47  sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
48  choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
49  choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
50  sizeFactor_=rhs.sizeFactor_;
51  memcpy(integerParameters_,rhs.integerParameters_,64*sizeof(int));
52  memcpy(doubleParameters_,rhs.doubleParameters_,64*sizeof(double));
53  rowCopy_ = rhs.rowCopy_->clone();
54  whichDense_ = NULL;
55  denseColumn_=NULL;
56  dense_=NULL;
57  denseThreshold_ = rhs.denseThreshold_;
58}
59
60
61//-------------------------------------------------------------------
62// Destructor
63//-------------------------------------------------------------------
64ClpCholeskyWssmp::~ClpCholeskyWssmp ()
65{
66  delete [] sparseFactor_;
67  delete [] choleskyStart_;
68  delete [] choleskyRow_;
69  delete rowCopy_;
70  delete [] whichDense_;
71  delete [] denseColumn_;
72  delete dense_;
73}
74
75//----------------------------------------------------------------
76// Assignment operator
77//-------------------------------------------------------------------
78ClpCholeskyWssmp &
79ClpCholeskyWssmp::operator=(const ClpCholeskyWssmp& rhs)
80{
81  if (this != &rhs) {
82    ClpCholeskyBase::operator=(rhs);
83    delete [] sparseFactor_;
84    delete [] choleskyStart_;
85    delete [] choleskyRow_;
86    delete [] whichDense_;
87    delete [] denseColumn_;
88    delete dense_;
89    sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
90    choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
91    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
92    sizeFactor_=rhs.sizeFactor_;
93    delete rowCopy_;
94    rowCopy_ = rhs.rowCopy_->clone();
95    whichDense_ = NULL;
96    denseColumn_=NULL;
97    dense_=NULL;
98    denseThreshold_ = rhs.denseThreshold_;
99  }
100  return *this;
101}
102//-------------------------------------------------------------------
103// Clone
104//-------------------------------------------------------------------
105ClpCholeskyBase * ClpCholeskyWssmp::clone() const
106{
107  return new ClpCholeskyWssmp(*this);
108}
109// At present I can't get wssmp to work as my libraries seem to be out of sync
110// so I have linked in ekkwssmp which is an older version
111#if 0
112  extern "C" void wssmp(int * n,
113                        int * columnStart , int * rowIndex , double * element,
114                        double * diagonal , int * perm , int * invp ,
115                        double * rhs , int * ldb , int * nrhs ,
116                        double * aux , int * naux ,
117                        int   * mrp , int * iparm , double * dparm);
118#else
119/* minimum needed for user */
120typedef struct EKKModel EKKModel;
121typedef struct EKKContext EKKContext;
122
123
124extern "C"{
125   EKKContext *  ekk_initializeContext();
126   void ekk_endContext(EKKContext * context);
127   EKKModel *  ekk_newModel(EKKContext * env,const char * name);
128   int ekk_deleteModel(EKKModel * model);
129}
130static  EKKModel * model=NULL;
131static  EKKContext * context=NULL;
132extern "C" void ekkwssmp(EKKModel *, int * n,
133                         int * columnStart , int * rowIndex , double * element,
134                         double * diagonal , int * perm , int * invp ,
135                         double * rhs , int * ldb , int * nrhs ,
136                         double * aux , int * naux ,
137                         int   * mrp , int * iparm , double * dparm);
138static void wssmp( int *n, int *ia, int *ja,
139                   double *avals, double *diag, int *perm, int *invp,
140                   double *b, int *ldb, int *nrhs, double *aux, int *
141                   naux, int *mrp, int *iparm, double *dparm)
142{
143  if (!context) {
144    /* initialize OSL environment */
145    context=ekk_initializeContext();
146    model=ekk_newModel(context,"");
147  }
148  ekkwssmp(model,n, ia, ja,
149           avals, diag, perm, invp,
150           b, ldb, nrhs, aux,
151           naux, mrp, iparm, dparm);
152  //ekk_deleteModel(model);
153  //ekk_endContext(context);
154}
155#endif
156/* Orders rows and saves pointer to matrix.and model */
157int 
158ClpCholeskyWssmp::order(ClpInterior * model) 
159{
160  numberRows_ = model->numberRows();
161  rowsDropped_ = new char [numberRows_];
162  memset(rowsDropped_,0,numberRows_);
163  numberRowsDropped_=0;
164  model_=model;
165  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
166  // Space for starts
167  choleskyStart_ = new CoinBigIndex[numberRows_+1];
168  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
169  const int * columnLength = model_->clpMatrix()->getVectorLengths();
170  const int * row = model_->clpMatrix()->getIndices();
171  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
172  const int * rowLength = rowCopy_->getVectorLengths();
173  const int * column = rowCopy_->getIndices();
174  // We need two arrays for counts
175  int * which = new int [numberRows_];
176  int * used = new int[numberRows_+1];
177  CoinZeroN(used,numberRows_);
178  int iRow;
179  sizeFactor_=0;
180  int numberColumns = model->numberColumns();
181  int numberDense=0;
182  if (denseThreshold_>0) {
183    delete [] whichDense_;
184    delete [] denseColumn_;
185    delete dense_;
186    whichDense_ = new char[numberColumns];
187    int iColumn;
188    used[numberRows_]=0;
189    for (iColumn=0;iColumn<numberColumns;iColumn++) {
190      int length = columnLength[iColumn];
191      used[length] += 1;
192    }
193    int nLong=0;
194    int stop = max(denseThreshold_/2,100);
195    for (iRow=numberRows_;iRow>=stop;iRow--) {
196      if (used[iRow]) 
197        printf("%d columns are of length %d\n",used[iRow],iRow);
198      nLong += used[iRow];
199      if (nLong>50||nLong>(numberColumns>>2))
200        break;
201    }
202    CoinZeroN(used,numberRows_);
203    for (iColumn=0;iColumn<numberColumns;iColumn++) {
204      if (columnLength[iColumn]<denseThreshold_) {
205        whichDense_[iColumn]=0;
206      } else {
207        whichDense_[iColumn]=1;
208        numberDense++;
209      }
210    }
211    if (!numberDense||numberDense>100) {
212      // free
213      delete [] whichDense_;
214      whichDense_=NULL;
215      denseColumn_=NULL;
216      dense_=NULL;
217    } else {
218      // space for dense columns
219      denseColumn_ = new double [numberDense*numberRows_];
220      // dense cholesky
221      dense_ = new ClpCholeskyDense();
222      dense_->reserveSpace(numberDense);
223      printf("Taking %d columns as dense\n",numberDense);
224    }
225  }
226  for (iRow=0;iRow<numberRows_;iRow++) {
227    int number=1;
228    // make sure diagonal exists
229    which[0]=iRow;
230    used[iRow]=1;
231    if (!rowsDropped_[iRow]) {
232      CoinBigIndex startRow=rowStart[iRow];
233      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
234      for (CoinBigIndex k=startRow;k<endRow;k++) {
235        int iColumn=column[k];
236        if (!whichDense_||!whichDense_[iColumn]) {
237          CoinBigIndex start=columnStart[iColumn];
238          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
239          for (CoinBigIndex j=start;j<end;j++) {
240            int jRow=row[j];
241            if (jRow>=iRow&&!rowsDropped_[jRow]) {
242              if (!used[jRow]) {
243                used[jRow]=1;
244                which[number++]=jRow;
245              }
246            }
247          }
248        }
249      }
250      sizeFactor_ += number;
251      int j;
252      for (j=0;j<number;j++)
253        used[which[j]]=0;
254    }
255  }
256  delete [] which;
257  // Now we have size - create arrays and fill in
258  choleskyRow_ = new int [sizeFactor_];
259  sparseFactor_ = new double[sizeFactor_];
260  sizeFactor_=0;
261  which = choleskyRow_;
262  for (iRow=0;iRow<numberRows_;iRow++) {
263    int number=1;
264    // make sure diagonal exists
265    which[0]=iRow;
266    used[iRow]=1;
267    choleskyStart_[iRow]=sizeFactor_;
268    if (!rowsDropped_[iRow]) {
269      CoinBigIndex startRow=rowStart[iRow];
270      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
271      for (CoinBigIndex k=startRow;k<endRow;k++) {
272        int iColumn=column[k];
273        if (!whichDense_||!whichDense_[iColumn]) {
274          CoinBigIndex start=columnStart[iColumn];
275          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
276          for (CoinBigIndex j=start;j<end;j++) {
277            int jRow=row[j];
278            if (jRow>=iRow&&!rowsDropped_[jRow]) {
279              if (!used[jRow]) {
280                used[jRow]=1;
281                which[number++]=jRow;
282              }
283            }
284          }
285        }
286      }
287      sizeFactor_ += number;
288      int j;
289      for (j=0;j<number;j++)
290        used[which[j]]=0;
291      // Sort
292      std::sort(which,which+number);
293      // move which on
294      which += number;
295    }
296  }
297  choleskyStart_[numberRows_]=sizeFactor_;
298  delete [] used;
299  permuteIn_ = new int [numberRows_];
300  permuteOut_ = new int[numberRows_];
301  integerParameters_[0]=0;
302  int i0=0;
303  int i1=1;
304  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
305         NULL,permuteOut_,permuteIn_,0,&numberRows_,&i1,
306         NULL,&i0,NULL,integerParameters_,doubleParameters_);
307  integerParameters_[1]=1;//order and symbolic
308  integerParameters_[2]=2;
309  integerParameters_[3]=0;//CSR
310  integerParameters_[4]=0;//C style
311  integerParameters_[13]=1;//reuse initial factorization space
312  integerParameters_[15+0]=1;//ordering
313  integerParameters_[15+1]=0;
314  integerParameters_[15+2]=1;
315  integerParameters_[15+3]=0;
316  integerParameters_[15+4]=1;
317  doubleParameters_[10]=1.0e-20;
318  doubleParameters_[11]=1.0e-15;
319  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
320         NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
321         NULL,&i0,NULL,integerParameters_,doubleParameters_);
322  //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
323  if (integerParameters_[63]) {
324    std::cout<<"wssmp returning error code of "<<integerParameters_[63]<<std::endl;
325    return 1;
326  }
327  std::cout<<integerParameters_[23]<<" elements in sparse Cholesky"<<std::endl;
328  if (!integerParameters_[23]) {
329    for (int iRow=0;iRow<numberRows_;iRow++) {
330      permuteIn_[iRow]=iRow;
331      permuteOut_[iRow]=iRow;
332    }
333    std::cout<<"wssmp says no elements - fully dense? - switching to dense"<<std::endl;
334    integerParameters_[1]=2;
335    integerParameters_[2]=2;
336    integerParameters_[7]=1; // no permute
337    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
338          NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
339          NULL,&i0,NULL,integerParameters_,doubleParameters_);
340    std::cout<<integerParameters_[23]<<" elements in dense Cholesky"<<std::endl;
341  }
342  return 0;
343}
344/* Factorize - filling in rowsDropped and returning number dropped */
345int 
346ClpCholeskyWssmp::factorize(const double * diagonal, int * rowsDropped) 
347{
348  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
349  const int * columnLength = model_->clpMatrix()->getVectorLengths();
350  const int * row = model_->clpMatrix()->getIndices();
351  const double * element = model_->clpMatrix()->getElements();
352  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
353  const int * rowLength = rowCopy_->getVectorLengths();
354  const int * column = rowCopy_->getIndices();
355  const double * elementByRow = rowCopy_->getElements();
356  int numberColumns=model_->clpMatrix()->getNumCols();
357  int iRow;
358  double * work = new double[numberRows_];
359  CoinZeroN(work,numberRows_);
360  const double * diagonalSlack = diagonal + numberColumns;
361  int newDropped=0;
362  double largest;
363  double smallest;
364  int numberDense=0;
365  if (dense_)
366    numberDense=dense_->numberRows();
367  //perturbation
368  double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
369  perturbation=perturbation*perturbation;
370  if (perturbation>1.0) {
371    //if (model_->model()->logLevel()&4)
372      std::cout <<"large perturbation "<<perturbation<<std::endl;
373    perturbation=sqrt(perturbation);;
374    perturbation=1.0;
375  }
376  if (whichDense_) {
377    longDouble * denseBlob = dense_->aMatrix();
378    CoinZeroN(denseBlob,numberDense*numberDense);
379    double * dense = denseColumn_;
380    longDouble * blob = denseBlob;
381    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
382      if (whichDense_[iColumn]) {
383        blob[0]=1.0/diagonal[iColumn];
384        CoinZeroN(dense,numberRows_);
385        CoinBigIndex start=columnStart[iColumn];
386        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
387        for (CoinBigIndex j=start;j<end;j++) {
388          int jRow=row[j];
389          dense[jRow] = element[j];
390        }
391        dense += numberRows_;
392        blob += numberDense+1;
393      }
394    }
395  }
396  for (iRow=0;iRow<numberRows_;iRow++) {
397    double * put = sparseFactor_+choleskyStart_[iRow];
398    int * which = choleskyRow_+choleskyStart_[iRow];
399    int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
400    if (!rowLength[iRow])
401      rowsDropped_[iRow]=1;
402    if (!rowsDropped_[iRow]) {
403      CoinBigIndex startRow=rowStart[iRow];
404      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
405      work[iRow] = diagonalSlack[iRow];
406      for (CoinBigIndex k=startRow;k<endRow;k++) {
407        int iColumn=column[k];
408        if (!whichDense_||!whichDense_[iColumn]) {
409          CoinBigIndex start=columnStart[iColumn];
410          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
411          double multiplier = diagonal[iColumn]*elementByRow[k];
412          for (CoinBigIndex j=start;j<end;j++) {
413            int jRow=row[j];
414            if (jRow>=iRow&&!rowsDropped_[jRow]) {
415              double value=element[j]*multiplier;
416              work[jRow] += value;
417            }
418          }
419        }
420      }
421      int j;
422      for (j=0;j<number;j++) {
423        int jRow = which[j];
424        put[j]=work[jRow];
425        work[jRow]=0.0;
426      }
427    } else {
428      // dropped
429      int j;
430      for (j=1;j<number;j++) {
431        put[j]=0.0;
432      }
433      put[0]=1.0;
434    }
435  }
436  //check sizes
437  double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
438  largest2*=1.0e-19;
439  largest = min (largest2,1.0e-11);
440  int numberDroppedBefore=0;
441  for (iRow=0;iRow<numberRows_;iRow++) {
442    int dropped=rowsDropped_[iRow];
443    // Move to int array
444    rowsDropped[iRow]=dropped;
445    if (!dropped) {
446      CoinBigIndex start = choleskyStart_[iRow];
447      double diagonal = sparseFactor_[start];
448      if (diagonal>largest2) {
449        sparseFactor_[start]=diagonal+perturbation;
450      } else {
451        sparseFactor_[start]=diagonal+perturbation;
452        rowsDropped[iRow]=2;
453        numberDroppedBefore++;
454      } 
455    } 
456  } 
457  int i1=1;
458  int i0=0;
459  integerParameters_[1]=3;
460  integerParameters_[2]=3;
461  integerParameters_[10]=2;
462  //integerParameters_[11]=1;
463  integerParameters_[12]=2;
464  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
465        NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
466        NULL,&i0,rowsDropped,integerParameters_,doubleParameters_);
467  //    NULL,&i0,(int *) diagonal,integerParameters_,doubleParameters_);
468  //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
469  if (integerParameters_[9]) {
470    std::cout<<"scaling applied"<<std::endl;
471  } 
472  newDropped=integerParameters_[20]+numberDroppedBefore;
473  if (integerParameters_[20]) 
474    std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
475  largest=doubleParameters_[3];
476  smallest=doubleParameters_[4];
477  delete [] work;
478  if (model_->messageHandler()->logLevel()>1) 
479    std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
480  choleskyCondition_=largest/smallest;
481  if (whichDense_) {
482    // Update dense columns (just L)
483    // Zero out dropped rows
484    int i;
485    for (i=0;i<numberDense;i++) {
486      double * a = denseColumn_+i*numberRows_;
487      for (int j=0;j<numberRows_;j++) {
488        if (rowsDropped[j])
489          a[j]=0.0;
490      }
491    }
492    integerParameters_[29]=1;
493    int i0=0;
494    integerParameters_[1]=4;
495    integerParameters_[2]=4;
496    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
497       NULL,permuteOut_,permuteIn_,denseColumn_,&numberRows_,&numberDense,
498       NULL,&i0,NULL,integerParameters_,doubleParameters_);
499    integerParameters_[29]=0;
500    dense_->resetRowsDropped();
501    longDouble * denseBlob = dense_->aMatrix();
502    // Update dense matrix
503    for (i=0;i<numberDense;i++) {
504      const double * a = denseColumn_+i*numberRows_;
505      for (int j=0;j<numberDense;j++) {
506        double value = denseBlob[i+j*numberDense];
507        const double * b = denseColumn_+j*numberRows_;
508        for (int k=0;k<numberRows_;k++) 
509          value += a[k]*b[k];
510        denseBlob[i+j*numberDense]=value;
511      }
512    }
513    // dense cholesky (? long double)
514    int * dropped = new int [numberDense];
515    dense_->factorizePart2(dropped);
516    delete [] dropped;
517  }
518  bool cleanCholesky;
519  if (model_->numberIterations()<200) 
520    cleanCholesky=true;
521  else 
522    cleanCholesky=false;
523  if (cleanCholesky) {
524    //drop fresh makes some formADAT easier
525    int oldDropped=numberRowsDropped_;
526    if (newDropped||numberRowsDropped_) {
527      std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
528          newDropped<<" dropped)";
529      if (newDropped>oldDropped) 
530        std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
531      std::cout<<std::endl;
532      newDropped=0;
533      for (int i=0;i<numberRows_;i++) {
534        char dropped = rowsDropped[i];
535        rowsDropped_[i]=dropped;
536        if (dropped==2) {
537          //dropped this time
538          rowsDropped[newDropped++]=i;
539          rowsDropped_[i]=0;
540        } 
541      } 
542      numberRowsDropped_=newDropped;
543      newDropped=-(1+newDropped);
544    } 
545  } else {
546    if (newDropped) {
547      newDropped=0;
548      for (int i=0;i<numberRows_;i++) {
549        char dropped = rowsDropped[i];
550        rowsDropped_[i]=dropped;
551        if (dropped==2) {
552          //dropped this time
553          rowsDropped[newDropped++]=i;
554          rowsDropped_[i]=1;
555        } 
556      } 
557    } 
558    numberRowsDropped_+=newDropped;
559    if (numberRowsDropped_) {
560      std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
561          numberRowsDropped_<<" dropped)";
562      if (newDropped) {
563        std::cout<<" ( "<<newDropped<<" dropped this time)";
564      } 
565      std::cout<<std::endl;
566    } 
567  }
568  status_=0;
569  return newDropped;
570}
571/* Uses factorization to solve. */
572void 
573ClpCholeskyWssmp::solve (double * region) 
574{
575  int i1=1;
576  int i0=0;
577  integerParameters_[1]=4;
578  integerParameters_[2]=4;
579#if 0
580  integerParameters_[5]=3;
581  doubleParameters_[5]=1.0e-10;
582  integerParameters_[6]=6;
583#endif
584  if (!whichDense_) {
585    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
586          NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
587          NULL,&i0,NULL,integerParameters_,doubleParameters_);
588  } else {
589    // dense columns
590    integerParameters_[29]=1;
591    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
592          NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
593          NULL,&i0,NULL,integerParameters_,doubleParameters_);
594    // do change;
595    int numberDense = dense_->numberRows();
596    double * change = new double[numberDense];
597    int i;
598    for (i=0;i<numberDense;i++) {
599      const double * a = denseColumn_+i*numberRows_;
600      double value =0.0;
601      for (int iRow=0;iRow<numberRows_;iRow++) 
602        value += a[iRow]*region[iRow];
603      change[i]=value;
604    }
605    // solve
606    dense_->solve(change);
607    for (i=0;i<numberDense;i++) {
608      const double * a = denseColumn_+i*numberRows_;
609      double value = change[i];
610      for (int iRow=0;iRow<numberRows_;iRow++) 
611        region[iRow] -= value*a[iRow];
612    }
613    delete [] change;
614    // and finish off
615    integerParameters_[29]=2;
616    integerParameters_[1]=4;
617    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
618          NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
619          NULL,&i0,NULL,integerParameters_,doubleParameters_);
620    integerParameters_[29]=0;
621  }
622#if 0
623  if (integerParameters_[5]) {
624    std::cout<<integerParameters_[5]<<" refinements ";
625  }
626  std::cout<<doubleParameters_[6]<<std::endl;
627#endif
628}
629#endif
Note: See TracBrowser for help on using the repository browser.