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