source: trunk/Clp/src/ClpCholeskyWssmp.cpp @ 754

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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