source: trunk/ClpCholeskyWssmp.cpp @ 284

Last change on this file since 284 was 284, checked in by jpfasano, 16 years ago

Modified to compile with MS Visual Studio Version 6

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.5 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 "ClpMessage.hpp"
13
14//#############################################################################
15// Constructors / Destructor / Assignment
16//#############################################################################
17
18//-------------------------------------------------------------------
19// Default Constructor
20//-------------------------------------------------------------------
21ClpCholeskyWssmp::ClpCholeskyWssmp () 
22  : ClpCholeskyBase(),
23    sparseFactor_(NULL),
24    choleskyStart_(NULL),
25    choleskyRow_(NULL),
26    sizeFactor_(0),
27    rowCopy_(NULL)
28{
29  type_=12;
30  memset(integerParameters_,0,64*sizeof(int));
31  memset(doubleParameters_,0,64*sizeof(double));
32}
33
34//-------------------------------------------------------------------
35// Copy constructor
36//-------------------------------------------------------------------
37ClpCholeskyWssmp::ClpCholeskyWssmp (const ClpCholeskyWssmp & rhs) 
38: ClpCholeskyBase(rhs)
39{
40  type_=rhs.type_;
41  sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
42  choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
43  choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
44  sizeFactor_=rhs.sizeFactor_;
45  memcpy(integerParameters_,rhs.integerParameters_,64*sizeof(int));
46  memcpy(doubleParameters_,rhs.doubleParameters_,64*sizeof(double));
47  rowCopy_ = rhs.rowCopy_->clone();
48}
49
50
51//-------------------------------------------------------------------
52// Destructor
53//-------------------------------------------------------------------
54ClpCholeskyWssmp::~ClpCholeskyWssmp ()
55{
56  delete [] sparseFactor_;
57  delete [] choleskyStart_;
58  delete [] choleskyRow_;
59  delete rowCopy_;
60}
61
62//----------------------------------------------------------------
63// Assignment operator
64//-------------------------------------------------------------------
65ClpCholeskyWssmp &
66ClpCholeskyWssmp::operator=(const ClpCholeskyWssmp& rhs)
67{
68  if (this != &rhs) {
69    ClpCholeskyBase::operator=(rhs);
70    delete [] sparseFactor_;
71    delete [] choleskyStart_;
72    delete [] choleskyRow_;
73    sparseFactor_ = ClpCopyOfArray(rhs.sparseFactor_,rhs.sizeFactor_);
74    choleskyStart_ = ClpCopyOfArray(rhs.choleskyStart_,numberRows_+1);
75    choleskyRow_ = ClpCopyOfArray(rhs.choleskyRow_,rhs.sizeFactor_);
76    sizeFactor_=rhs.sizeFactor_;
77    delete rowCopy_;
78    rowCopy_ = rhs.rowCopy_->clone();
79  }
80  return *this;
81}
82//-------------------------------------------------------------------
83// Clone
84//-------------------------------------------------------------------
85ClpCholeskyBase * ClpCholeskyWssmp::clone() const
86{
87  return new ClpCholeskyWssmp(*this);
88}
89// At present I can't get wssmp to work as my libraries seem to be out of sync
90// so I have linked in ekkwssmp which is an older version
91#if 0
92  extern "C" void wssmp(int * n,
93                        int * columnStart , int * rowIndex , double * element,
94                        double * diagonal , int * perm , int * invp ,
95                        double * rhs , int * ldb , int * nrhs ,
96                        double * aux , int * naux ,
97                        int   * mrp , int * iparm , double * dparm);
98#else
99/* minimum needed for user */
100typedef struct EKKModel EKKModel;
101typedef struct EKKContext EKKContext;
102
103
104extern "C"{
105   EKKContext *  ekk_initializeContext();
106   void ekk_endContext(EKKContext * context);
107   EKKModel *  ekk_newModel(EKKContext * env,const char * name);
108   int ekk_deleteModel(EKKModel * model);
109}
110static  EKKModel * model=NULL;
111static  EKKContext * context=NULL;
112extern "C" void ekkwssmp(EKKModel *, 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);
118static void wssmp( int *n, int *ia, int *ja,
119                   double *avals, double *diag, int *perm, int *invp,
120                   double *b, int *ldb, int *nrhs, double *aux, int *
121                   naux, int *mrp, int *iparm, double *dparm)
122{
123  if (!context) {
124    /* initialize OSL environment */
125    context=ekk_initializeContext();
126    model=ekk_newModel(context,"");
127  }
128  ekkwssmp(model,n, ia, ja,
129           avals, diag, perm, invp,
130           b, ldb, nrhs, aux,
131           naux, mrp, iparm, dparm);
132  //ekk_deleteModel(model);
133  //ekk_endContext(context);
134}
135#endif
136/* Orders rows and saves pointer to matrix.and model */
137int 
138ClpCholeskyWssmp::order(ClpInterior * model) 
139{
140  numberRows_ = model->numberRows();
141  rowsDropped_ = new char [numberRows_];
142  memset(rowsDropped_,0,numberRows_);
143  numberRowsDropped_=0;
144  model_=model;
145  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
146  // Space for starts
147  choleskyStart_ = new CoinBigIndex[numberRows_+1];
148  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
149  const int * columnLength = model_->clpMatrix()->getVectorLengths();
150  const int * row = model_->clpMatrix()->getIndices();
151  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
152  const int * rowLength = rowCopy_->getVectorLengths();
153  const int * column = rowCopy_->getIndices();
154  // We need two arrays for counts
155  int * which = new int [numberRows_];
156  int * used = new int[numberRows_];
157  CoinZeroN(used,numberRows_);
158  int iRow;
159  sizeFactor_=0;
160  for (iRow=0;iRow<numberRows_;iRow++) {
161    int number=0;
162    if (!rowsDropped_[iRow]) {
163      CoinBigIndex startRow=rowStart[iRow];
164      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
165      for (CoinBigIndex k=startRow;k<endRow;k++) {
166        int iColumn=column[k];
167        CoinBigIndex start=columnStart[iColumn];
168        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
169        for (CoinBigIndex j=start;j<end;j++) {
170          int jRow=row[j];
171          if (jRow>=iRow&&!rowsDropped_[jRow]) {
172            if (!used[jRow]) {
173              used[jRow]=1;
174              which[number++]=jRow;
175            }
176          }
177        }
178      }
179      sizeFactor_ += number;
180      int j;
181      for (j=0;j<number;j++)
182        used[which[j]]=0;
183    }
184  }
185  delete [] which;
186  // Now we have size - create arrays and fill in
187  choleskyRow_ = new int [sizeFactor_];
188  sparseFactor_ = new double[sizeFactor_];
189  sizeFactor_=0;
190  which = choleskyRow_;
191  for (iRow=0;iRow<numberRows_;iRow++) {
192    int number=0;
193    choleskyStart_[iRow]=sizeFactor_;
194    if (!rowsDropped_[iRow]) {
195      CoinBigIndex startRow=rowStart[iRow];
196      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
197      for (CoinBigIndex k=startRow;k<endRow;k++) {
198        int iColumn=column[k];
199        CoinBigIndex start=columnStart[iColumn];
200        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
201        for (CoinBigIndex j=start;j<end;j++) {
202          int jRow=row[j];
203          if (jRow>=iRow&&!rowsDropped_[jRow]) {
204            if (!used[jRow]) {
205              used[jRow]=1;
206              which[number++]=jRow;
207            }
208          }
209        }
210      }
211      sizeFactor_ += number;
212      int j;
213      for (j=0;j<number;j++)
214        used[which[j]]=0;
215      // Sort
216      std::sort(which,which+number);
217      // move which on
218      which += number;
219    }
220  }
221  choleskyStart_[numberRows_]=sizeFactor_;
222  delete [] used;
223  permuteIn_ = new int [numberRows_];
224  permuteOut_ = new int[numberRows_];
225  integerParameters_[0]=0;
226  int i0=0;
227  int i1=1;
228  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
229         NULL,permuteOut_,permuteIn_,0,&numberRows_,&i1,
230         NULL,&i0,NULL,integerParameters_,doubleParameters_);
231  integerParameters_[1]=1;//order and symbolic
232  integerParameters_[2]=2;
233  integerParameters_[3]=0;//CSR
234  integerParameters_[4]=0;//C style
235  integerParameters_[13]=1;//reuse initial factorization space
236  integerParameters_[15+0]=1;//ordering
237  integerParameters_[15+1]=0;
238  integerParameters_[15+2]=1;
239  integerParameters_[15+3]=0;
240  integerParameters_[15+4]=1;
241  doubleParameters_[10]=1.0e-20;
242  doubleParameters_[11]=1.0e-15;
243  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
244         NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
245         NULL,&i0,NULL,integerParameters_,doubleParameters_);
246  //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
247  if (integerParameters_[63]) {
248    std::cout<<"wssmp returning error code of "<<integerParameters_[63]<<std::endl;
249    abort();
250  }
251  std::cout<<integerParameters_[23]<<" elements in sparse Cholesky"<<std::endl;
252  if (!integerParameters_[23]) {
253    std::cout<<"wssmp says no elements - fully dense? - switching to dense"<<std::endl;
254    integerParameters_[1]=2;
255    integerParameters_[2]=2;
256    integerParameters_[7]=1; // no permute
257    wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
258          NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
259          NULL,&i0,NULL,integerParameters_,doubleParameters_);
260    std::cout<<integerParameters_[23]<<" elements in dense Cholesky"<<std::endl;
261  }
262  return 0;
263}
264/* Factorize - filling in rowsDropped and returning number dropped */
265int 
266ClpCholeskyWssmp::factorize(const double * diagonal, int * rowsDropped) 
267{
268  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
269  const int * columnLength = model_->clpMatrix()->getVectorLengths();
270  const int * row = model_->clpMatrix()->getIndices();
271  const double * element = model_->clpMatrix()->getElements();
272  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
273  const int * rowLength = rowCopy_->getVectorLengths();
274  const int * column = rowCopy_->getIndices();
275  const double * elementByRow = rowCopy_->getElements();
276  int numberColumns=model_->clpMatrix()->getNumCols();
277  int iRow;
278  double * work = new double[numberRows_];
279  CoinZeroN(work,numberRows_);
280  const double * diagonalSlack = diagonal + numberColumns;
281  int newDropped=0;
282  double largest;
283  double smallest;
284  //perturbation
285  double perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
286  perturbation=perturbation*perturbation;
287  if (perturbation>1.0) {
288    //if (model_->model()->logLevel()&4)
289      std::cout <<"large perturbation "<<perturbation<<std::endl;
290    perturbation=sqrt(perturbation);;
291    perturbation=1.0;
292  } 
293  for (iRow=0;iRow<numberRows_;iRow++) {
294    double * put = sparseFactor_+choleskyStart_[iRow];
295    int * which = choleskyRow_+choleskyStart_[iRow];
296    int number = choleskyStart_[iRow+1]-choleskyStart_[iRow];
297    if (!rowsDropped_[iRow]) {
298      CoinBigIndex startRow=rowStart[iRow];
299      CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
300      work[iRow] = diagonalSlack[iRow];
301      for (CoinBigIndex k=startRow;k<endRow;k++) {
302        int iColumn=column[k];
303        CoinBigIndex start=columnStart[iColumn];
304        CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
305        double multiplier = diagonal[iColumn]*elementByRow[k];
306        for (CoinBigIndex j=start;j<end;j++) {
307          int jRow=row[j];
308          if (jRow>=iRow&&!rowsDropped_[jRow]) {
309            double value=element[j]*multiplier;
310            work[jRow] += value;
311          }
312        }
313      }
314      int j;
315      for (j=0;j<number;j++) {
316        int jRow = which[j];
317        put[j]=work[jRow];
318        work[jRow]=0.0;
319      }
320    } else {
321      // dropped
322      int j;
323      for (j=0;j<number;j++) {
324        put[j]=0.0;
325      }
326    }
327  }
328  //check sizes
329  double largest2=maximumAbsElement(sparseFactor_,sizeFactor_);
330  largest2*=1.0e-19;
331  largest = min (largest2,1.0e-11);
332  int numberDroppedBefore=0;
333  for (iRow=0;iRow<numberRows_;iRow++) {
334    int dropped=rowsDropped_[iRow];
335    // Move to int array
336    rowsDropped[iRow]=dropped;
337    if (!dropped) {
338      CoinBigIndex start = choleskyStart_[iRow];
339      double diagonal = sparseFactor_[start];
340      if (diagonal>largest2) {
341        sparseFactor_[start]=diagonal+perturbation;
342      } else {
343        sparseFactor_[start]=diagonal+perturbation;
344        rowsDropped[iRow]=2;
345        numberDroppedBefore++;
346      } 
347    } 
348  } 
349  int i1=1;
350  int i0=0;
351  integerParameters_[1]=3;
352  integerParameters_[2]=3;
353  integerParameters_[10]=2;
354  //integerParameters_[11]=1;
355  integerParameters_[12]=2;
356  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
357        NULL,permuteOut_,permuteIn_,NULL,&numberRows_,&i1,
358        NULL,&i0,rowsDropped,integerParameters_,doubleParameters_);
359  //    NULL,&i0,(int *) diagonal,integerParameters_,doubleParameters_);
360  //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
361  if (integerParameters_[9]) {
362    std::cout<<"scaling applied"<<std::endl;
363  } 
364  newDropped=integerParameters_[20]+numberDroppedBefore;
365  if (integerParameters_[20]) 
366    std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
367  largest=doubleParameters_[3];
368  smallest=doubleParameters_[4];
369  delete [] work;
370  if (model_->messageHandler()->logLevel()>1) 
371    std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
372  choleskyCondition_=largest/smallest;
373  bool cleanCholesky;
374  if (model_->numberIterations()<10) 
375    cleanCholesky=true;
376  else 
377    cleanCholesky=false;
378  if (cleanCholesky) {
379    //drop fresh makes some formADAT easier
380    int oldDropped=numberRowsDropped_;
381    if (newDropped||numberRowsDropped_) {
382      std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
383          newDropped<<" dropped)";
384      if (newDropped>oldDropped) 
385        std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
386      std::cout<<std::endl;
387      newDropped=0;
388      for (int i=0;i<numberRows_;i++) {
389        char dropped = rowsDropped[i];
390        rowsDropped_[i]=dropped;
391        if (dropped==2) {
392          //dropped this time
393          rowsDropped[newDropped++]=i;
394          rowsDropped_[i]=0;
395        } 
396      } 
397      numberRowsDropped_=newDropped;
398      newDropped=-(1+newDropped);
399    } 
400  } else {
401    if (newDropped) {
402      newDropped=0;
403      for (int i=0;i<numberRows_;i++) {
404        char dropped = rowsDropped[i];
405        rowsDropped_[i]=dropped;
406        if (dropped==2) {
407          //dropped this time
408          rowsDropped[newDropped++]=i;
409          rowsDropped_[i]=1;
410        } 
411      } 
412    } 
413    numberRowsDropped_+=newDropped;
414    if (numberRowsDropped_) {
415      std::cout <<"Rank "<<numberRows_-numberRowsDropped_<<" ( "<<
416          numberRowsDropped_<<" dropped)";
417      if (newDropped) {
418        std::cout<<" ( "<<newDropped<<" dropped this time)";
419      } 
420      std::cout<<std::endl;
421    } 
422  } 
423  status_=0;
424  return newDropped;
425}
426/* Uses factorization to solve. */
427void 
428ClpCholeskyWssmp::solve (double * region) 
429{
430  int i1=1;
431  int i0=0;
432  integerParameters_[1]=4;
433  integerParameters_[2]=4;
434#if 0
435  integerParameters_[5]=3;
436  doubleParameters_[5]=1.0e-10;
437  integerParameters_[6]=6;
438#endif
439  wssmp(&numberRows_,choleskyStart_,choleskyRow_,sparseFactor_,
440       NULL,permuteOut_,permuteIn_,region,&numberRows_,&i1,
441       NULL,&i0,NULL,integerParameters_,doubleParameters_);
442#if 0
443  if (integerParameters_[5]) {
444    std::cout<<integerParameters_[5]<<" refinements ";
445  }
446  std::cout<<doubleParameters_[6]<<std::endl;
447#endif
448}
Note: See TracBrowser for help on using the repository browser.