29 #ifndef BASKER_DEF_HPP
30 #define BASKER_DEF_HPP
32 #include "basker_decl.hpp"
33 #include "basker_scalartraits.hpp"
45 namespace BaskerClassicNS{
47 template <
class Int,
class Entry>
48 BaskerClassic<Int, Entry>::BaskerClassic()
52 A =
new basker_matrix<Int,Entry>;
55 L =
new basker_matrix<Int, Entry>;
59 U =
new basker_matrix<Int,Entry>;
70 template <
class Int,
class Entry>
71 BaskerClassic<Int, Entry>::BaskerClassic(Int nnzL, Int nnzU)
75 A =
new basker_matrix<Int, Entry>;
77 L =
new basker_matrix<Int, Entry>;
80 U =
new basker_matrix<Int, Entry>;
91 template <
class Int,
class Entry>
92 BaskerClassic<Int, Entry>::~BaskerClassic()
114 template <
class Int,
class Entry>
115 int BaskerClassic<Int,Entry>:: basker_dfs
131 Int start, end, done, *store ;
136 bool has_elements =
true;
153 BASKERASSERT (color[j] == 1) ;
161 for ( i1 = start ; i1 < end ; i1++ )
175 pattern[--*top] = j ;
179 has_elements =
false;
188 std::cout <<
"Out of DFS: " << j << std::endl;
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)
200 BASKERASSERT(nrow > 0);
201 BASKERASSERT(ncol > 0);
202 BASKERASSERT(nnz > 0);
209 A->col_ptr = col_ptr;
210 A->row_idx = row_idx;
230 L->col_ptr =
new Int[ncol+1]();
232 L->row_idx =
new Int[L->nnz]();
234 L->val =
new Entry[L->nnz]();
243 U->col_ptr =
new Int[ncol+1]();
245 U->row_idx =
new Int[U->nnz]();
247 U->val =
new Entry[U->nnz]();
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))
258 Int *color, *pattern, *stack;
260 color =
new Int[ncol]();
261 pattern =
new Int[nrow]();
262 stack =
new Int[2*nrow]();
264 X =
new Entry[2*nrow]();
266 pinv =
new Int[ncol+1]();
269 if( (color ==
nullptr) || (pattern ==
nullptr) || (stack ==
nullptr) || (X ==
nullptr) || (pinv ==
nullptr) )
279 Int top, top1, maxindex, t;
280 Int lnnz, unnz, xnnz, lcnt, ucnt;
281 Int cu_ltop, cu_utop;
284 Entry pivot, value, xj;
294 for(k = 0 ; k < ncol; k++)
300 for (k = 0; k < ncol; k++)
304 cout <<
"k = " << k << endl;
316 BASKERASSERT (top == ncol);
318 for(i = 0; i < nrow; i++)
320 BASKERASSERT(X[i] == (Entry)0);
322 for(i = 0; i < ncol; i++)
324 BASKERASSERT(color[i] == 0);
328 for( i = col_ptr[k]; i < col_ptr[k+1]; i++)
336 basker_dfs(nrow, j, L->row_idx, L->col_ptr, color, pattern, &top, pinv, stack);
345 cout << ncol << endl;
346 cout << xnnz << endl;
351 for(pp = 0; pp < xnnz; pp++)
360 p2 = L->col_ptr[t+1];
361 for(p = L->col_ptr[t]+1; p < p2; p++)
363 X[L->row_idx[p]] -= L->val[p] * xj;
371 for(i = top; i < nrow; i++)
378 absv = BASKER_ScalarTraits<Entry>::approxABS(value);
383 if( BASKER_ScalarTraits<Entry>::gt(absv , maxv))
392 ucnt = nrow - top - lcnt + 1;
394 if(maxindex == ncol || pivot == ((Entry)0))
396 cout <<
"Matrix is singular at index: " << maxindex <<
" pivot: " << pivot << endl;
405 cout <<
"Permuting pivot: " << k <<
" for row: " << maxindex << endl;
409 if(lnnz + lcnt >= L->nnz)
412 newsize = L->nnz * 1.1 + 2*nrow + 1;
414 cout <<
"Out of memory -- Reallocating. Old Size: " << L->nnz <<
" New Size: " << newsize << endl;
417 L->row_idx = int_realloc(L->row_idx , L->nnz, newsize);
420 cout <<
"WARNING: Cannot Realloc Memory" << endl;
425 L->val = entry_realloc(L->val, L->nnz, newsize);
428 cout <<
"WARNING: Cannot Realloc Memory" << endl;
436 if(unnz + ucnt >= U->nnz)
438 newsize = U->nnz*1.1 + 2*nrow + 1;
440 cout <<
"Out of memory -- Reallocating. Old Size: " << U->nnz <<
" New Size: " << newsize << endl;
443 U->row_idx = int_realloc(U->row_idx, U->nnz, newsize);
446 cout <<
"WARNING: Cannot Realloc Memory" << endl;
452 U->val = entry_realloc(U->val, U->nnz, newsize);
455 cout <<
"WARNING: Cannot Realloc Memory" << endl;
463 L->row_idx[lnnz] = maxindex;
467 Entry last_v_temp = 0;
469 for(i = top; i < nrow; i++)
477 if(X[j] != ((Entry)0))
484 cout <<
"BASKER: Insufficent memory for U" << endl;
491 U->row_idx[unnz] = pinv[j];
507 cout <<
"BASKER: Insufficent memory for L" << endl;
512 L->row_idx[lnnz] = j;
514 L->val[lnnz] = BASKER_ScalarTraits<Entry>::divide(X[j],pivot);
526 U->row_idx[unnz] = k;
527 U->val[unnz] = last_v_temp;
533 L->col_ptr[k] = cu_ltop;
534 L->col_ptr[k+1] = lnnz;
537 U->col_ptr[k] = cu_utop;
538 U->col_ptr[k+1] = unnz;
545 for(k = 0; k < lnnz; k++)
547 printf(
"L[%d]=%g" , k , L->val[k]);
550 for(k = 0; k < lnnz; k++)
552 printf(
"Li[%d]=%d", k, L->row_idx[k]);
555 for(k = 0; k < nrow; k++)
557 printf(
"p[%d]=%d", k, pinv[k]);
562 for(k = 0; k < ncol; k++)
564 printf(
"Up[%d]=%d", k, U->col_ptr[k]);
568 for(k = 0; k < unnz; k++)
570 printf(
"U[%d]=%g" , k , U->val[k]);
573 for(k = 0; k < unnz; k++)
575 printf(
"Ui[%d]=%d", k, U->row_idx[k]);
582 for( i = 0; i < ncol; i++)
584 for(k = L->col_ptr[i]; k < L->col_ptr[i+1]; k++)
594 cout <<
"After Permuting" << endl;
595 for(k = 0; k < lnnz; k++)
597 printf(
"Li[%d]=%d", k, L->row_idx[k]);
616 template <
class Int,
class Entry>
617 Int BaskerClassic<Int, Entry>::get_NnzL()
622 template <
class Int,
class Entry>
623 Int BaskerClassic<Int, Entry>::get_NnzU()
628 template <
class Int,
class Entry>
629 Int BaskerClassic<Int, Entry>::get_NnzLU()
631 return (actual_lnnz + actual_unnz);
634 template <
class Int,
class Entry>
635 int BaskerClassic<Int, Entry>::returnL(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
644 *col_ptr =
new Int[L->nrow+1];
646 *row_idx =
new Int[L->nnz];
648 *val =
new Entry[L->nnz];
650 if( (*col_ptr ==
nullptr) || (*row_idx ==
nullptr) || (*val ==
nullptr) )
655 for(i = 0; i < L->nrow+1; i++)
657 (*col_ptr)[i] = L->col_ptr[i];
660 for(i = 0; i < actual_lnnz; i++)
662 (*row_idx)[i] = pinv[L->row_idx[i]];
663 (*val)[i] = L->val[i];
669 template <
class Int,
class Entry>
670 int BaskerClassic<Int, Entry>::returnU(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
677 *col_ptr =
new Int[U->nrow+1];
679 *row_idx =
new Int[U->nnz];
681 *val =
new Entry[U->nnz];
683 if( (*col_ptr ==
nullptr) || (*row_idx ==
nullptr) || (*val ==
nullptr) )
688 for(i = 0; i < U->nrow+1; i++)
690 (*col_ptr)[i] = U->col_ptr[i];
692 for(i = 0; i < actual_unnz; i++)
694 (*row_idx)[i] = U->row_idx[i];
695 (*val)[i] = U->val[i];
700 template <
class Int,
class Entry>
701 int BaskerClassic<Int, Entry>::returnP(Int** p)
705 *p =
new Int[A->nrow];
707 if( (*p ==
nullptr ) )
712 for(i = 0; i < A->nrow; i++)
719 template <
class Int,
class Entry>
720 void BaskerClassic<Int, Entry>::free_factor()
740 template <
class Int,
class Entry>
741 void BaskerClassic<Int, Entry>::free_perm_matrix()
748 template <
class Int,
class Entry>
749 int BaskerClassic<Int, Entry>::solveMultiple(Int nrhs, Entry *b, Entry *x)
752 for(i = 0; i < nrhs; i++)
755 int result = solve(&(b[k]), &(x[k]));
758 cout <<
"Error in Solving \n";
766 template <
class Int,
class Entry>
767 int BaskerClassic<Int, Entry>::solve(Entry* b, Entry* x)
775 Entry* temp =
new Entry[A->nrow]();
778 for(i = 0 ; i < A->ncol; i++)
784 result = low_tri_solve_csc(L->nrow, L->col_ptr, L->row_idx, L->val, temp, x);
787 result = up_tri_solve_csc(U->nrow, U->col_ptr, U->row_idx, U->val, x, temp);
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)
801 for(i = 0; i < n ; i++)
804 BASKERASSERT(val[col_ptr[i]] != (Entry)0);
806 if(val[col_ptr[i]] == (Entry) 0)
811 x[i] = BASKER_ScalarTraits<Entry>::divide(b[i], val[col_ptr[i]]);
813 for(j = col_ptr[i]+1; j < (col_ptr[i+1]); j++)
815 b[pinv[row_idx[j]]] = b[pinv[row_idx[j]]] - (val[j]*x[i]);
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)
826 for(i = n; i > 1 ; i--)
830 BASKERASSERT(val[col_ptr[i]-1] != (Entry)0);
832 if(val[col_ptr[i]-1] == (Entry) 0)
834 cout <<
"Dig(" << i <<
") = " << val[col_ptr[i]-1] << endl;
839 x[ii] = BASKER_ScalarTraits<Entry>::divide(b[ii],val[col_ptr[i]-1]);
841 for(j = (col_ptr[i]-2); j >= (col_ptr[ii]); j--)
843 b[row_idx[j]] = b[row_idx[j]] - (val[j]*x[ii]);
847 x[0] = BASKER_ScalarTraits<Entry>::divide(b[0],val[col_ptr[1]-1]);
851 template <
class Int,
class Entry>
852 int BaskerClassic<Int, Entry>::preorder(Int *row_perm, Int *col_perm)
855 basker_matrix <Int, Entry> *B;
856 B =
new basker_matrix<Int, Entry>;
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));
864 if( (B->col_ptr ==
nullptr) || (B->row_idx ==
nullptr) || (B->val ==
nullptr) )
870 (void) permute_column(col_perm, B);
871 (void) permute_row(row_perm, B);
875 A->col_ptr = B->col_ptr;
876 A->row_idx = B->row_idx;
884 template <
class Int,
class Entry>
885 int BaskerClassic <Int, Entry>::permute_column(Int *p, basker_matrix<Int,Entry> *B)
891 for(j=0; j < B->ncol; j++)
894 B->col_ptr[i+1] = A->col_ptr[j+1] - A->col_ptr[j];
898 for(j=0; j < B->ncol; j++)
900 B->col_ptr[j+1] = B->col_ptr[j+1] + B->col_ptr[j];
905 for(ii = 0 ; ii < B->ncol; ii++)
907 ko = B->col_ptr(p[ii]);
908 for(k = A->col_ptr[ii]; k < A->col_ptr[ii+1]; k++)
910 B->row_index[ko] = A->row_index[k];
911 B->val[ko] = A->val[ko];
918 template <
class Int,
class Entry>
919 int BaskerClassic <Int, Entry>::permute_row(Int *p, basker_matrix<Int,Entry> *B)
922 for(k=0; k < A->nnz; k++)
924 B->row_idx[k] = p[A->row_idx[k]];
929 template <
class Int,
class Entry>
930 int BaskerClassic <Int, Entry>::sort_factors()
937 for(i = 0 ; i < L->ncol; i++)
942 for(j = L->col_ptr[i]+1; j < (L->col_ptr[i+1]); j++)
944 if(L->row_idx[j] < val)
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;
960 for(i = 0 ; i < U->ncol; i++)
962 p = U->col_ptr[i+1]-1;
965 for(j = U->col_ptr[i]; j < (U->col_ptr[i+1]-1); j++)
967 if(U->row_idx[j] > val)
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;
984 template <
class Int,
class Entry>
985 Entry* BaskerClassic <Int, Entry>::entry_realloc(Entry *old , Int old_size, Int new_size)
987 Entry *new_entry =
new Entry[new_size];
988 for(Int i = 0; i < old_size; i++)
991 new_entry[i] = old[i];
998 template <
class Int,
class Entry>
999 Int* BaskerClassic <Int, Entry>::int_realloc(Int *old, Int old_size, Int new_size)
1001 Int *new_int =
new Int[new_size];
1002 for(Int i =0; i < old_size; i++)
1005 new_int[i] = old[i];