Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
basker_def.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Basker: A Direct Linear Solver package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, with Sandia Corporation, the
8 // U.S. Government retains certain rights in this software.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
23 // USA
24 // Questions? Contact Mike A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef BASKER_DEF_HPP
30 #define BASKER_DEF_HPP
31 
32 #include "basker_decl.hpp"
33 #include "basker_scalartraits.hpp"
34 //#include "basker.hpp"
35 
36 //#include <cassert>
37 #include <iostream>
38 #include <cstdio>
39 
40 using namespace std;
41 
42 //#define BASKER_DEBUG 1
43 //#undef UDEBUG
44 
45 namespace BaskerClassicNS{
46 
47  template <class Int, class Entry>
48  BaskerClassic<Int, Entry>::BaskerClassic()
49  {
50 
51  //A = (basker_matrix<Int,Entry> *) malloc(sizeof(basker_matrix<Int,Entry>));
52  A = new basker_matrix<Int,Entry>;
53 
54  //L = (basker_matrix<Int,Entry> *) malloc(sizeof(basker_matrix<Int,Entry>));
55  L = new basker_matrix<Int, Entry>;
56  L->nnz = 0;
57 
58  //U = (basker_matrix<Int,Entry> *) malloc(sizeof(basker_matrix<Int,Entry>));
59  U = new basker_matrix<Int,Entry>;
60  U->nnz = 0;
61 
62  actual_lnnz = Int(0);
63  actual_unnz = Int(0);
64 
65  been_fact = false;
66  perm_flag = false;
67  }
68 
69 
70  template <class Int, class Entry>
71  BaskerClassic<Int, Entry>::BaskerClassic(Int nnzL, Int nnzU)
72  {
73 
74  //A = (basker_matrix<Int, Entry> *) malloc(sizeof(basker_matrix<Int,Entry>));
75  A = new basker_matrix<Int, Entry>;
76  //L = (basker_matrix<Int, Entry> *) malloc(sizeof(basker_matrix<Int,Entry>));
77  L = new basker_matrix<Int, Entry>;
78  L->nnz = nnzL;
79  //U = (basker_matrix<Int, Entry> *) malloc(sizeof(basker_matrix<Int,Entry>));
80  U = new basker_matrix<Int, Entry>;
81  U->nnz = nnzU;
82 
83  actual_lnnz = Int(0);
84  actual_unnz = Int(0);
85 
86  been_fact = false;
87  perm_flag = false;
88  }
89 
90 
91  template <class Int, class Entry>
92  BaskerClassic<Int, Entry>::~BaskerClassic()
93  {
94  //free factor
95  if(been_fact)
96  {
97  free_factor();
98  //BASKERFREE(pinv);
99  delete pinv;
100  }
101  if(perm_flag)
102  {
103  //free_perm_matrix();
104  }
105  //BASKERFREE(A);
106  delete A;
107  //BASKERFREE(L);
108  delete L;
109  //BASKERFREE(U);
110  delete U;
111  }
112 
113 
114  template <class Int, class Entry>
115  int BaskerClassic<Int,Entry>:: basker_dfs
116  (
117  Int n,
118  Int j,
119  Int *Li,
120  Int *Lp,
121  Int *color,
122  Int *pattern, /* o/p */
123  Int *top, /* o/p */
124  //Int k,
125  Int *tpinv,
126  Int *stack
127  )
128  {
129 
130  Int i, t, i1, head ;
131  Int start, end, done, *store ;
132 
133  store = stack + n ;
134  head = 0;
135  stack[head] = j;
136  bool has_elements = true;
137 
138  while(has_elements)
139  {
140  j = stack[head] ;
141 #ifdef BASKER_DEBUG
142  //std::cout << "DFS: " << j << "COLOR: " << color[j] << std::endl;
143 #endif
144  t = tpinv [j] ;
145  if (color[j] == 0)
146  {
147  /* Seeing this column for first time */
148  color[j] = 1 ;
149  start = Lp[t] ;
150  }
151  else
152  {
153  BASKERASSERT (color[j] == 1) ; /* color cannot be 2 when we are here */
154  start = store[j];
155  }
156  done = 1;
157 
158  if ( t != n )
159  {
160  end = Lp[t+1] ;
161  for ( i1 = start ; i1 < end ; i1++ )
162  {
163  i = Li[i1] ;
164  if ( color[i] == 0 )
165  {
166  stack[++head] = i;
167  store[j] = i1+1;
168  done = 0;
169  break;
170  }
171  }
172  }
173  if (done)
174  {
175  pattern[--*top] = j ;
176  color[j] = 2 ;
177  if(head == 0)
178  {
179  has_elements = false;
180  }
181  else
182  {
183  head--;
184  }
185  }
186  }
187 #ifdef BASKER_DEBUG
188  std::cout << "Out of DFS: " << j << std::endl;
189 #endif
190  return 0;
191  } //End dfs
192 
193  template <class Int, class Entry>
194  int BaskerClassic<Int,Entry>::factor(Int nrow, Int ncol , Int nnz, Int *col_ptr, Int *row_idx, Entry *val)
195  {
196  int ierr = 0;
197  /*Initalize A basker matrix struc */
198 #ifdef BASKER_DEBUG
199 
200  BASKERASSERT(nrow > 0);
201  BASKERASSERT(ncol > 0);
202  BASKERASSERT(nnz > 0);
203 
204 #endif
205 
206  A->nrow = nrow;
207  A->ncol = ncol;
208  A->nnz = nnz;
209  A->col_ptr = col_ptr;
210  A->row_idx = row_idx;
211  A->val = val;
212  /*End initalize A*/
213 
214  /*Creating space for L and U*/
215  L->nrow = nrow;
216  L->ncol = ncol;
217  if(L->nnz == 0)
218  {
219  L->nnz = 2*A->nnz;
220  }
221  //L->col_ptr = (Int *) BASKERCALLOC(ncol+1, sizeof(Int));
222  L->col_ptr = new Int[ncol+1]();
223  //L->row_idx = (Int *) BASKERCALLOC(L->nnz, sizeof(Int));
224  L->row_idx = new Int[L->nnz]();
225  //L->val = (Entry *) BASKERCALLOC(L->nnz, sizeof(Entry));
226  L->val = new Entry[L->nnz]();
227 
228  U->nrow = nrow;
229  U->ncol = ncol;
230  if(U->nnz == 0)
231  {
232  U->nnz = 2*A->nnz;
233  }
234  //U->col_ptr = (Int *) BASKERCALLOC(ncol+1, sizeof(Int));
235  U->col_ptr = new Int[ncol+1]();
236  //U->row_idx = (Int *) BASKERCALLOC(U->nnz, sizeof(Int));
237  U->row_idx = new Int[U->nnz]();
238  //U->val = (Entry *) BASKERCALLOC(U->nnz, sizeof(Entry));
239  U->val = new Entry[U->nnz]();
240 
241  if((L->col_ptr == NULL) || (L->row_idx == NULL) || (L->val == NULL) ||
242  (U->col_ptr == NULL) || (U->row_idx == NULL) || (U->val == NULL))
243  {
244  ierr = -1;
245  return ierr;
246  }
247  /*End creating space for L and U*/
248 
249  /*Creating working space*/
250  Int *tptr;
251  Entry *X;
252  //tptr = (Int *) BASKERCALLOC( (ncol)+(4*nrow), sizeof(Int));
253  tptr = new Int[(ncol)+(4*nrow)]();
254  //X = (Entry *) BASKERCALLOC(2*nrow, sizeof(Entry));
255  X = new Entry[2*nrow]();
256  //pinv = (Int * ) BASKERCALLOC(ncol+1, sizeof(Int)); //Note extra pad
257  pinv = new Int[ncol+1]();
258 
259 
260  if( (tptr == NULL) || (X == NULL) || (pinv == NULL) )
261  {
262  ierr = -2;
263  return ierr;
264  }
265 
266  /*End creating working space */
267 
268  /*Defining Variables Used*/
269  Int i, j, k;
270  Int *color, *pattern, *stack; // pointers into the work space
271  Int top, top1, maxindex, t; // j1, j2;
272  Int lnnz, unnz, xnnz, lcnt, ucnt;
273  Int cu_ltop, cu_utop;
274  Int pp, p2, p;
275  Int newsize;
276  Entry pivot, value, xj;
277  Entry absv, maxv;
278 
279  color = tptr;
280  tptr += ncol;
281 
282  pattern = tptr;
283  tptr += nrow;
284 
285  stack = tptr;
286  tptr += 2*(nrow);
287 
288  cu_ltop = 0;
289  cu_utop = 0;
290  top = ncol;
291  top1 = ncol;
292  lnnz = 0; //real found lnnz
293  unnz = 0; //real found unnz
294 
295  for(k = 0 ; k < ncol; k++)
296  {
297  pinv[k] = ncol;
298  }
299 
300  /*For all columns in A .... */
301  for (k = 0; k < ncol; k++)
302  {
303 
304 #ifdef BASKER_DEBUG
305  cout << "k = " << k << endl;
306 #endif
307 
308  value = 0.0;
309  pivot = 0.0;
310  maxindex = ncol;
311  //j1 = 0;
312  //j2 = 0;
313  lcnt = 0;
314  ucnt = 0;
315 
316 #ifdef BASKER_DEBUG
317  BASKERASSERT (top == ncol);
318 
319  for(i = 0; i < nrow; i++)
320  {
321  BASKERASSERT(X[i] == (Entry)0);
322  }
323  for(i = 0; i < ncol; i++)
324  {
325  BASKERASSERT(color[i] == 0);
326  }
327 #endif
328  /* Reachability for every nonzero in Ak */
329  for( i = col_ptr[k]; i < col_ptr[k+1]; i++)
330  {
331  j = row_idx[i];
332  X[j] = val[i];
333 
334  if(color[j] == 0)
335  {
336  //do dfs
337  basker_dfs(nrow, j, L->row_idx, L->col_ptr, color, pattern, &top, pinv, stack);
338 
339  }
340 
341  }//end reachable
342 
343  xnnz = ncol - top;
344 #ifdef BASKER_DEBUG
345  cout << top << endl;
346  cout << ncol << endl;
347  cout << xnnz << endl;
348  //BASKERASSERT(xnnz <= nrow);
349 #endif
350  /*Lx = b where x will be the column k in L and U*/
351  top1 = top;
352  for(pp = 0; pp < xnnz; pp++)
353  {
354  j = pattern[top1++];
355  color[j] = 0;
356  t = pinv[j];
357 
358  if(t!=ncol) //it has not been assigned
359  {
360  xj = X[j];
361  p2 = L->col_ptr[t+1];
362  for(p = L->col_ptr[t]+1; p < p2; p++)
363  {
364  X[L->row_idx[p]] -= L->val[p] * xj;
365  }//over all rows
366  }
367 
368  }
369 
370  /*get the pivot*/
371  maxv = 0.0;
372  for(i = top; i < nrow; i++)
373  {
374  j = pattern[i];
375  t = pinv[j];
376  value = X[j];
377  /*note may want to change this to traits*/
378  //absv = (value < 0.0 ? -value : value);
379  absv = BASKER_ScalarTraits<Entry>::approxABS(value);
380 
381  if(t == ncol)
382  {
383  lcnt++;
384  if( BASKER_ScalarTraits<Entry>::gt(absv , maxv))
385  //if(absv > BASKER_ScalarTraits<Entry>::approxABS(maxv))
386  {
387  maxv = absv;
388  pivot = value;
389  maxindex= j;
390  }
391  }
392  }
393  ucnt = nrow - top - lcnt + 1;
394 
395  if(maxindex == ncol || pivot == ((Entry)0))
396  {
397  cout << "Matrix is singular at index: " << maxindex << " pivot: " << pivot << endl;
398  ierr = maxindex;
399  return ierr;
400  }
401 
402  pinv[maxindex] = k;
403 #ifdef BASKER_DEBUG
404  if(maxindex != k )
405  {
406  cout << "Permuting pivot: " << k << " for row: " << maxindex << endl;
407  }
408 #endif
409 
410  if(lnnz + lcnt >= L->nnz)
411  {
412 
413  newsize = L->nnz * 1.1 + 2*nrow + 1;
414 #ifdef BASKER_DEBUG
415  cout << "Out of memory -- Reallocating. Old Size: " << L->nnz << " New Size: " << newsize << endl;
416 #endif
417  //L->row_idx = (Int *) BASKERREALLOC(L->row_idx, newsize*sizeof(Int));
418  L->row_idx = int_realloc(L->row_idx , L->nnz, newsize);
419  if(!(L->row_idx))
420  {
421  cout << "WARNING: Cannot Realloc Memory" << endl;
422  ierr = -3;
423  return ierr;
424  }
425  //L->val = (Entry *) BASKERREALLOC(L->val, newsize*sizeof(Entry));
426  L->val = entry_realloc(L->val, L->nnz, newsize);
427  if(!(L->val))
428  {
429  cout << "WARNING: Cannot Realloc Memory" << endl;
430  ierr = -3;
431  return ierr;
432  }
433  L->nnz = newsize;
434 
435  }//realloc if L is out of memory
436 
437  if(unnz + ucnt >= U->nnz)
438  {
439  newsize = U->nnz*1.1 + 2*nrow + 1;
440 #ifdef BASKER_DEBUG
441  cout << "Out of memory -- Reallocating. Old Size: " << U->nnz << " New Size: " << newsize << endl;
442 #endif
443  //U->row_idx = (Int *) BASKERREALLOC(U->row_idx, newsize*sizeof(Int));
444  U->row_idx = int_realloc(U->row_idx, U->nnz, newsize);
445  if(!(U->row_idx))
446  {
447  cout << "WARNING: Cannot Realloc Memory" << endl;
448  ierr = -3;
449  return ierr;
450  }
451 
452  //U->val = (Entry *) BASKERREALLOC(U->val, newsize*sizeof(Entry));
453  U->val = entry_realloc(U->val, U->nnz, newsize);
454  if(!(U->val))
455  {
456  cout << "WARNING: Cannot Realloc Memory" << endl;
457  ierr = -3;
458  return ierr;
459  }
460  U->nnz = newsize;
461  }//realloc if U is out of memory
462 
463  //L->col_ptr[lnnz] = maxindex;
464  L->row_idx[lnnz] = maxindex;
465  L->val[lnnz] = 1.0;
466  lnnz++;
467 
468  Entry last_v_temp = 0;
469 
470  for(i = top; i < nrow; i++)
471  {
472  j = pattern[i];
473  t = pinv[j];
474 
475  /* check for numerical cancellations */
476 
477 
478  if(X[j] != ((Entry)0))
479  {
480 
481  if(t != ncol)
482  {
483  if(unnz >= U->nnz)
484  {
485  cout << "BASKER: Insufficent memory for U" << endl;
486  ierr = -3;
487  return ierr;
488  }
489  if(t < k)
490  //if(true)
491  {
492  U->row_idx[unnz] = pinv[j];
493  U->val[unnz] = X[j];
494  unnz++;
495  }
496  else
497  {
498 
499  last_v_temp = X[j];
500  //cout << "Called. t: " << t << "Val: " << last_v_temp << endl;
501  }
502 
503  }
504  else if (t == ncol)
505  {
506  if(lnnz >= L->nnz)
507  {
508  cout << "BASKER: Insufficent memory for L" << endl;
509  ierr = -3;
510  return ierr;
511  }
512 
513  L->row_idx[lnnz] = j;
514  //L->val[lnnz] = X[j]/pivot;
515  L->val[lnnz] = BASKER_ScalarTraits<Entry>::divide(X[j],pivot);
516  lnnz++;
517 
518  }
519 
520  }
521 
522 
523  X[j] = 0;
524 
525  }
526  //cout << "Value added at end: " << last_v_temp << endl;
527  U->row_idx[unnz] = k;
528  U->val[unnz] = last_v_temp;
529  unnz++;
530 
531  xnnz = 0;
532  top = ncol;
533 
534  L->col_ptr[k] = cu_ltop;
535  L->col_ptr[k+1] = lnnz;
536  cu_ltop = lnnz;
537 
538  U->col_ptr[k] = cu_utop;
539  U->col_ptr[k+1] = unnz;
540  cu_utop = unnz;
541 
542  } //end for every column
543 
544 #ifdef BASKER_DEBUG
545  /*Print out found L and U*/
546  for(k = 0; k < lnnz; k++)
547  {
548  printf("L[%d]=%g" , k , L->val[k]);
549  }
550  cout << endl;
551  for(k = 0; k < lnnz; k++)
552  {
553  printf("Li[%d]=%d", k, L->row_idx[k]);
554  }
555  cout << endl;
556  for(k = 0; k < nrow; k++)
557  {
558  printf("p[%d]=%d", k, pinv[k]);
559  }
560  cout << endl;
561  cout << endl;
562 
563  for(k = 0; k < ncol; k++)
564  {
565  printf("Up[%d]=%d", k, U->col_ptr[k]);
566  }
567  cout << endl;
568 
569  for(k = 0; k < unnz; k++)
570  {
571  printf("U[%d]=%g" , k , U->val[k]);
572  }
573  cout << endl;
574  for(k = 0; k < unnz; k++)
575  {
576  printf("Ui[%d]=%d", k, U->row_idx[k]);
577  }
578  cout << endl;
579 
580 
581 #endif
582  /* Repermute */
583  for( i = 0; i < ncol; i++)
584  {
585  for(k = L->col_ptr[i]; k < L->col_ptr[i+1]; k++)
586  {
587  //L->row_idx[k] = pinv[L->row_idx[k]];
588  }
589  }
590  //Max sure correct location of min in L and max in U for CSC format//
591  //Speeds up tri-solve//
592  //sort_factors();
593 
594 #ifdef BASKER_DEBUG
595  cout << "After Permuting" << endl;
596  for(k = 0; k < lnnz; k++)
597  {
598  printf("Li[%d]=%d", k, L->row_idx[k]);
599  }
600  cout << endl;
601 #endif
602 
603  //BASKERFREE(X);
604  //BASKERFREE(tptr);
605 
606  actual_lnnz = lnnz;
607  actual_unnz = unnz;
608 
609  been_fact = true;
610  return 0;
611  }//end factor
612 
613 
614  template <class Int, class Entry>
615  Int BaskerClassic<Int, Entry>::get_NnzL()
616  {
617  return actual_lnnz;
618  }
619 
620  template <class Int, class Entry>
621  Int BaskerClassic<Int, Entry>::get_NnzU()
622  {
623  return actual_unnz;
624  }
625 
626  template <class Int, class Entry>
627  Int BaskerClassic<Int, Entry>::get_NnzLU()
628  {
629  return (actual_lnnz + actual_unnz);
630  }
631 
632  template <class Int, class Entry>
633  int BaskerClassic<Int, Entry>::returnL(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
634  {
635  int i;
636  *dim = L->nrow;
637  *nnz = L->nnz;
638 
639  /*Does a bad copy*/
640 
641  //*col_ptr = (Int *) BASKERCALLOC(L->nrow+1, sizeof(Int));
642  *col_ptr = new Int[L->nrow+1];
643  //*row_idx = (Int *) BASKERCALLOC(L->nnz, sizeof(Int));
644  *row_idx = new Int[L->nnz];
645  //*val = (Entry *) BASKERCALLOC(L->nnz, sizeof(Entry));
646  *val = new Entry[L->nnz];
647 
648  if( (*col_ptr == NULL) || (*row_idx == NULL) || (*val == NULL) )
649  {
650  return -1;
651  }
652 
653  for(i = 0; i < L->nrow+1; i++)
654  {
655  (*col_ptr)[i] = L->col_ptr[i];
656  }
657 
658  for(i = 0; i < actual_lnnz; i++)
659  {
660  (*row_idx)[i] = pinv[L->row_idx[i]];
661  (*val)[i] = L->val[i];
662  }
663  return 0;
664 
665  }
666 
667  template <class Int, class Entry>
668  int BaskerClassic<Int, Entry>::returnU(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
669  {
670  int i;
671  *dim = U->nrow;
672  *nnz = U->nnz;
673  /*Does a bad copy*/
674  //*col_ptr = (Int *) BASKERCALLOC(U->nrow+1, sizeof(Int));
675  *col_ptr = new Int[U->nrow+1];
676  //*row_idx = (Int *) BASKERCALLOC(U->nnz, sizeof(Int));
677  *row_idx = new Int[U->nnz];
678  //*val = (Entry *) BASKERCALLOC(U->nnz, sizeof(Entry));
679  *val = new Entry[U->nnz];
680 
681  if( (*col_ptr == NULL) || (*row_idx == NULL) || (*val == NULL) )
682  {
683  return -1;
684  }
685 
686  for(i = 0; i < U->nrow+1; i++)
687  {
688  (*col_ptr)[i] = U->col_ptr[i];
689  }
690  for(i = 0; i < actual_unnz; i++)
691  {
692  (*row_idx)[i] = U->row_idx[i];
693  (*val)[i] = U->val[i];
694  }
695  return 0;
696  }
697 
698  template <class Int, class Entry>
699  int BaskerClassic<Int, Entry>::returnP(Int** p)
700  {
701  Int i;
702  //*p = (Int *) BASKERCALLOC(A->nrow, sizeof(Int));
703  *p = new Int[A->nrow];
704 
705  if( (*p == NULL ) )
706  {
707  return -1;
708  }
709 
710  for(i = 0; i < A->nrow; i++)
711  {
712  (*p)[pinv[i]] = i; //Matlab perm-style
713  }
714  return 0;
715  }
716 
717  template <class Int, class Entry>
718  void BaskerClassic<Int, Entry>::free_factor()
719  {
720  //BASKERFREE L
721  //BASKERFREE(L->col_ptr);
722  delete[] L->col_ptr;
723  //BASKERFREE(L->row_idx);
724  delete[] L->row_idx;
725  //BASKERFREE(L->val);
726  delete[] L->val;
727 
728  //BASKERFREE U
729  //BASKERFREE(U->col_ptr);
730  delete[] U->col_ptr;
731  //BASKERFREE(U->row_idx);
732  delete[] U->row_idx;
733  //BASKERFREE(U->val);
734  delete[] U->val;
735 
736  }
737  template <class Int, class Entry>
738  void BaskerClassic<Int, Entry>::free_perm_matrix()
739  {
740  //BASKERFREE(A->col_ptr);
741  //BASKERFREE(A->row_idx);
742  //BASKERFREE(A->val);
743  }
744 
745  template <class Int, class Entry>
746  int BaskerClassic<Int, Entry>::solveMultiple(Int nrhs, Entry *b, Entry *x)
747  {
748  Int i;
749  for(i = 0; i < nrhs; i++)
750  {
751  Int k = i*A->nrow;
752  int result = solve(&(b[k]), &(x[k]));
753  if(result != 0)
754  {
755  cout << "Error in Solving \n";
756  return result;
757  }
758  }
759  return 0;
760  }
761 
762 
763  template <class Int, class Entry>
764  int BaskerClassic<Int, Entry>::solve(Entry* b, Entry* x)
765  {
766 
767  if(!been_fact)
768  {
769  return -10;
770  }
771  //Entry* temp = (Entry *)BASKERCALLOC(A->nrow, sizeof(Entry));
772  Entry* temp = new Entry[A->nrow]();
773  Int i;
774  int result = 0;
775  for(i = 0 ; i < A->ncol; i++)
776  {
777  Int k = pinv[i];
778  x[k] = b[i];
779  }
780 
781  result = low_tri_solve_csc(L->nrow, L->col_ptr, L->row_idx, L->val, temp, x);
782  if(result == 0)
783  {
784  result = up_tri_solve_csc(U->nrow, U->col_ptr, U->row_idx, U->val, x, temp);
785  }
786 
787 
788  //BASKERFREE(temp);
789  delete[] temp;
790  return 0;
791  }
792 
793  template < class Int, class Entry>
794  int BaskerClassic<Int, Entry>::low_tri_solve_csc( Int n, Int *col_ptr, Int *row_idx, Entry* val, Entry *x, Entry *b)
795  {
796  Int i, j;
797  /*for each column*/
798  for(i = 0; i < n ; i++)
799  {
800 #ifdef BASKER_DEBUG
801  BASKERASSERT(val[col_ptr[i]] != (Entry)0);
802 #else
803  if(val[col_ptr[i]] == (Entry) 0)
804  {
805  return i;
806  }
807 #endif
808  x[i] = BASKER_ScalarTraits<Entry>::divide(b[i], val[col_ptr[i]]);
809 
810  for(j = col_ptr[i]+1; j < (col_ptr[i+1]); j++) //update all rows
811  {
812  b[pinv[row_idx[j]]] = b[pinv[row_idx[j]]] - (val[j]*x[i]);
813  }
814  }
815  return 0;
816  }
817 
818  template < class Int, class Entry>
819  int BaskerClassic<Int, Entry>::up_tri_solve_csc( Int n, Int *col_ptr, Int *row_idx, Entry *val, Entry *x, Entry *b)
820  {
821  Int i, j;
822  /*for each column*/
823  for(i = n; i > 1 ; i--)
824  {
825  int ii = i-1;
826 #ifdef BASKER_DEBUG
827  BASKERASSERT(val[col_ptr[i]-1] != (Entry)0);
828 #else
829  if(val[col_ptr[i]-1] == (Entry) 0)
830  {
831  cout << "Dig(" << i << ") = " << val[col_ptr[i]-1] << endl;
832  return i;
833  }
834 #endif
835  //x[ii] = b[ii]/val[col_ptr[i]-1]; //diag
836  x[ii] = BASKER_ScalarTraits<Entry>::divide(b[ii],val[col_ptr[i]-1]);
837 
838  for(j = (col_ptr[i]-2); j >= (col_ptr[ii]); j--)
839  {
840  b[row_idx[j]] = b[row_idx[j]] - (val[j]*x[ii]);
841  }
842  }
843  //x[0] = b[0]/val[col_ptr[1]-1];
844  x[0] = BASKER_ScalarTraits<Entry>::divide(b[0],val[col_ptr[1]-1]);
845  return 0;
846  }
847 
848  template <class Int, class Entry>
849  int BaskerClassic<Int, Entry>::preorder(Int *row_perm, Int *col_perm)
850  {
851 
852  basker_matrix <Int, Entry> *B;
853  B = new basker_matrix<Int, Entry>;
854  B->nrow = A->nrow;
855  B->ncol = A->ncol;
856  B->nnz = A->nnz;
857  B->col_ptr = (Int *) BASKERCALLOC(A->ncol + 1, sizeof(Int));
858  B->row_idx = (Int *) BASKERCALLOC(A->nnz, sizeof(Int));
859  B->val = (Entry *) BASKERCALLOC(A->val, sizeof(Int));
860 
861  if( (B->col_ptr == NULL) || (B->row_idx == NULL) || (B->val == NULL) )
862  {
863  perm_flag = false;
864  return -1;
865  }
866 
867  /* int resultcol = (unused) */ (void) permute_column(col_perm, B);
868  /* int resultrow = (unused) */ (void) permute_row(row_perm, B);
869 
870  /*Note: the csc matrices of A are the problem of the user
871  therefore we will not free them*/
872  A->col_ptr = B->col_ptr;
873  A->row_idx = B->row_idx;
874  A->val = A->val;
875 
876  perm_flag = true; /*Now we will free A at the end*/
877 
878  return 0;
879  }
880 
881  template <class Int, class Entry>
882  int BaskerClassic <Int, Entry>::permute_column(Int *p, basker_matrix<Int,Entry> *B)
883  {
884  /*p(i) contains the destination of row i in the permuted matrix*/
885  Int i,j, ii, jj;
886 
887  /*Determine column pointer of output matrix*/
888  for(j=0; j < B->ncol; j++)
889  {
890  i = p[j];
891  B->col_ptr[i+1] = A->col_ptr[j+1] - A->col_ptr[j];
892  }
893  /*get pointers from lengths*/
894  B->col_ptr[0] = 0;
895  for(j=0; j < B->ncol; j++)
896  {
897  B->col_ptr[j+1] = B->col_ptr[j+1] + B->col_ptr[j];
898  }
899 
900  /*copy idxs*/
901  Int k, ko;
902  for(ii = 0 ; ii < B->ncol; ii++)
903  {// old colum ii new column p[ii] k->pointer
904  ko = B->col_ptr(p[ii]);
905  for(k = A->col_ptr[ii]; k < A->col_ptr[ii+1]; k++)
906  {
907  B->row_index[ko] = A->row_index[k];
908  B->val[ko] = A->val[ko];
909  ko++;
910  }
911  }
912  return 0;
913  }
914 
915  template <class Int, class Entry>
916  int BaskerClassic <Int, Entry>::permute_row(Int *p, basker_matrix<Int,Entry> *B)
917  {
918  Int k,i;
919  for(k=0; k < A->nnz; k++)
920  {
921  B->row_idx[k] = p[A->row_idx[k]];
922  }
923  return 0;
924  }
925 
926  template <class Int, class Entry>
927  int BaskerClassic <Int, Entry>::sort_factors()
928  {
929 
930  /*Sort CSC of L - just make sure min_index is in lowest position*/
931  Int i, j;
932  Int p;
933  Int val;
934  for(i = 0 ; i < L->ncol; i++)
935  {
936  p = L->col_ptr[i];
937  val = L->row_idx[p];
938 
939  for(j = L->col_ptr[i]+1; j < (L->col_ptr[i+1]); j++)
940  {
941  if(L->row_idx[j] < val)
942  {
943  p = j;
944  val = L->row_idx[p];
945  }
946  }
947  Int temp_index = L->row_idx[L->col_ptr[i]];
948  Entry temp_entry = L->val[L->col_ptr[i]];
949  L->row_idx[L->col_ptr[i]] = val;
950  L->val[L->col_ptr[i]] = L->val[p];
951  L->row_idx[p] = temp_index;
952  L->val[p] = temp_entry;
953  }//end for all columns
954 
955 
956  /* Sort CSC U --- just make sure max is in right location*/
957  for(i = 0 ; i < U->ncol; i++)
958  {
959  p = U->col_ptr[i+1]-1;
960  val = U->row_idx[p];
961 
962  for(j = U->col_ptr[i]; j < (U->col_ptr[i+1]-1); j++)
963  {
964  if(U->row_idx[j] > val)
965  {
966  p = j;
967  val = U->row_idx[p];
968  }
969  }
970  Int temp_index = U->row_idx[U->col_ptr[i+1]-1];
971  Entry temp_entry = U->val[U->col_ptr[i+1]-1];
972  U->row_idx[U->col_ptr[i+1]-1] = val;
973  U->val[U->col_ptr[i+1]-1] = U->val[p];
974  U->row_idx[p] = temp_index;
975  U->val[p] = temp_entry;
976  }//end for all columns
977 
978  return 0;
979  }
980 
981  template <class Int, class Entry>
982  Entry* BaskerClassic <Int, Entry>::entry_realloc(Entry *old , Int old_size, Int new_size)
983  {
984  Entry *new_entry = new Entry[new_size];
985  for(Int i = 0; i < old_size; i++)
986  {
987  /*Assumption that Entry was own defined copy constructor*/
988  new_entry[i] = old[i];
989  }
990  delete[] old;
991  return new_entry;
992 
993 
994  }
995  template <class Int, class Entry>
996  Int* BaskerClassic <Int, Entry>::int_realloc(Int *old, Int old_size, Int new_size)
997  {
998  Int *new_int = new Int[new_size];
999  for(Int i =0; i < old_size; i++)
1000  {
1001  /*Assumption that Int was own defined copy constructor*/
1002  new_int[i] = old[i];
1003  }
1004  delete[] old;
1005  return new_int;
1006 
1007  }
1008 
1009 
1010 }//end namespace
1011 #endif