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

Last change on this file since 994 was 994, checked in by andreasw, 13 years ago

reran autotools; added support for WSMP (--with-wsmp as configure option)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.2 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#include "ClpHelperFunctions.hpp"
9
10#include "ClpInterior.hpp"
11#include "ClpCholeskyWssmp.hpp"
12#include "ClpCholeskyDense.hpp"
13#include "ClpMessage.hpp"
14
15#ifdef WSSMP_BARRIER
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#ifdef COIN_DEVELOP
385    //if (model_->model()->logLevel()&4)
386      std::cout <<"large perturbation "<<perturbation<<std::endl;
387#endif
388    perturbation=sqrt(perturbation);;
389    perturbation=1.0;
390  }
391  if (whichDense_) {
392    double * denseDiagonal = dense_->diagonal();
393    double * dense = denseColumn_;
394    int iDense=0;
395    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
396      if (whichDense_[iColumn]) {
397        denseDiagonal[iDense++]=1.0/diagonal[iColumn];
398        CoinZeroN(dense,numberRows_);
399        CoinBigIndex start=columnStart[iColumn];
400        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
401        for (CoinBigIndex j=start;j<end;j++) {
402          int jRow=row[j];
403          dense[jRow] = element[j];
404        }
405        dense += numberRows_;
406      }
407    }
408  }
409  double delta2 = model_->delta(); // add delta*delta to diagonal
410  delta2 *= delta2;
411  for (iRow=0;iRow<numberRows_;iRow++) {
412    double * put = sparseFactor_+choleskyStart_[iRow];
413    int * which = choleskyRow_+choleskyStart_[iRow];
414    int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
415    if (!rowLength[iRow])
416      rowsDropped_[iRow]=1;
417    if (!rowsDropped_[iRow]) {
418      CoinBigIndex startRow=rowStart[iRow];
419      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
420      work[iRow] = diagonalSlack[iRow]+delta2;
421      for (CoinBigIndex k=startRow;k<endRow;k++) {
422        int iColumn=column[k];
423        if (!whichDense_||!whichDense_[iColumn]) {
424          CoinBigIndex start=columnStart[iColumn];
425          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
426          double multiplier = diagonal[iColumn]*elementByRow[k];
427          for (CoinBigIndex j=start;j<end;j++) {
428            int jRow=row[j];
429            if (jRow>=iRow&&!rowsDropped_[jRow]) {
430              double value=element[j]*multiplier;
431              work[jRow] += value;
432            }
433          }
434        }
435      }
436      int j;
437      for (j=0;j<number;j++) {
438        int jRow = which[j];
439        put[j]=work[jRow];
440        work[jRow]=0.0;
441      }
442    } else {
443      // dropped
444      int j;
445      for (j=1;j<number;j++) {
446        put[j]=0.0;
447      }
448      put[0]=1.0;
449    }
450  }
451  //check sizes
452  double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
453  largest2*=1.0e-20;
454  largest = CoinMin(largest2,1.0e-11);
455  int numberDroppedBefore=0;
456  for (iRow=0;iRow<numberRows_;iRow++) {
457    int dropped=rowsDropped_[iRow];
458    // Move to int array
459    rowsDropped[iRow]=dropped;
460    if (!dropped) {
461      CoinBigIndex start = choleskyStart_[iRow];
462      double diagonal = sparseFactor_[start];
463      if (diagonal>largest2) {
464        sparseFactor_[start]=diagonal+perturbation;
465      } else {
466        sparseFactor_[start]=diagonal+perturbation;
467        rowsDropped[iRow]=2;
468        numberDroppedBefore++;
469      } 
470    } 
471  } 
472  int i1=1;
473  int i0=0;
474  integerParameters_[1]=3;
475  integerParameters_[2]=3;
476  integerParameters_[10]=2;
477  //integerParameters_[11]=1;
478  integerParameters_[12]=2;
479  doubleParameters_[10]=CoinMax(1.0e-20,largest);
480  if (doubleParameters_[10]>1.0e-3)
481    integerParameters_[9]=1;
482  else
483    integerParameters_[9]=0;
484  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
485        NULL,permute_,permuteInverse_,NULL,&numberRows_,&i1,
486        NULL,&i0,rowsDropped,integerParameters_,doubleParameters_);
487  //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
488  if (integerParameters_[9]) {
489    std::cout<<"scaling applied"<<std::endl;
490  } 
491  newDropped=integerParameters_[20]+numberDroppedBefore;
492  //if (integerParameters_[20])
493  //std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
494  largest=doubleParameters_[3];
495  smallest=doubleParameters_[4];
496  delete [] work;
497  if (model_->messageHandler()->logLevel()>1) 
498    std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
499  choleskyCondition_=largest/smallest;
500  if (integerParameters_[63]<0)
501    return -1; // out of memory
502  if (whichDense_) {
503    // Update dense columns (just L)
504    // Zero out dropped rows
505    int i;
506    for (i=0;i<numberDense;i++) {
507      double * a = denseColumn_+i*numberRows_;
508      for (int j=0;j<numberRows_;j++) {
509        if (rowsDropped[j])
510          a[j]=0.0;
511      }
512    }
513    integerParameters_[29]=1;
514    int i0=0;
515    integerParameters_[1]=4;
516    integerParameters_[2]=4;
517    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
518       NULL,permute_,permuteInverse_,denseColumn_,&numberRows_,&numberDense,
519       NULL,&i0,NULL,integerParameters_,doubleParameters_);
520    integerParameters_[29]=0;
521    dense_->resetRowsDropped();
522    longDouble * denseBlob = dense_->aMatrix();
523    double * denseDiagonal = dense_->diagonal();
524    // Update dense matrix
525    for (i=0;i<numberDense;i++) {
526      const double * a = denseColumn_+i*numberRows_;
527      // do diagonal
528      double value = denseDiagonal[i];
529      const double * b = denseColumn_+i*numberRows_;
530      for (int k=0;k<numberRows_;k++) 
531        value += a[k]*b[k];
532      denseDiagonal[i]=value;
533      for (int j=i+1;j<numberDense;j++) {
534        double value = 0.0;
535        const double * b = denseColumn_+j*numberRows_;
536        for (int k=0;k<numberRows_;k++) 
537          value += a[k]*b[k];
538        *denseBlob=value;
539        denseBlob++;
540      }
541    }
542    // dense cholesky (? long double)
543    int * dropped = new int [numberDense];
544    dense_->factorizePart2(dropped);
545    delete [] dropped;
546  }
547  bool cleanCholesky;
548  if (model_->numberIterations()<2000) 
549    cleanCholesky=true;
550  else 
551    cleanCholesky=false;
552  if (cleanCholesky) {
553    //drop fresh makes some formADAT easier
554    //int oldDropped=numberRowsDropped_;
555    if (newDropped||numberRowsDropped_) {
556      //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
557      //  newDropped<<" dropped)";
558      //if (newDropped>oldDropped)
559      //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
560      //std::cout<<std::endl;
561      newDropped=0;
562      for (int i=0;i<numberRows_;i++) {
563        char dropped = rowsDropped[i];
564        rowsDropped_[i]=dropped;
565        if (dropped==2) {
566          //dropped this time
567          rowsDropped[newDropped++]=i;
568          rowsDropped_[i]=0;
569        } 
570      } 
571      numberRowsDropped_=newDropped;
572      newDropped=-(2+newDropped);
573    } 
574  } else {
575    if (newDropped) {
576      newDropped=0;
577      for (int i=0;i<numberRows_;i++) {
578        char dropped = rowsDropped[i];
579        rowsDropped_[i]=dropped;
580        if (dropped==2) {
581          //dropped this time
582          rowsDropped[newDropped++]=i;
583          rowsDropped_[i]=1;
584        } 
585      } 
586    } 
587    numberRowsDropped_+=newDropped;
588    if (numberRowsDropped_&&0) {
589      std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
590          numberRowsDropped_<<" dropped)";
591      if (newDropped) {
592        std::cout<<" ( "<<newDropped<<" dropped this time)";
593      } 
594      std::cout<<std::endl;
595    } 
596  }
597  status_=0;
598  return newDropped;
599}
600/* Uses factorization to solve. */
601void 
602ClpCholeskyWssmp::solve (double * region) 
603{
604  int i1=1;
605  int i0=0;
606  integerParameters_[1]=4;
607  integerParameters_[2]=4;
608#if 0
609  integerParameters_[5]=3;
610  doubleParameters_[5]=1.0e-10;
611  integerParameters_[6]=6;
612#endif
613  if (!whichDense_) {
614    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
615          NULL,permute_,permuteInverse_,region,&numberRows_,&i1,
616          NULL,&i0,NULL,integerParameters_,doubleParameters_);
617  } else {
618    // dense columns
619    integerParameters_[29]=1;
620    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
621          NULL,permute_,permuteInverse_,region,&numberRows_,&i1,
622          NULL,&i0,NULL,integerParameters_,doubleParameters_);
623    // do change;
624    int numberDense = dense_->numberRows();
625    double * change = new double[numberDense];
626    int i;
627    for (i=0;i<numberDense;i++) {
628      const double * a = denseColumn_+i*numberRows_;
629      double value =0.0;
630      for (int iRow=0;iRow<numberRows_;iRow++) 
631        value += a[iRow]*region[iRow];
632      change[i]=value;
633    }
634    // solve
635    dense_->solve(change);
636    for (i=0;i<numberDense;i++) {
637      const double * a = denseColumn_+i*numberRows_;
638      double value = change[i];
639      for (int iRow=0;iRow<numberRows_;iRow++) 
640        region[iRow] -= value*a[iRow];
641    }
642    delete [] change;
643    // and finish off
644    integerParameters_[29]=2;
645    integerParameters_[1]=4;
646    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
647          NULL,permute_,permuteInverse_,region,&numberRows_,&i1,
648          NULL,&i0,NULL,integerParameters_,doubleParameters_);
649    integerParameters_[29]=0;
650  }
651#if 0
652  if (integerParameters_[5]) {
653    std::cout<<integerParameters_[5]<<" refinements ";
654  }
655  std::cout<<doubleParameters_[6]<<std::endl;
656#endif
657}
658#endif
Note: See TracBrowser for help on using the repository browser.