MOOSE - Multiscale Object Oriented Simulation Environment
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
standaloneTestFastElim.cpp
Go to the documentation of this file.
1 // To compile:
2 // g++ standaloneTestFastElim.cpp -lgsl -lgslcblas
3 
4 
5 #include <vector>
6 #include <algorithm>
7 #include <cassert>
8 #include <functional>
9 #include <iostream>
10 #include <iomanip>
11 // #include <gsl/gsl_linalg.h>
12 #include "/usr/include/gsl/gsl_linalg.h"
13 using namespace std;
14 #include "../basecode/SparseMatrix.h"
15 
16 const unsigned int SM_MAX_ROWS = 200000;
17 const unsigned int SM_MAX_COLUMNS = 200000;
18 const unsigned int SM_RESERVE = 8;
19 
20 void sortByColumn(
21  vector< unsigned int >& col, vector< double >& entry );
22 void testSorting();
23 
24 class Unroll
25 {
26  public:
27  Unroll( double diag, double off, unsigned int i, unsigned int j )
28  :
29  diagVal( diag ),
30  offDiagVal( off ),
31  row( i ),
32  col( j )
33  {;}
34  double diagVal;
35  double offDiagVal;
36  unsigned int row; // On which the diagonal is located
37  unsigned int col; // Col on which the offDiagVal is located.
38 };
39 
40 class FastElim: public SparseMatrix< double >
41 {
42  public:
43  void makeTestMatrix( const double* test, unsigned int numCompts );
44  /*
45  void rowElim( unsigned int row1, unsigned int row2,
46  vector< double >& rhs );
47  */
48  void buildForwardElim( vector< unsigned int >& diag,
49  vector< Triplet< double > >& fops );
50  void buildBackwardSub( vector< unsigned int >& diag,
51  vector< Triplet< double > >& bops, vector< double >& diagVal );
53  // Here we do stuff to set up the Hines ordering of the matrix.
55  bool hinesReorder( const vector< unsigned int >& parentVoxel );
56  const double* allEntries() const;
57  void shuffleRows(
58  const vector< unsigned int >& lookupOldRowFromNew );
59  /*
60  bool hinesReorder();
61  void extractTwig( unsigned int i,
62  vector< unsigned int >& rowReorder,
63  vector< bool >& extracted );
64  void findClosedEnds(
65  vector< unsigned int >& rowReorder,
66  vector< bool >& extracted );
67  void extractClosedEnds( unsigned int i,
68  vector< unsigned int >& rowReorder,
69  vector< bool >& extracted );
70  */
71 };
72 
73 const double* FastElim::allEntries() const {
74  return &N_[0];
75 }
76 
77 //
78 // static unsigned int parents[] = { 1,6,3,6,5,8,7,8,9,10,-1};
79 // unsigned int numKids[] = {0,1,0,1,0,2,
80 
85 bool FastElim::hinesReorder( const vector< unsigned int >& parentVoxel )
86 {
87  // First we fill in the vector that specifies the old row number
88  // assigned to each row of the reordered matrix.
89  assert( parentVoxel.size() == nrows_ );
90  vector< unsigned int > numKids( nrows_, 0 );
91  vector< unsigned int > lookupOldRowFromNew;
92  vector< bool > rowPending( nrows_, true );
93  unsigned int numDone = 0;
94  for ( unsigned int i = 0; i < nrows_; ++i ) {
95  if ( parentVoxel[i] != -1 )
96  numKids[ parentVoxel[i] ]++;
97  }
98  while ( numDone < nrows_ ) {
99  for ( unsigned int i = 0; i < nrows_; ++i ) {
100  if ( rowPending[i] && numKids[i] == 0 ) {
101  lookupOldRowFromNew.push_back( i );
102  rowPending[i] = false;
103  numDone++;
104  unsigned int pa = parentVoxel[i];
105  // Unsure what the root parent is. Assume it is -1
106  while ( pa != -1 && numKids[pa] == 1 ) {
107  assert( rowPending[pa] );
108  rowPending[pa] = false;
109  numDone++;
110  lookupOldRowFromNew.push_back( pa );
111  pa = parentVoxel[pa];
112  }
113  if ( pa != -1 ) {
114  assert( numKids[pa] > 0 );
115  numKids[pa]--;
116  }
117  }
118  }
119  }
120 
121  cout << setprecision(4);
122  cout << "oldRowFromNew= {" ;
123  for ( int i = 0; i < nrows_; ++i )
124  cout << lookupOldRowFromNew[i] << ", ";
125  cout << "}\n";
126  // Then we fill in the reordered matrix. Note we need to reorder
127  // columns too.
128  shuffleRows( lookupOldRowFromNew );
129 }
130 
131 // Fill in the reordered matrix. Note we need to reorder columns too.
132 void FastElim::shuffleRows(
133  const vector< unsigned int >& lookupOldRowFromNew )
134 {
135  vector< unsigned int > lookupNewRowFromOld( nrows_ );
136  for ( unsigned int i = 0; i < nrows_; ++i )
137  lookupNewRowFromOld[ lookupOldRowFromNew[i] ] = i;
138 
139  FastElim temp = *this;
140  clear();
141  setSize( temp.nrows_, temp.nrows_ );
142  for ( unsigned int i = 0; i < lookupOldRowFromNew.size(); ++i ) {
143  vector< unsigned int > c;
144  vector< double > e;
145  unsigned int num = temp.getRow( lookupOldRowFromNew[i], e, c );
146  vector< unsigned int > newc( num );
147  vector< double > newe( num );
148  for ( unsigned int j = 0; j < num; ++j ) {
149  newc[j] = lookupNewRowFromOld[ c[j] ];
150  newe[j] = e[j];
151  }
152  // Now we need to sort the new row entries in increasing col order.
153  /*
154  sortByColumn( newc, newe );
155  addRow( i, newe, newc );
156  */
157  sortByColumn( newc, e );
158  addRow( i, e, newc );
159  }
160 }
161 
162 void sortByColumn( vector< unsigned int >& col, vector< double >& entry )
163 {
164  unsigned int num = col.size();
165  assert( num == entry.size() );
166  // Stupid bubble sort, as we only have up to 5 entries and need to
167  // sort both the col and reorder the entries by the same sequence.
168  for ( unsigned int i = 0; i < num; ++i ) {
169  for ( unsigned int j = 1; j < num; ++j ) {
170  if ( col[j] < col[j-1] ) {
171  unsigned int temp = col[j];
172  col[j] = col[j-1];
173  col[j-1] = temp;
174  double v = entry[j];
175  entry[j] = entry[j-1];
176  entry[j-1] = v;
177  }
178  }
179  }
180 }
181 
182 /*
183 bool FastElim::hinesReorder()
184 {
185  vector< unsigned int > rowReorder;
186  vector< bool > extracted( nrows_, false );
187  for ( unsigned int i = 0; i < nrows_; ++i ) {
188  if ( rowStart_[i+1] - rowStart_[i] == 2 )
189  extractTwig( i, rowReorder, extracted );
190  }
191  // List of rows that now have all sub-branches extracted
192  vector< unsigned int > closedEnds;
193  do {
194  closedEnds.resize( 0 );
195  findClosedEnds( closedEnds, extracted );
196  for ( vector< unsigned int >::iterator
197  i = closedEnds.begin(); i != closedEnds.end(); ++i ) {
198  extractClosedEnds( *i, rowReorder, extracted );
199  }
200  } while ( closedEnds.size() > 0 );
201 }
202 */
203 
228 /*
229 void FastElim::rowElim( unsigned int row1, unsigned int row2,
230  vector< double >& rhs )
231 {
232  unsigned int rs1 = rowStart_[row1];
233  unsigned int rs2 = rowStart_[row2];
234  unsigned int n1 = rowStart_[row1+1] - rs1;
235  if ( n1 < 2 ) return;
236  assert( colIndex_[rs2] == row1 );
237  double diag1 = N_[diagIndex_[row1]];
238  double temp = N_[rs2];
239  double r = temp/diag1;
240 
241  const double* p1 = &(N_[ diagIndex_[row1] + 1 ] );
242  double* p2 = N_ +
243  for ( unsigned int i =
244  N_[rs2+1]
245 
246  double* v1 = N_ + rs1;
247  double* v2 = N_ + rs2;
248 
249 
250  v2' - v1'*v2/v1
251 
252 }
253 */
254 
255 void FastElim::makeTestMatrix( const double* test, unsigned int numCompts )
256 {
257  setSize( numCompts, numCompts );
258  vector< double > row( numCompts, ~0 );
259  unsigned int i = 1;
260  for ( unsigned int i = 0; i < numCompts; ++i ) {
261  for ( unsigned int j = 0; j < numCompts; ++j ) {
262  unsigned int k = i * numCompts + j;
263  if ( test[k] < 0.1 ) {
264  } else {
265  N_.push_back( test[k] );
266  colIndex_.push_back( j );
267  }
268  }
269  rowStart_[i + 1] = N_.size();
270  }
271 }
272 
273 /*
274 I need an outer function to fill the vector of ops for forward elim.
275 Then I need another outer function to fill another set of ops for
276 back-substitution.
277 */
278 
286 void FastElim::buildForwardElim( vector< unsigned int >& diag,
287  vector< Triplet< double > >& fops )
288 {
289  vector< vector< unsigned int > > rowsToElim( nrows_ );
290  diag.clear();
291  for ( unsigned int i = 0; i < nrows_; ++i ) {
292  unsigned int rs = rowStart_[i];
293  unsigned int re = rowStart_[i+1];
294  for ( unsigned int j = rs; j < re; ++j ) {
295  unsigned int k = colIndex_[j];
296  if ( k == i ) {
297  diag.push_back(j);
298  } else if ( k > i ) {
299  rowsToElim[ i ].push_back( k );
300  }
301  }
302  }
303  for ( unsigned int i = 0; i < nrows_; ++i ) {
304  double d = N_[diag[i]];
305  unsigned int diagend = rowStart_[ i + 1 ];
306  assert( diag[i] < diagend );
307  vector< unsigned int >& elim = rowsToElim[i];
308  for ( unsigned int j = 0; j < elim.size(); ++j ) {
309  unsigned int erow = elim[j];
310  if ( erow == i ) continue;
311  unsigned int rs = rowStart_[ erow ];
312  unsigned int re = rowStart_[ erow+1 ];
313  // assert( colIndex_[rs] == i );
314  double ratio = get( erow, i ) / d;
315  // double ratio = N_[rs]/N_[diag[i]];
316  for ( unsigned int k = diag[i]+1; k < diagend; ++k ) {
317  unsigned int col = colIndex_[k];
318  // findElimEntry, subtract it out.
319  for ( unsigned int q = rs; q < re; ++q ) {
320  if ( colIndex_[q] == col ) {
321  N_[q] -= N_[k] * ratio;
322  }
323  }
324  }
325  fops.push_back( Triplet< double >( ratio, i, erow) );
326  }
327  }
328  for ( unsigned int i = 0; i < rowsToElim.size(); ++i ) {
329  cout << i << " : ";
330  for ( unsigned int j = 0; j < rowsToElim[i].size(); ++j ) {
331  cout << rowsToElim[i][j] << " ";
332  }
333  cout << endl;
334  }
335  for ( unsigned int i = 0; i < fops.size(); ++i ) {
336  cout << "fops[" << i << "]= " << fops[i].b_ << " " << fops[i].c_ <<
337  " " << fops[i].a_ << endl;
338  }
339  /*
340  */
341 }
342 
352 void FastElim::buildBackwardSub( vector< unsigned int >& diag,
353  vector< Triplet< double > >& bops, vector< double >& diagVal )
354 {
355  // This vec tells the routine which rows below have to be back-subbed.
356  // This includes the rows if any in the tridiagonal band and also
357  // rows, if any, on branches.
358  vector< vector< unsigned int > > rowsToSub( nrows_ );
359 
360  for ( unsigned int i = 0; i < nrows_; ++i ) {
361  unsigned int d = diag[i] + 1;
362  unsigned int re = rowStart_[i+1];
363  for ( unsigned int j = d; j < re; ++j ) {
364  unsigned int k = colIndex_[j];
365  // At this point the row to sub is at (i, k). We need to go down
366  // to the (k,k) diagonal to sub it out.
367  rowsToSub[ k ].push_back( i );
368  }
369  }
370  for ( unsigned int i = 0; i < rowsToSub.size(); ++i ) {
371  cout << i << " : ";
372  for ( unsigned int j = 0; j < rowsToSub[i].size(); ++j ) {
373  cout << rowsToSub[i][j] << " ";
374  }
375  cout << endl;
376  }
377 
378  diagVal.resize( 0 );
379  // Fill in the diagonal terms. Here we do all entries.
380  for ( unsigned int i = 0; i != nrows_ ; ++i ) {
381  diagVal.push_back( 1.0 / N_[diag[i]] );
382  }
383 
384  // Fill in the back-sub operations. Note we don't need to check zero.
385  for ( unsigned int i = nrows_-1; i != 0 ; --i ) {
386  for ( unsigned int j = rowsToSub[i].size() - 1; j != -1; --j ) {
387  unsigned int k = rowsToSub[i][j];
388  double val = get( k, i ); //k is the row to go, i is the diag.
389  bops.push_back( Triplet< double >( val * diagVal[i], i, k ) );
390  }
391  }
392 
393  for ( unsigned int i = 0; i < bops.size(); ++i ) {
394  cout << i << ": " << bops[i].a_ << " " <<
395  bops[i].b_ << " " << // diagonal index
396  bops[i].c_ << " " << // off-diagonal index
397  1.0 / diagVal[bops[i].b_] << // diagonal value.
398  endl;
399  }
400 }
401 
402 void advance( vector< double >& y,
403  const vector< Triplet< double > >& ops, // has both fops and bops.
404  const vector< double >& diagVal )
405 {
406  for ( vector< Triplet< double > >::const_iterator
407  i = ops.begin(); i != ops.end(); ++i )
408  y[i->c_] -= y[i->b_] * i->a_;
409 
410  assert( y.size() == diagVal.size() );
411  vector< double >::iterator iy = y.begin();
412  for ( vector< double >::const_iterator
413  i = diagVal.begin(); i != diagVal.end(); ++i )
414  *iy++ *= *i;
415 }
416 
417 double checkAns(
418  const double* m, unsigned int numCompts,
419  const double* ans, const double* rhs )
420 {
421  vector< double > check( numCompts, 0.0 );
422  for ( unsigned int i = 0; i < numCompts; ++i ) {
423  for ( unsigned int j = 0; j < numCompts; ++j )
424  check[i] += m[i*numCompts + j] * ans[j];
425  }
426  double ret = 0.0;
427  for ( unsigned int i = 0; i < numCompts; ++i )
428  ret += (check[i] - rhs[i]) * (check[i] - rhs[i] );
429  return ret;
430 }
431 
432 
433 
434 main()
435 {
436 
437 
438 /*
439 2 11
440  1 4
441  3 10
442  9 5
443  6
444  7
445  8
446 
447  1 2 3 4 5 6 7 8 9 10 11
448 1 x x x x
449 2 x x
450 3 x x x x
451 4 x x x x
452 5 x x x x
453 6 x x x x
454 7 x x x
455 8 x x
456 9 x x x x
457 10 x x
458 11 x x
459  static double test[] = {
460  1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0,
461  5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0,
462  7, 0, 8, 9, 0, 0, 0, 0, 10, 0, 0,
463  11, 0, 12, 13, 0, 0, 0, 0, 0, 0, 14,
464  0, 0, 0, 0, 15, 16, 0, 0, 17, 18, 0,
465  0, 0, 0, 0, 19, 20, 21, 0, 22, 0, 0,
466  0, 0, 0, 0, 0, 23, 24, 25, 0, 0, 0,
467  0, 0, 0, 0, 0, 0, 26, 27, 0, 0, 0,
468  0, 0, 28, 0, 29, 30, 0, 0, 31, 0, 0,
469  0, 0, 0, 0, 32, 0, 0, 0, 0, 33, 0,
470  0, 0, 0, 34, 0, 0, 0, 0, 0, 0, 35,
471  };
472  const unsigned int numCompts = 11;
473 // static unsigned int parents[] = { 3,1,9,3,6,7,8,-1,6,5,4 };
474  static unsigned int parents[] = { 2,0,8,2,5,6,7,-1,5,4,3 };
475 */
476 
477 /*
478 1 3
479  2 4
480  7 5
481  8 6
482  9
483  10
484  11
485 
486  1 2 3 4 5 6 7 8 9 10 11
487 1 x x
488 2 x x x
489 3 x x
490 4 x x x
491 5 x x
492 6 x x x
493 7 x x x x
494 8 x x x
495 9 x x x x
496 10 x x x
497 11 x x
498  static double test[] = {
499  1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
500  3, 4, 0, 0, 0, 0, 5, 0, 0, 0, 0,
501  0, 0, 6, 7, 0, 0, 0, 0, 0, 0, 0,
502  0, 0, 8, 9, 0, 0, 10, 0, 0, 0, 0,
503  0, 0, 0, 0, 11, 12, 0, 0, 0, 0, 0,
504  0, 0, 0, 0, 13, 14, 0, 0, 15, 0, 0,
505  0, 16, 0, 17, 0, 0, 18, 19, 0, 0, 0,
506  0, 0, 0, 0, 0, 0, 20, 21, 22, 0, 0,
507  0, 0, 0, 0, 0, 23, 0, 24, 25, 26, 0,
508  0, 0, 0, 0, 0, 0, 0, 0, 27, 28, 29,
509  0, 0, 0, 0, 0, 0, 0, 0, 0, 30, 31,
510  };
511  const unsigned int numCompts = 11;
512  static unsigned int parents[] = { 1,6,3,6,5,8,7,8,9,10,-1};
513 */
514 
515 /*
516 Linear cable, 12 segments.
517 */
518 
519  static double test[] = {
520  1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
521  3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0,
522  0, 6, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0,
523  0, 0, 9, 10, 11, 0, 0, 0, 0, 0, 0, 0,
524  0, 0, 0, 12, 13, 14, 0, 0, 0, 0, 0, 0,
525  0, 0, 0, 0, 15, 16, 17, 0, 0, 0, 0, 0,
526  0, 0, 0, 0, 0, 18, 19, 20, 0, 0, 0, 0,
527  0, 0, 0, 0, 0, 0, 21, 22, 23, 0, 0, 0,
528  0, 0, 0, 0, 0, 0, 0, 24, 25, 26, 0, 0,
529  0, 0, 0, 0, 0, 0, 0, 0, 27, 28, 29, 0,
530  0, 0, 0, 0, 0, 0, 0, 0, 0, 30, 31, 32,
531  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 33, 34,
532  };
533  const unsigned int numCompts = 12;
534  static unsigned int parents[] = { 1,2,3,4,5,6,7,8,9,10,11,-1};
535 
536  /*
537  static double test[] = {
538  1, 2,
539  3, 4
540  };
541  const unsigned int numCompts = 2;
542  static double test[] = {
543  1, 2, 0, 0,
544  3, 4, 5, 0,
545  0, 6, 7, 8,
546  0, 0, 9, 10
547  };
548  const unsigned int numCompts = 4;
549  static double test[] = {
550  1, 2, 0, 0, 0, 0,
551  3, 4, 5, 0, 0, 0,
552  0, 6, 7, 8, 0, 0,
553  0, 0, 9, 10, 11, 0,
554  0, 0, 0, 12, 13, 14,
555  0, 0, 0, 0, 15, 16,
556  };
557  const unsigned int numCompts = 6;
558  static double test[] = {
559  1, 2, 0, 0, 0, 0,
560  3, 4, 0, 0, 1, 0,
561  0, 0, 7, 8, 0, 0,
562  0, 0, 9, 10, 11, 0,
563  0, 1, 0, 12, 13, 14,
564  0, 0, 0, 0, 15, 16,
565  };
566  const unsigned int numCompts = 6;
567  */
568  // testSorting(); // seems to work fine.
569  FastElim fe;
570  vector< Triplet< double > > fops;
571  vector< Unroll > bops;
572  fe.makeTestMatrix( test, numCompts );
573  fe.print();
574  cout << endl << endl;
575  vector< unsigned int > parentVoxel;
576  parentVoxel.insert( parentVoxel.begin(), &parents[0], &parents[numCompts] );
577  fe.hinesReorder( parentVoxel );
578  /*
579  */
580  /*
581  vector< unsigned int > shuf;
582  for ( unsigned int i = 0; i < numCompts; ++i )
583  shuf.push_back( i );
584  shuf[0] = 1;
585  shuf[1] = 0;
586  fe.shuffleRows( shuf );
587  */
588  fe.print();
589  cout << endl << endl;
590  FastElim foo = fe;
591 
592  vector< unsigned int > diag;
593  vector< double > diagVal;
594  fe.buildForwardElim( diag, fops );
595  fe.print();
596  fe.buildBackwardSub( diag, fops, diagVal );
597  vector< double > y( numCompts, 1.0 );
598  vector< double > ones( numCompts, 1.0 );
599  advance( y, fops, diagVal );
600  for ( int i = 0; i < numCompts; ++i )
601  cout << "y" << i << "]= " << y[i] << endl;
602 
603  // Here we verify the answer
604 
605  vector< double > alle;
606  for( unsigned int i = 0; i < numCompts; ++i ) {
607  for( unsigned int j = 0; j < numCompts; ++j ) {
608  alle.push_back( foo.get( i, j ) );
609  }
610  }
611  cout << "myCode: " <<
612  checkAns( &alle[0], numCompts, &y[0], &ones[0] ) << endl;
613 
614 
615 
617  // Here we do the gsl test.
618  vector< double > temp( &test[0], &test[numCompts*numCompts] );
619  gsl_matrix_view m = gsl_matrix_view_array( &temp[0], numCompts, numCompts );
620 
621  vector< double > z( numCompts, 1.0 );
622  gsl_vector_view b = gsl_vector_view_array( &z[0], numCompts );
623  gsl_vector* x = gsl_vector_alloc( numCompts );
624  int s;
625  gsl_permutation* p = gsl_permutation_alloc( numCompts );
626  gsl_linalg_LU_decomp( &m.matrix, p, &s );
627  gsl_linalg_LU_solve( &m.matrix, p, &b.vector, x);
628  vector< double > gslAns( numCompts );
629  for ( int i = 0; i < numCompts; ++i ) {
630  gslAns[i] = gsl_vector_get( x, i );
631  cout << "x[" << i << "]= " << gslAns[i] << endl;
632  }
633  /*
634  */
635  cout << "GSL: " << checkAns( test, numCompts, &gslAns[0], &ones[0] ) << endl;
636  gsl_vector_free( x );
637 
638 
639 }
640 
641 void testSorting()
642 {
643  static unsigned int k[] = {20,40,60,80,100,10,30,50,70,90};
644  static double d[] = {1,2,3,4,5,6,7,8,9,10};
645  vector< unsigned int > col;
646  col.insert( col.begin(), k, k+10);
647  vector< double > entry;
648  entry.insert( entry.begin(), d, d+10);
649  sortByColumn( col, entry );
650  cout << "testing sorting\n";
651  for ( int i = 0; i < col.size(); ++i ) {
652  cout << "d[" << i << "]= " << k[i] <<
653  ", col[" << i << "]= " << col[i] << ", e=" << entry[i] << endl;
654  }
655  cout << endl;
656 }
FastElim::allEntries
const double * allEntries() const
Definition: standaloneTestFastElim.cpp:73
Unroll::col
unsigned int col
Definition: standaloneTestFastElim.cpp:37
SparseMatrix::get
T get(unsigned int row, unsigned int column) const
Definition: SparseMatrix.h:253
FastElim::hinesReorder
bool hinesReorder(const vector< unsigned int > &parentVoxel)
Definition: standaloneTestFastElim.cpp:85
SparseMatrix::getRow
unsigned int getRow(unsigned int row, const T **entry, const unsigned int **colIndex) const
Definition: SparseMatrix.h:288
FastElim::buildBackwardSub
void buildBackwardSub(vector< unsigned int > &diag, vector< Triplet< double > > &bops, vector< double > &diagVal)
Definition: standaloneTestFastElim.cpp:352
FastElim::shuffleRows
void shuffleRows(const vector< unsigned int > &lookupOldRowFromNew)
Definition: standaloneTestFastElim.cpp:132
advance
void advance(vector< double > &y, const vector< Triplet< double > > &ops, const vector< double > &diagVal)
Definition: standaloneTestFastElim.cpp:402
testSorting
void testSorting()
Definition: standaloneTestFastElim.cpp:641
Unroll::Unroll
Unroll(double diag, double off, unsigned int i, unsigned int j)
Definition: standaloneTestFastElim.cpp:27
SM_MAX_COLUMNS
const unsigned int SM_MAX_COLUMNS
Definition: standaloneTestFastElim.cpp:17
FastElim::makeTestMatrix
void makeTestMatrix(const double *test, unsigned int numCompts)
Definition: standaloneTestFastElim.cpp:255
Unroll::row
unsigned int row
Definition: standaloneTestFastElim.cpp:36
SparseMatrix::nrows_
unsigned int nrows_
Definition: SparseMatrix.h:701
sortByColumn
void sortByColumn(vector< unsigned int > &col, vector< double > &entry)
Definition: standaloneTestFastElim.cpp:162
checkAns
double checkAns(const double *m, unsigned int numCompts, const double *ans, const double *rhs)
Definition: standaloneTestFastElim.cpp:417
FastElim::buildForwardElim
void buildForwardElim(vector< unsigned int > &diag, vector< Triplet< double > > &fops)
Definition: standaloneTestFastElim.cpp:286
SparseMatrix::print
void print() const
Definition: SparseMatrix.h:650
SM_RESERVE
const unsigned int SM_RESERVE
Definition: standaloneTestFastElim.cpp:18
SM_MAX_ROWS
const unsigned int SM_MAX_ROWS
Definition: standaloneTestFastElim.cpp:16