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  //free factor
215  if(been_fact)
216  {
217  free_factor();
218  //BASKERFREE(pinv);
219  delete [] pinv;
220  }
221 
222  /*Creating space for L and U*/
223  L->nrow = nrow;
224  L->ncol = ncol;
225  if(L->nnz == 0)
226  {
227  L->nnz = 2*A->nnz;
228  }
229  //L->col_ptr = (Int *) BASKERCALLOC(ncol+1, sizeof(Int));
230  L->col_ptr = new Int[ncol+1]();
231  //L->row_idx = (Int *) BASKERCALLOC(L->nnz, sizeof(Int));
232  L->row_idx = new Int[L->nnz]();
233  //L->val = (Entry *) BASKERCALLOC(L->nnz, sizeof(Entry));
234  L->val = new Entry[L->nnz]();
235 
236  U->nrow = nrow;
237  U->ncol = ncol;
238  if(U->nnz == 0)
239  {
240  U->nnz = 2*A->nnz;
241  }
242  //U->col_ptr = (Int *) BASKERCALLOC(ncol+1, sizeof(Int));
243  U->col_ptr = new Int[ncol+1]();
244  //U->row_idx = (Int *) BASKERCALLOC(U->nnz, sizeof(Int));
245  U->row_idx = new Int[U->nnz]();
246  //U->val = (Entry *) BASKERCALLOC(U->nnz, sizeof(Entry));
247  U->val = new Entry[U->nnz]();
248 
249  if((L->col_ptr == nullptr) || (L->row_idx == nullptr) || (L->val == nullptr) ||
250  (U->col_ptr == nullptr) || (U->row_idx == nullptr) || (U->val == nullptr))
251  {
252  ierr = -1;
253  return ierr;
254  }
255  /*End creating space for L and U*/
256 
257  /*Creating working space*/
258  Int *color, *pattern, *stack;
259  Entry *X;
260  color = new Int[ncol]();
261  pattern = new Int[nrow]();
262  stack = new Int[2*nrow]();
263  //X = (Entry *) BASKERCALLOC(2*nrow, sizeof(Entry));
264  X = new Entry[2*nrow]();
265  //pinv = (Int * ) BASKERCALLOC(ncol+1, sizeof(Int)); //Note extra pad
266  pinv = new Int[ncol+1]();
267 
268 
269  if( (color == nullptr) || (pattern == nullptr) || (stack == nullptr) || (X == nullptr) || (pinv == nullptr) )
270  {
271  ierr = -2;
272  return ierr;
273  }
274 
275  /*End creating working space */
276 
277  /*Defining Variables Used*/
278  Int i, j, k;
279  Int top, top1, maxindex, t; // j1, j2;
280  Int lnnz, unnz, xnnz, lcnt, ucnt;
281  Int cu_ltop, cu_utop;
282  Int pp, p2, p;
283  Int newsize;
284  Entry pivot, value, xj;
285  Entry absv, maxv;
286 
287  cu_ltop = 0;
288  cu_utop = 0;
289  top = ncol;
290  top1 = ncol;
291  lnnz = 0; //real found lnnz
292  unnz = 0; //real found unnz
293 
294  for(k = 0 ; k < ncol; k++)
295  {
296  pinv[k] = ncol;
297  }
298 
299  /*For all columns in A .... */
300  for (k = 0; k < ncol; k++)
301  {
302 
303 #ifdef BASKER_DEBUG
304  cout << "k = " << k << endl;
305 #endif
306 
307  value = 0.0;
308  pivot = 0.0;
309  maxindex = ncol;
310  //j1 = 0;
311  //j2 = 0;
312  lcnt = 0;
313  ucnt = 0;
314 
315 #ifdef BASKER_DEBUG
316  BASKERASSERT (top == ncol);
317 
318  for(i = 0; i < nrow; i++)
319  {
320  BASKERASSERT(X[i] == (Entry)0);
321  }
322  for(i = 0; i < ncol; i++)
323  {
324  BASKERASSERT(color[i] == 0);
325  }
326 #endif
327  /* Reachability for every nonzero in Ak */
328  for( i = col_ptr[k]; i < col_ptr[k+1]; i++)
329  {
330  j = row_idx[i];
331  X[j] = val[i];
332 
333  if(color[j] == 0)
334  {
335  //do dfs
336  basker_dfs(nrow, j, L->row_idx, L->col_ptr, color, pattern, &top, pinv, stack);
337 
338  }
339 
340  }//end reachable
341 
342  xnnz = ncol - top;
343 #ifdef BASKER_DEBUG
344  cout << top << endl;
345  cout << ncol << endl;
346  cout << xnnz << endl;
347  //BASKERASSERT(xnnz <= nrow);
348 #endif
349  /*Lx = b where x will be the column k in L and U*/
350  top1 = top;
351  for(pp = 0; pp < xnnz; pp++)
352  {
353  j = pattern[top1++];
354  color[j] = 0;
355  t = pinv[j];
356 
357  if(t!=ncol) //it has not been assigned
358  {
359  xj = X[j];
360  p2 = L->col_ptr[t+1];
361  for(p = L->col_ptr[t]+1; p < p2; p++)
362  {
363  X[L->row_idx[p]] -= L->val[p] * xj;
364  }//over all rows
365  }
366 
367  }
368 
369  /*get the pivot*/
370  maxv = 0.0;
371  for(i = top; i < nrow; i++)
372  {
373  j = pattern[i];
374  t = pinv[j];
375  value = X[j];
376  /*note may want to change this to traits*/
377  //absv = (value < 0.0 ? -value : value);
378  absv = BASKER_ScalarTraits<Entry>::approxABS(value);
379 
380  if(t == ncol)
381  {
382  lcnt++;
383  if( BASKER_ScalarTraits<Entry>::gt(absv , maxv))
384  //if(absv > BASKER_ScalarTraits<Entry>::approxABS(maxv))
385  {
386  maxv = absv;
387  pivot = value;
388  maxindex= j;
389  }
390  }
391  }
392  ucnt = nrow - top - lcnt + 1;
393 
394  if(maxindex == ncol || pivot == ((Entry)0))
395  {
396  cout << "Matrix is singular at index: " << maxindex << " pivot: " << pivot << endl;
397  ierr = maxindex;
398  return ierr;
399  }
400 
401  pinv[maxindex] = k;
402 #ifdef BASKER_DEBUG
403  if(maxindex != k )
404  {
405  cout << "Permuting pivot: " << k << " for row: " << maxindex << endl;
406  }
407 #endif
408 
409  if(lnnz + lcnt >= L->nnz)
410  {
411 
412  newsize = L->nnz * 1.1 + 2*nrow + 1;
413 #ifdef BASKER_DEBUG
414  cout << "Out of memory -- Reallocating. Old Size: " << L->nnz << " New Size: " << newsize << endl;
415 #endif
416  //L->row_idx = (Int *) BASKERREALLOC(L->row_idx, newsize*sizeof(Int));
417  L->row_idx = int_realloc(L->row_idx , L->nnz, newsize);
418  if(!(L->row_idx))
419  {
420  cout << "WARNING: Cannot Realloc Memory" << endl;
421  ierr = -3;
422  return ierr;
423  }
424  //L->val = (Entry *) BASKERREALLOC(L->val, newsize*sizeof(Entry));
425  L->val = entry_realloc(L->val, L->nnz, newsize);
426  if(!(L->val))
427  {
428  cout << "WARNING: Cannot Realloc Memory" << endl;
429  ierr = -3;
430  return ierr;
431  }
432  L->nnz = newsize;
433 
434  }//realloc if L is out of memory
435 
436  if(unnz + ucnt >= U->nnz)
437  {
438  newsize = U->nnz*1.1 + 2*nrow + 1;
439 #ifdef BASKER_DEBUG
440  cout << "Out of memory -- Reallocating. Old Size: " << U->nnz << " New Size: " << newsize << endl;
441 #endif
442  //U->row_idx = (Int *) BASKERREALLOC(U->row_idx, newsize*sizeof(Int));
443  U->row_idx = int_realloc(U->row_idx, U->nnz, newsize);
444  if(!(U->row_idx))
445  {
446  cout << "WARNING: Cannot Realloc Memory" << endl;
447  ierr = -3;
448  return ierr;
449  }
450 
451  //U->val = (Entry *) BASKERREALLOC(U->val, newsize*sizeof(Entry));
452  U->val = entry_realloc(U->val, U->nnz, newsize);
453  if(!(U->val))
454  {
455  cout << "WARNING: Cannot Realloc Memory" << endl;
456  ierr = -3;
457  return ierr;
458  }
459  U->nnz = newsize;
460  }//realloc if U is out of memory
461 
462  //L->col_ptr[lnnz] = maxindex;
463  L->row_idx[lnnz] = maxindex;
464  L->val[lnnz] = 1.0;
465  lnnz++;
466 
467  Entry last_v_temp = 0;
468 
469  for(i = top; i < nrow; i++)
470  {
471  j = pattern[i];
472  t = pinv[j];
473 
474  /* check for numerical cancellations */
475 
476 
477  if(X[j] != ((Entry)0))
478  {
479 
480  if(t != ncol)
481  {
482  if(unnz >= U->nnz)
483  {
484  cout << "BASKER: Insufficent memory for U" << endl;
485  ierr = -3;
486  return ierr;
487  }
488  if(t < k)
489  //if(true)
490  {
491  U->row_idx[unnz] = pinv[j];
492  U->val[unnz] = X[j];
493  unnz++;
494  }
495  else
496  {
497 
498  last_v_temp = X[j];
499  //cout << "Called. t: " << t << "Val: " << last_v_temp << endl;
500  }
501 
502  }
503  else if (t == ncol)
504  {
505  if(lnnz >= L->nnz)
506  {
507  cout << "BASKER: Insufficent memory for L" << endl;
508  ierr = -3;
509  return ierr;
510  }
511 
512  L->row_idx[lnnz] = j;
513  //L->val[lnnz] = X[j]/pivot;
514  L->val[lnnz] = BASKER_ScalarTraits<Entry>::divide(X[j],pivot);
515  lnnz++;
516 
517  }
518 
519  }
520 
521 
522  X[j] = 0;
523 
524  }
525  //cout << "Value added at end: " << last_v_temp << endl;
526  U->row_idx[unnz] = k;
527  U->val[unnz] = last_v_temp;
528  unnz++;
529 
530  xnnz = 0;
531  top = ncol;
532 
533  L->col_ptr[k] = cu_ltop;
534  L->col_ptr[k+1] = lnnz;
535  cu_ltop = lnnz;
536 
537  U->col_ptr[k] = cu_utop;
538  U->col_ptr[k+1] = unnz;
539  cu_utop = unnz;
540 
541  } //end for every column
542 
543 #ifdef BASKER_DEBUG
544  /*Print out found L and U*/
545  for(k = 0; k < lnnz; k++)
546  {
547  printf("L[%d]=%g" , k , L->val[k]);
548  }
549  cout << endl;
550  for(k = 0; k < lnnz; k++)
551  {
552  printf("Li[%d]=%d", k, L->row_idx[k]);
553  }
554  cout << endl;
555  for(k = 0; k < nrow; k++)
556  {
557  printf("p[%d]=%d", k, pinv[k]);
558  }
559  cout << endl;
560  cout << endl;
561 
562  for(k = 0; k < ncol; k++)
563  {
564  printf("Up[%d]=%d", k, U->col_ptr[k]);
565  }
566  cout << endl;
567 
568  for(k = 0; k < unnz; k++)
569  {
570  printf("U[%d]=%g" , k , U->val[k]);
571  }
572  cout << endl;
573  for(k = 0; k < unnz; k++)
574  {
575  printf("Ui[%d]=%d", k, U->row_idx[k]);
576  }
577  cout << endl;
578 
579 
580 #endif
581  /* Repermute */
582  for( i = 0; i < ncol; i++)
583  {
584  for(k = L->col_ptr[i]; k < L->col_ptr[i+1]; k++)
585  {
586  //L->row_idx[k] = pinv[L->row_idx[k]];
587  }
588  }
589  //Max sure correct location of min in L and max in U for CSC format//
590  //Speeds up tri-solve//
591  //sort_factors();
592 
593 #ifdef BASKER_DEBUG
594  cout << "After Permuting" << endl;
595  for(k = 0; k < lnnz; k++)
596  {
597  printf("Li[%d]=%d", k, L->row_idx[k]);
598  }
599  cout << endl;
600 #endif
601 
602  // Cleanup workspace allocations
603  delete [] X;
604  delete [] color;
605  delete [] pattern;
606  delete [] stack;
607 
608  actual_lnnz = lnnz;
609  actual_unnz = unnz;
610 
611  been_fact = true;
612  return 0;
613  }//end factor
614 
615 
616  template <class Int, class Entry>
617  Int BaskerClassic<Int, Entry>::get_NnzL()
618  {
619  return actual_lnnz;
620  }
621 
622  template <class Int, class Entry>
623  Int BaskerClassic<Int, Entry>::get_NnzU()
624  {
625  return actual_unnz;
626  }
627 
628  template <class Int, class Entry>
629  Int BaskerClassic<Int, Entry>::get_NnzLU()
630  {
631  return (actual_lnnz + actual_unnz);
632  }
633 
634  template <class Int, class Entry>
635  int BaskerClassic<Int, Entry>::returnL(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
636  {
637  int i;
638  *dim = L->nrow;
639  *nnz = L->nnz;
640 
641  /*Does a bad copy*/
642 
643  //*col_ptr = (Int *) BASKERCALLOC(L->nrow+1, sizeof(Int));
644  *col_ptr = new Int[L->nrow+1];
645  //*row_idx = (Int *) BASKERCALLOC(L->nnz, sizeof(Int));
646  *row_idx = new Int[L->nnz];
647  //*val = (Entry *) BASKERCALLOC(L->nnz, sizeof(Entry));
648  *val = new Entry[L->nnz];
649 
650  if( (*col_ptr == nullptr) || (*row_idx == nullptr) || (*val == nullptr) )
651  {
652  return -1;
653  }
654 
655  for(i = 0; i < L->nrow+1; i++)
656  {
657  (*col_ptr)[i] = L->col_ptr[i];
658  }
659 
660  for(i = 0; i < actual_lnnz; i++)
661  {
662  (*row_idx)[i] = pinv[L->row_idx[i]];
663  (*val)[i] = L->val[i];
664  }
665  return 0;
666 
667  }
668 
669  template <class Int, class Entry>
670  int BaskerClassic<Int, Entry>::returnU(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
671  {
672  int i;
673  *dim = U->nrow;
674  *nnz = U->nnz;
675  /*Does a bad copy*/
676  //*col_ptr = (Int *) BASKERCALLOC(U->nrow+1, sizeof(Int));
677  *col_ptr = new Int[U->nrow+1];
678  //*row_idx = (Int *) BASKERCALLOC(U->nnz, sizeof(Int));
679  *row_idx = new Int[U->nnz];
680  //*val = (Entry *) BASKERCALLOC(U->nnz, sizeof(Entry));
681  *val = new Entry[U->nnz];
682 
683  if( (*col_ptr == nullptr) || (*row_idx == nullptr) || (*val == nullptr) )
684  {
685  return -1;
686  }
687 
688  for(i = 0; i < U->nrow+1; i++)
689  {
690  (*col_ptr)[i] = U->col_ptr[i];
691  }
692  for(i = 0; i < actual_unnz; i++)
693  {
694  (*row_idx)[i] = U->row_idx[i];
695  (*val)[i] = U->val[i];
696  }
697  return 0;
698  }
699 
700  template <class Int, class Entry>
701  int BaskerClassic<Int, Entry>::returnP(Int** p)
702  {
703  Int i;
704  //*p = (Int *) BASKERCALLOC(A->nrow, sizeof(Int));
705  *p = new Int[A->nrow];
706 
707  if( (*p == nullptr ) )
708  {
709  return -1;
710  }
711 
712  for(i = 0; i < A->nrow; i++)
713  {
714  (*p)[pinv[i]] = i; //Matlab perm-style
715  }
716  return 0;
717  }
718 
719  template <class Int, class Entry>
720  void BaskerClassic<Int, Entry>::free_factor()
721  {
722  //BASKERFREE L
723  //BASKERFREE(L->col_ptr);
724  delete[] L->col_ptr;
725  //BASKERFREE(L->row_idx);
726  delete[] L->row_idx;
727  //BASKERFREE(L->val);
728  delete[] L->val;
729 
730  //BASKERFREE U
731  //BASKERFREE(U->col_ptr);
732  delete[] U->col_ptr;
733  //BASKERFREE(U->row_idx);
734  delete[] U->row_idx;
735  //BASKERFREE(U->val);
736  delete[] U->val;
737 
738  been_fact = false;
739  }
740  template <class Int, class Entry>
741  void BaskerClassic<Int, Entry>::free_perm_matrix()
742  {
743  //BASKERFREE(A->col_ptr);
744  //BASKERFREE(A->row_idx);
745  //BASKERFREE(A->val);
746  }
747 
748  template <class Int, class Entry>
749  int BaskerClassic<Int, Entry>::solveMultiple(Int nrhs, Entry *b, Entry *x)
750  {
751  Int i;
752  for(i = 0; i < nrhs; i++)
753  {
754  Int k = i*A->nrow;
755  int result = solve(&(b[k]), &(x[k]));
756  if(result != 0)
757  {
758  cout << "Error in Solving \n";
759  return result;
760  }
761  }
762  return 0;
763  }
764 
765 
766  template <class Int, class Entry>
767  int BaskerClassic<Int, Entry>::solve(Entry* b, Entry* x)
768  {
769 
770  if(!been_fact)
771  {
772  return -10;
773  }
774  //Entry* temp = (Entry *)BASKERCALLOC(A->nrow, sizeof(Entry));
775  Entry* temp = new Entry[A->nrow]();
776  Int i;
777  int result = 0;
778  for(i = 0 ; i < A->ncol; i++)
779  {
780  Int k = pinv[i];
781  x[k] = b[i];
782  }
783 
784  result = low_tri_solve_csc(L->nrow, L->col_ptr, L->row_idx, L->val, temp, x);
785  if(result == 0)
786  {
787  result = up_tri_solve_csc(U->nrow, U->col_ptr, U->row_idx, U->val, x, temp);
788  }
789 
790 
791  //BASKERFREE(temp);
792  delete[] temp;
793  return 0;
794  }
795 
796  template < class Int, class Entry>
797  int BaskerClassic<Int, Entry>::low_tri_solve_csc( Int n, Int *col_ptr, Int *row_idx, Entry* val, Entry *x, Entry *b)
798  {
799  Int i, j;
800  /*for each column*/
801  for(i = 0; i < n ; i++)
802  {
803 #ifdef BASKER_DEBUG
804  BASKERASSERT(val[col_ptr[i]] != (Entry)0);
805 #else
806  if(val[col_ptr[i]] == (Entry) 0)
807  {
808  return i;
809  }
810 #endif
811  x[i] = BASKER_ScalarTraits<Entry>::divide(b[i], val[col_ptr[i]]);
812 
813  for(j = col_ptr[i]+1; j < (col_ptr[i+1]); j++) //update all rows
814  {
815  b[pinv[row_idx[j]]] = b[pinv[row_idx[j]]] - (val[j]*x[i]);
816  }
817  }
818  return 0;
819  }
820 
821  template < class Int, class Entry>
822  int BaskerClassic<Int, Entry>::up_tri_solve_csc( Int n, Int *col_ptr, Int *row_idx, Entry *val, Entry *x, Entry *b)
823  {
824  Int i, j;
825  /*for each column*/
826  for(i = n; i > 1 ; i--)
827  {
828  int ii = i-1;
829 #ifdef BASKER_DEBUG
830  BASKERASSERT(val[col_ptr[i]-1] != (Entry)0);
831 #else
832  if(val[col_ptr[i]-1] == (Entry) 0)
833  {
834  cout << "Dig(" << i << ") = " << val[col_ptr[i]-1] << endl;
835  return i;
836  }
837 #endif
838  //x[ii] = b[ii]/val[col_ptr[i]-1]; //diag
839  x[ii] = BASKER_ScalarTraits<Entry>::divide(b[ii],val[col_ptr[i]-1]);
840 
841  for(j = (col_ptr[i]-2); j >= (col_ptr[ii]); j--)
842  {
843  b[row_idx[j]] = b[row_idx[j]] - (val[j]*x[ii]);
844  }
845  }
846  //x[0] = b[0]/val[col_ptr[1]-1];
847  x[0] = BASKER_ScalarTraits<Entry>::divide(b[0],val[col_ptr[1]-1]);
848  return 0;
849  }
850 
851  template <class Int, class Entry>
852  int BaskerClassic<Int, Entry>::preorder(Int *row_perm, Int *col_perm)
853  {
854 
855  basker_matrix <Int, Entry> *B;
856  B = new basker_matrix<Int, Entry>;
857  B->nrow = A->nrow;
858  B->ncol = A->ncol;
859  B->nnz = A->nnz;
860  B->col_ptr = (Int *) BASKERCALLOC(A->ncol + 1, sizeof(Int));
861  B->row_idx = (Int *) BASKERCALLOC(A->nnz, sizeof(Int));
862  B->val = (Entry *) BASKERCALLOC(A->val, sizeof(Int));
863 
864  if( (B->col_ptr == nullptr) || (B->row_idx == nullptr) || (B->val == nullptr) )
865  {
866  perm_flag = false;
867  return -1;
868  }
869 
870  /* int resultcol = (unused) */ (void) permute_column(col_perm, B);
871  /* int resultrow = (unused) */ (void) permute_row(row_perm, B);
872 
873  /*Note: the csc matrices of A are the problem of the user
874  therefore we will not free them*/
875  A->col_ptr = B->col_ptr;
876  A->row_idx = B->row_idx;
877  A->val = A->val;
878 
879  perm_flag = true; /*Now we will free A at the end*/
880 
881  return 0;
882  }
883 
884  template <class Int, class Entry>
885  int BaskerClassic <Int, Entry>::permute_column(Int *p, basker_matrix<Int,Entry> *B)
886  {
887  /*p(i) contains the destination of row i in the permuted matrix*/
888  Int i,j, ii, jj;
889 
890  /*Determine column pointer of output matrix*/
891  for(j=0; j < B->ncol; j++)
892  {
893  i = p[j];
894  B->col_ptr[i+1] = A->col_ptr[j+1] - A->col_ptr[j];
895  }
896  /*get pointers from lengths*/
897  B->col_ptr[0] = 0;
898  for(j=0; j < B->ncol; j++)
899  {
900  B->col_ptr[j+1] = B->col_ptr[j+1] + B->col_ptr[j];
901  }
902 
903  /*copy idxs*/
904  Int k, ko;
905  for(ii = 0 ; ii < B->ncol; ii++)
906  {// old colum ii new column p[ii] k->pointer
907  ko = B->col_ptr(p[ii]);
908  for(k = A->col_ptr[ii]; k < A->col_ptr[ii+1]; k++)
909  {
910  B->row_index[ko] = A->row_index[k];
911  B->val[ko] = A->val[ko];
912  ko++;
913  }
914  }
915  return 0;
916  }
917 
918  template <class Int, class Entry>
919  int BaskerClassic <Int, Entry>::permute_row(Int *p, basker_matrix<Int,Entry> *B)
920  {
921  Int k,i;
922  for(k=0; k < A->nnz; k++)
923  {
924  B->row_idx[k] = p[A->row_idx[k]];
925  }
926  return 0;
927  }
928 
929  template <class Int, class Entry>
930  int BaskerClassic <Int, Entry>::sort_factors()
931  {
932 
933  /*Sort CSC of L - just make sure min_index is in lowest position*/
934  Int i, j;
935  Int p;
936  Int val;
937  for(i = 0 ; i < L->ncol; i++)
938  {
939  p = L->col_ptr[i];
940  val = L->row_idx[p];
941 
942  for(j = L->col_ptr[i]+1; j < (L->col_ptr[i+1]); j++)
943  {
944  if(L->row_idx[j] < val)
945  {
946  p = j;
947  val = L->row_idx[p];
948  }
949  }
950  Int temp_index = L->row_idx[L->col_ptr[i]];
951  Entry temp_entry = L->val[L->col_ptr[i]];
952  L->row_idx[L->col_ptr[i]] = val;
953  L->val[L->col_ptr[i]] = L->val[p];
954  L->row_idx[p] = temp_index;
955  L->val[p] = temp_entry;
956  }//end for all columns
957 
958 
959  /* Sort CSC U --- just make sure max is in right location*/
960  for(i = 0 ; i < U->ncol; i++)
961  {
962  p = U->col_ptr[i+1]-1;
963  val = U->row_idx[p];
964 
965  for(j = U->col_ptr[i]; j < (U->col_ptr[i+1]-1); j++)
966  {
967  if(U->row_idx[j] > val)
968  {
969  p = j;
970  val = U->row_idx[p];
971  }
972  }
973  Int temp_index = U->row_idx[U->col_ptr[i+1]-1];
974  Entry temp_entry = U->val[U->col_ptr[i+1]-1];
975  U->row_idx[U->col_ptr[i+1]-1] = val;
976  U->val[U->col_ptr[i+1]-1] = U->val[p];
977  U->row_idx[p] = temp_index;
978  U->val[p] = temp_entry;
979  }//end for all columns
980 
981  return 0;
982  }
983 
984  template <class Int, class Entry>
985  Entry* BaskerClassic <Int, Entry>::entry_realloc(Entry *old , Int old_size, Int new_size)
986  {
987  Entry *new_entry = new Entry[new_size];
988  for(Int i = 0; i < old_size; i++)
989  {
990  /*Assumption that Entry was own defined copy constructor*/
991  new_entry[i] = old[i];
992  }
993  delete[] old;
994  return new_entry;
995 
996 
997  }
998  template <class Int, class Entry>
999  Int* BaskerClassic <Int, Entry>::int_realloc(Int *old, Int old_size, Int new_size)
1000  {
1001  Int *new_int = new Int[new_size];
1002  for(Int i =0; i < old_size; i++)
1003  {
1004  /*Assumption that Int was own defined copy constructor*/
1005  new_int[i] = old[i];
1006  }
1007  delete[] old;
1008  return new_int;
1009 
1010  }
1011 
1012 
1013 }//end namespace
1014 #endif