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