source: trunk/ClpCholeskyWssmp.cpp @ 328

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

For dynamic matrices

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.1 KB
RevLine 
[313]1#ifdef REAL_BARRIER
[266]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"
[320]13#include "ClpCholeskyDense.hpp"
[266]14#include "ClpMessage.hpp"
15
16//#############################################################################
17// Constructors / Destructor / Assignment
18//#############################################################################
19
20//-------------------------------------------------------------------
21// Default Constructor
22//-------------------------------------------------------------------
[320]23ClpCholeskyWssmp::ClpCholeskyWssmp (int denseThreshold) 
[266]24  : ClpCholeskyBase(),
25    sparseFactor_(NULL),
26    choleskyStart_(NULL),
27    choleskyRow_(NULL),
28    sizeFactor_(0),
[320]29    rowCopy_(NULL),
30    whichDense_(NULL),
31    denseColumn_(NULL),
32    dense_(NULL),
33    denseThreshold_(denseThreshold)
[266]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();
[320]54  whichDense_ = NULL;
55  denseColumn_=NULL;
56  dense_=NULL;
57  denseThreshold_ = rhs.denseThreshold_;
[266]58}
59
60
61//-------------------------------------------------------------------
62// Destructor
63//-------------------------------------------------------------------
64ClpCholeskyWssmp::~ClpCholeskyWssmp ()
65{
66  delete [] sparseFactor_;
67  delete [] choleskyStart_;
68  delete [] choleskyRow_;
69  delete rowCopy_;
[320]70  delete [] whichDense_;
71  delete [] denseColumn_;
72  delete dense_;
[266]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_;
[320]86    delete [] whichDense_;
87    delete [] denseColumn_;
88    delete dense_;
[266]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();
[320]95    whichDense_ = NULL;
96    denseColumn_=NULL;
97    dense_=NULL;
98    denseThreshold_ = rhs.denseThreshold_;
[266]99  }
100  return *this;
101}
102//-------------------------------------------------------------------
103// Clone
104//-------------------------------------------------------------------
105ClpCholeskyBase * ClpCholeskyWssmp::clone() const
106{
107  return new ClpCholeskyWssmp(*this);
108}
[275]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
[266]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);
[275]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
[266]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];
[275]168  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
169  const int * columnLength = model_->clpMatrix()->getVectorLengths();
170  const int * row = model_->clpMatrix()->getIndices();
[266]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_];
[320]176  int * used = new int[numberRows_+1];
[266]177  CoinZeroN(used,numberRows_);
178  int iRow;
179  sizeFactor_=0;
[320]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  }
[266]226  for (iRow=0;iRow<numberRows_;iRow++) {
[312]227    int number=1;
228    // make sure diagonal exists
229    which[0]=iRow;
230    used[iRow]=1;
[266]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];
[320]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              }
[266]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
[328]258  try { 
259    choleskyRow_ = new int [sizeFactor_];
260  }
261  catch (...) {
262    // no memory
263    delete [] choleskyStart_;
264    choleskyStart_=NULL;
265    return -1;
266  }
267  try {
268    sparseFactor_ = new double[sizeFactor_];
269  }
270  catch (...) {
271    // no memory
272    delete [] choleskyRow_;
273    choleskyRow_=NULL;
274    delete [] choleskyStart_;
275    choleskyStart_=NULL;
276    return -1;
277  }
278 
[266]279  sizeFactor_=0;
280  which = choleskyRow_;
281  for (iRow=0;iRow<numberRows_;iRow++) {
[312]282    int number=1;
283    // make sure diagonal exists
284    which[0]=iRow;
285    used[iRow]=1;
[266]286    choleskyStart_[iRow]=sizeFactor_;
287    if (!rowsDropped_[iRow]) {
288      CoinBigIndex startRow=rowStart[iRow];
289      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
290      for (CoinBigIndex k=startRow;k<endRow;k++) {
291        int iColumn=column[k];
[320]292        if (!whichDense_||!whichDense_[iColumn]) {
293          CoinBigIndex start=columnStart[iColumn];
294          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
295          for (CoinBigIndex j=start;j<end;j++) {
296            int jRow=row[j];
297            if (jRow>=iRow&&!rowsDropped_[jRow]) {
298              if (!used[jRow]) {
299                used[jRow]=1;
300                which[number++]=jRow;
301              }
[266]302            }
303          }
304        }
305      }
306      sizeFactor_ += number;
307      int j;
308      for (j=0;j<number;j++)
309        used[which[j]]=0;
[275]310      // Sort
311      std::sort(which,which+number);
[266]312      // move which on
313      which += number;
314    }
315  }
316  choleskyStart_[numberRows_]=sizeFactor_;
317  delete [] used;
318  permuteIn_ = new int [numberRows_];
319  permuteOut_ = new int[numberRows_];
320  integerParameters_[0]=0;
321  int i0=0;
322  int i1=1;
323  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
324         NULL,permuteOut_,permuteIn_,0,&numberRows_,&i1,
325         NULL,&i0,NULL,integerParameters_,doubleParameters_);
326  integerParameters_[1]=1;//order and symbolic
327  integerParameters_[2]=2;
328  integerParameters_[3]=0;//CSR
329  integerParameters_[4]=0;//C style
330  integerParameters_[13]=1;//reuse initial factorization space
331  integerParameters_[15+0]=1;//ordering
332  integerParameters_[15+1]=0;
333  integerParameters_[15+2]=1;
334  integerParameters_[15+3]=0;
335  integerParameters_[15+4]=1;
336  doubleParameters_[10]=1.0e-20;
337  doubleParameters_[11]=1.0e-15;
338  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
339         NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
340         NULL,&i0,NULL,integerParameters_,doubleParameters_);
[275]341  //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
342  if (integerParameters_[63]) {
343    std::cout<<"wssmp returning error code of "<<integerParameters_[63]<<std::endl;
[315]344    return 1;
[275]345  }
[266]346  std::cout<<integerParameters_[23]<<" elements in sparse Cholesky"<<std::endl;
[275]347  if (!integerParameters_[23]) {
[316]348    for (int iRow=0;iRow<numberRows_;iRow++) {
349      permuteIn_[iRow]=iRow;
350      permuteOut_[iRow]=iRow;
351    }
[275]352    std::cout<<"wssmp says no elements - fully dense? - switching to dense"<<std::endl;
353    integerParameters_[1]=2;
354    integerParameters_[2]=2;
355    integerParameters_[7]=1; // no permute
356    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
357          NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
358          NULL,&i0,NULL,integerParameters_,doubleParameters_);
359    std::cout<<integerParameters_[23]<<" elements in dense Cholesky"<<std::endl;
360  }
[266]361  return 0;
362}
363/* Factorize - filling in rowsDropped and returning number dropped */
364int 
365ClpCholeskyWssmp::factorize(const double * diagonal, int * rowsDropped) 
366{
[275]367  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
368  const int * columnLength = model_->clpMatrix()->getVectorLengths();
369  const int * row = model_->clpMatrix()->getIndices();
370  const double * element = model_->clpMatrix()->getElements();
[266]371  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
372  const int * rowLength = rowCopy_->getVectorLengths();
373  const int * column = rowCopy_->getIndices();
374  const double * elementByRow = rowCopy_->getElements();
[275]375  int numberColumns=model_->clpMatrix()->getNumCols();
[266]376  int iRow;
377  double * work = new double[numberRows_];
378  CoinZeroN(work,numberRows_);
379  const double * diagonalSlack = diagonal + numberColumns;
380  int newDropped=0;
381  double largest;
382  double smallest;
[320]383  int numberDense=0;
384  if (dense_)
385    numberDense=dense_->numberRows();
[266]386  //perturbation
387  double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
388  perturbation=perturbation*perturbation;
389  if (perturbation>1.0) {
390    //if (model_->model()->logLevel()&4)
391      std::cout <<"large perturbation "<<perturbation<<std::endl;
392    perturbation=sqrt(perturbation);;
393    perturbation=1.0;
[320]394  }
395  if (whichDense_) {
396    longDouble * denseBlob = dense_->aMatrix();
397    CoinZeroN(denseBlob,numberDense*numberDense);
398    double * dense = denseColumn_;
399    longDouble * blob = denseBlob;
400    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
401      if (whichDense_[iColumn]) {
402        blob[0]=1.0/diagonal[iColumn];
403        CoinZeroN(dense,numberRows_);
404        CoinBigIndex start=columnStart[iColumn];
405        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
406        for (CoinBigIndex j=start;j<end;j++) {
407          int jRow=row[j];
408          dense[jRow] = element[j];
409        }
410        dense += numberRows_;
411        blob += numberDense+1;
412      }
413    }
414  }
[266]415  for (iRow=0;iRow<numberRows_;iRow++) {
416    double * put = sparseFactor_+choleskyStart_[iRow];
417    int * which = choleskyRow_+choleskyStart_[iRow];
418    int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
[312]419    if (!rowLength[iRow])
420      rowsDropped_[iRow]=1;
[266]421    if (!rowsDropped_[iRow]) {
422      CoinBigIndex startRow=rowStart[iRow];
423      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
424      work[iRow] = diagonalSlack[iRow];
425      for (CoinBigIndex k=startRow;k<endRow;k++) {
426        int iColumn=column[k];
[320]427        if (!whichDense_||!whichDense_[iColumn]) {
428          CoinBigIndex start=columnStart[iColumn];
429          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
430          double multiplier = diagonal[iColumn]*elementByRow[k];
431          for (CoinBigIndex j=start;j<end;j++) {
432            int jRow=row[j];
433            if (jRow>=iRow&&!rowsDropped_[jRow]) {
434              double value=element[j]*multiplier;
435              work[jRow] += value;
436            }
[266]437          }
438        }
439      }
440      int j;
441      for (j=0;j<number;j++) {
442        int jRow = which[j];
443        put[j]=work[jRow];
444        work[jRow]=0.0;
445      }
[275]446    } else {
447      // dropped
448      int j;
[286]449      for (j=1;j<number;j++) {
[275]450        put[j]=0.0;
451      }
[286]452      put[0]=1.0;
[266]453    }
454  }
455  //check sizes
456  double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
457  largest2*=1.0e-19;
458  largest = min (largest2,1.0e-11);
459  int numberDroppedBefore=0;
[284]460  for (iRow=0;iRow<numberRows_;iRow++) {
[266]461    int dropped=rowsDropped_[iRow];
462    // Move to int array
463    rowsDropped[iRow]=dropped;
464    if (!dropped) {
465      CoinBigIndex start = choleskyStart_[iRow];
466      double diagonal = sparseFactor_[start];
467      if (diagonal>largest2) {
468        sparseFactor_[start]=diagonal+perturbation;
469      } else {
470        sparseFactor_[start]=diagonal+perturbation;
471        rowsDropped[iRow]=2;
472        numberDroppedBefore++;
473      } 
474    } 
475  } 
476  int i1=1;
477  int i0=0;
478  integerParameters_[1]=3;
479  integerParameters_[2]=3;
480  integerParameters_[10]=2;
481  //integerParameters_[11]=1;
482  integerParameters_[12]=2;
483  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
484        NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
485        NULL,&i0,rowsDropped,integerParameters_,doubleParameters_);
[275]486  //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
[266]487  if (integerParameters_[9]) {
488    std::cout<<"scaling applied"<<std::endl;
489  } 
490  newDropped=integerParameters_[20]+numberDroppedBefore;
[328]491  //if (integerParameters_[20])
492  //std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
[266]493  largest=doubleParameters_[3];
494  smallest=doubleParameters_[4];
495  delete [] work;
[275]496  if (model_->messageHandler()->logLevel()>1) 
[266]497    std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
498  choleskyCondition_=largest/smallest;
[328]499  if (integerParameters_[63]<0)
500    return -1; // out of memory
[320]501  if (whichDense_) {
502    // Update dense columns (just L)
503    // Zero out dropped rows
504    int i;
505    for (i=0;i<numberDense;i++) {
506      double * a = denseColumn_+i*numberRows_;
507      for (int j=0;j<numberRows_;j++) {
508        if (rowsDropped[j])
509          a[j]=0.0;
510      }
511    }
512    integerParameters_[29]=1;
513    int i0=0;
514    integerParameters_[1]=4;
515    integerParameters_[2]=4;
516    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
517       NULL,permuteOut_,permuteIn_,denseColumn_,&numberRows_,&numberDense,
518       NULL,&i0,NULL,integerParameters_,doubleParameters_);
519    integerParameters_[29]=0;
520    dense_->resetRowsDropped();
521    longDouble * denseBlob = dense_->aMatrix();
522    // Update dense matrix
523    for (i=0;i<numberDense;i++) {
524      const double * a = denseColumn_+i*numberRows_;
525      for (int j=0;j<numberDense;j++) {
526        double value = denseBlob[i+j*numberDense];
527        const double * b = denseColumn_+j*numberRows_;
528        for (int k=0;k<numberRows_;k++) 
529          value += a[k]*b[k];
530        denseBlob[i+j*numberDense]=value;
531      }
532    }
533    // dense cholesky (? long double)
534    int * dropped = new int [numberDense];
535    dense_->factorizePart2(dropped);
536    delete [] dropped;
537  }
[266]538  bool cleanCholesky;
[286]539  if (model_->numberIterations()<200) 
[266]540    cleanCholesky=true;
541  else 
542    cleanCholesky=false;
543  if (cleanCholesky) {
544    //drop fresh makes some formADAT easier
[328]545    //int oldDropped=numberRowsDropped_;
[266]546    if (newDropped||numberRowsDropped_) {
[328]547      //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
548      //  newDropped<<" dropped)";
549      //if (newDropped>oldDropped)
550      //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
551      //std::cout<<std::endl;
[266]552      newDropped=0;
553      for (int i=0;i<numberRows_;i++) {
554        char dropped = rowsDropped[i];
555        rowsDropped_[i]=dropped;
556        if (dropped==2) {
557          //dropped this time
558          rowsDropped[newDropped++]=i;
559          rowsDropped_[i]=0;
560        } 
561      } 
562      numberRowsDropped_=newDropped;
[328]563      newDropped=-(2+newDropped);
[266]564    } 
565  } else {
566    if (newDropped) {
567      newDropped=0;
568      for (int i=0;i<numberRows_;i++) {
569        char dropped = rowsDropped[i];
570        rowsDropped_[i]=dropped;
571        if (dropped==2) {
572          //dropped this time
573          rowsDropped[newDropped++]=i;
574          rowsDropped_[i]=1;
575        } 
576      } 
577    } 
578    numberRowsDropped_+=newDropped;
[328]579    if (numberRowsDropped_&&0) {
[266]580      std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
581          numberRowsDropped_<<" dropped)";
582      if (newDropped) {
583        std::cout<<" ( "<<newDropped<<" dropped this time)";
584      } 
585      std::cout<<std::endl;
586    } 
[320]587  }
[266]588  status_=0;
589  return newDropped;
590}
591/* Uses factorization to solve. */
592void 
593ClpCholeskyWssmp::solve (double * region) 
594{
595  int i1=1;
596  int i0=0;
597  integerParameters_[1]=4;
598  integerParameters_[2]=4;
599#if 0
600  integerParameters_[5]=3;
601  doubleParameters_[5]=1.0e-10;
602  integerParameters_[6]=6;
603#endif
[320]604  if (!whichDense_) {
605    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
606          NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
607          NULL,&i0,NULL,integerParameters_,doubleParameters_);
608  } else {
609    // dense columns
610    integerParameters_[29]=1;
611    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
612          NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
613          NULL,&i0,NULL,integerParameters_,doubleParameters_);
614    // do change;
615    int numberDense = dense_->numberRows();
616    double * change = new double[numberDense];
617    int i;
618    for (i=0;i<numberDense;i++) {
619      const double * a = denseColumn_+i*numberRows_;
620      double value =0.0;
621      for (int iRow=0;iRow<numberRows_;iRow++) 
622        value += a[iRow]*region[iRow];
623      change[i]=value;
624    }
625    // solve
626    dense_->solve(change);
627    for (i=0;i<numberDense;i++) {
628      const double * a = denseColumn_+i*numberRows_;
629      double value = change[i];
630      for (int iRow=0;iRow<numberRows_;iRow++) 
631        region[iRow] -= value*a[iRow];
632    }
633    delete [] change;
634    // and finish off
635    integerParameters_[29]=2;
636    integerParameters_[1]=4;
637    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
638          NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
639          NULL,&i0,NULL,integerParameters_,doubleParameters_);
640    integerParameters_[29]=0;
641  }
[266]642#if 0
643  if (integerParameters_[5]) {
644    std::cout<<integerParameters_[5]<<" refinements ";
645  }
646  std::cout<<doubleParameters_[6]<<std::endl;
647#endif
648}
[313]649#endif
Note: See TracBrowser for help on using the repository browser.