29 #ifndef BASKER_DEF_HPP 30 #define BASKER_DEF_HPP 32 #include "basker_decl.hpp" 33 #include "basker_scalartraits.hpp" 47 template <
class Int,
class Entry>
48 Basker<Int, Entry>::Basker()
52 A =
new basker_matrix<Int,Entry>;
55 L =
new basker_matrix<Int, Entry>;
59 U =
new basker_matrix<Int,Entry>;
67 template <
class Int,
class Entry>
68 Basker<Int, Entry>::Basker(Int nnzL, Int nnzU)
72 A =
new basker_matrix<Int, Entry>;
74 L =
new basker_matrix<Int, Entry>;
77 U =
new basker_matrix<Int, Entry>;
85 template <
class Int,
class Entry>
86 Basker<Int, Entry>::~Basker()
109 template <
class Int,
class Entry>
110 int Basker<Int,Entry>:: basker_dfs
126 Int start, end, done, *store ;
131 bool has_elements =
true;
148 ASSERT (color[j] == 1) ;
156 for ( i1 = start ; i1 < end ; i1++ )
170 pattern[--*top] = j ;
174 has_elements =
false;
183 std::cout <<
"Out of DFS: " << j << std::endl;
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)
204 A->col_ptr = col_ptr;
205 A->row_idx = row_idx;
217 L->col_ptr =
new Int[ncol+1]();
219 L->row_idx =
new Int[L->nnz]();
221 L->val =
new Entry[L->nnz]();
230 U->col_ptr =
new Int[ncol+1]();
232 U->row_idx =
new Int[U->nnz]();
234 U->val =
new Entry[U->nnz]();
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))
248 tptr =
new Int[(ncol)+(4*nrow)]();
250 X =
new Entry[2*nrow]();
252 pinv =
new Int[ncol+1]();
255 if( (tptr == NULL) || (X == NULL) || (pinv == NULL) )
265 Int *color, *pattern, *stack;
266 Int top, top1, maxindex, t;
267 Int lnnz, unnz, xnnz, lcnt, ucnt;
268 Int cu_ltop, cu_utop;
271 Entry pivot, value, xj;
290 for(k = 0 ; k < ncol; k++)
296 for (k = 0; k < ncol; k++)
300 cout <<
"k = " << k << endl;
312 ASSERT (top == ncol);
314 for(i = 0; i < nrow; i++)
316 ASSERT(X[i] == (Entry)0);
318 for(i = 0; i < ncol; i++)
320 ASSERT(color[i] == 0);
324 for( i = col_ptr[k]; i < col_ptr[k+1]; i++)
332 basker_dfs(nrow, j, L->row_idx, L->col_ptr, color, pattern, &top, pinv, stack);
341 cout << ncol << endl;
342 cout << xnnz << endl;
347 for(pp = 0; pp < xnnz; pp++)
356 p2 = L->col_ptr[t+1];
357 for(p = L->col_ptr[t]+1; p < p2; p++)
359 X[L->row_idx[p]] -= L->val[p] * xj;
367 for(i = top; i < nrow; i++)
374 absv = BASKER_ScalarTraits<Entry>::approxABS(value);
379 if( BASKER_ScalarTraits<Entry>::gt(absv , maxv))
388 ucnt = nrow - top - lcnt + 1;
390 if(maxindex == ncol || pivot == ((Entry)0))
392 cout <<
"Matrix is singular at index: " << maxindex <<
" pivot: " << pivot << endl;
401 cout <<
"Permuting pivot: " << k <<
" for row: " << maxindex << endl;
405 if(lnnz + lcnt >= L->nnz)
408 newsize = L->nnz * 1.1 + 2*nrow + 1;
410 cout <<
"Out of memory -- Reallocating. Old Size: " << L->nnz <<
" New Size: " << newsize << endl;
413 L->row_idx = int_realloc(L->row_idx , L->nnz, newsize);
416 cout <<
"WARNING: Cannot Realloc Memory" << endl;
421 L->val = entry_realloc(L->val, L->nnz, newsize);
424 cout <<
"WARNING: Cannot Realloc Memory" << endl;
432 if(unnz + ucnt >= U->nnz)
434 newsize = U->nnz*1.1 + 2*nrow + 1;
436 cout <<
"Out of memory -- Reallocating. Old Size: " << U->nnz <<
" New Size: " << newsize << endl;
439 U->row_idx = int_realloc(U->row_idx, U->nnz, newsize);
442 cout <<
"WARNING: Cannot Realloc Memory" << endl;
448 U->val = entry_realloc(U->val, U->nnz, newsize);
451 cout <<
"WARNING: Cannot Realloc Memory" << endl;
459 L->row_idx[lnnz] = maxindex;
463 Entry last_v_temp = 0;
465 for(i = top; i < nrow; i++)
473 if(X[j] != ((Entry)0))
480 cout <<
"BASKER: Insufficent memory for U" << endl;
487 U->row_idx[unnz] = pinv[j];
503 cout <<
"BASKER: Insufficent memory for L" << endl;
508 L->row_idx[lnnz] = j;
510 L->val[lnnz] = BASKER_ScalarTraits<Entry>::divide(X[j],pivot);
522 U->row_idx[unnz] = k;
523 U->val[unnz] = last_v_temp;
529 L->col_ptr[k] = cu_ltop;
530 L->col_ptr[k+1] = lnnz;
533 U->col_ptr[k] = cu_utop;
534 U->col_ptr[k+1] = unnz;
541 for(k = 0; k < lnnz; k++)
543 printf(
"L[%d]=%g" , k , L->val[k]);
546 for(k = 0; k < lnnz; k++)
548 printf(
"Li[%d]=%d", k, L->row_idx[k]);
551 for(k = 0; k < nrow; k++)
553 printf(
"p[%d]=%d", k, pinv[k]);
558 for(k = 0; k < ncol; k++)
560 printf(
"Up[%d]=%d", k, U->col_ptr[k]);
564 for(k = 0; k < unnz; k++)
566 printf(
"U[%d]=%g" , k , U->val[k]);
569 for(k = 0; k < unnz; k++)
571 printf(
"Ui[%d]=%d", k, U->row_idx[k]);
578 for( i = 0; i < ncol; i++)
580 for(k = L->col_ptr[i]; k < L->col_ptr[i+1]; k++)
590 cout <<
"After Permuting" << endl;
591 for(k = 0; k < lnnz; k++)
593 printf(
"Li[%d]=%d", k, L->row_idx[k]);
604 template <
class Int,
class Entry>
605 int Basker<Int, Entry>::returnL(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
614 *col_ptr =
new Int[L->nrow+1];
616 *row_idx =
new Int[L->nnz];
618 *val =
new Entry[L->nnz];
620 if( (*col_ptr == NULL) || (*row_idx == NULL) || (*val == NULL) )
625 for(i = 0; i < L->nrow+1; i++)
627 (*col_ptr)[i] = L->col_ptr[i];
630 for(i = 0; i < L->nnz; i++)
632 (*row_idx)[i] = pinv[L->row_idx[i]];
633 (*val)[i] = L->val[i];
639 template <
class Int,
class Entry>
640 int Basker<Int, Entry>::returnU(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
647 *col_ptr =
new Int[U->nrow+1];
649 *row_idx =
new Int[U->nnz];
651 *val =
new Entry[U->nnz];
653 if( (*col_ptr == NULL) || (*row_idx == NULL) || (*val == NULL) )
658 for(i = 0; i < U->nrow+1; i++)
660 (*col_ptr)[i] = U->col_ptr[i];
662 for(i = 0; i < U->nnz; i++)
664 (*row_idx)[i] = U->row_idx[i];
665 (*val)[i] = U->val[i];
670 template <
class Int,
class Entry>
671 int Basker<Int, Entry>::returnP(Int** p)
675 *p =
new Int[A->nrow];
682 for(i = 0; i < A->nrow; i++)
689 template <
class Int,
class Entry>
690 void Basker<Int, Entry>::free_factor()
709 template <
class Int,
class Entry>
710 void Basker<Int, Entry>::free_perm_matrix()
717 template <
class Int,
class Entry>
718 int Basker<Int, Entry>::solveMultiple(Int nrhs, Entry *b, Entry *x)
721 for(i = 0; i < nrhs; i++)
724 int result = solve(&(b[k]), &(x[k]));
727 cout <<
"Error in Solving \n";
735 template <
class Int,
class Entry>
736 int Basker<Int, Entry>::solve(Entry* b, Entry* x)
744 Entry* temp =
new Entry[A->nrow]();
747 for(i = 0 ; i < A->ncol; i++)
753 result = low_tri_solve_csc(L->nrow, L->col_ptr, L->row_idx, L->val, temp, x);
756 result = up_tri_solve_csc(U->nrow, U->col_ptr, U->row_idx, U->val, x, temp);
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)
770 for(i = 0; i < n ; i++)
773 ASSERT(val[col_ptr[i]] != (Entry)0);
775 if(val[col_ptr[i]] == (Entry) 0)
780 x[i] = BASKER_ScalarTraits<Entry>::divide(b[i], val[col_ptr[i]]);
782 for(j = col_ptr[i]+1; j < (col_ptr[i+1]); j++)
784 b[pinv[row_idx[j]]] = b[pinv[row_idx[j]]] - (val[j]*x[i]);
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)
795 for(i = n; i > 1 ; i--)
799 ASSERT(val[col_ptr[i]-1] != (Entry)0);
801 if(val[col_ptr[i]-1] == (Entry) 0)
803 cout <<
"Dig(" << i <<
") = " << val[col_ptr[i]-1] << endl;
808 x[ii] = BASKER_ScalarTraits<Entry>::divide(b[ii],val[col_ptr[i]-1]);
810 for(j = (col_ptr[i]-2); j >= (col_ptr[ii]); j--)
812 b[row_idx[j]] = b[row_idx[j]] - (val[j]*x[ii]);
816 x[0] = BASKER_ScalarTraits<Entry>::divide(b[0],val[col_ptr[1]-1]);
820 template <
class Int,
class Entry>
821 int Basker<Int, Entry>::preorder(Int *row_perm, Int *col_perm)
824 basker_matrix <Int, Entry> *B;
825 B =
new basker_matrix<Int, Entry>;
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));
833 if( (B->col_ptr == NULL) || (B->row_idx == NULL) || (B->val == NULL) )
839 (void) permute_column(col_perm, B);
840 (void) permute_row(row_perm, B);
844 A->col_ptr = B->col_ptr;
845 A->row_idx = B->row_idx;
853 template <
class Int,
class Entry>
854 int Basker <Int, Entry>::permute_column(Int *p, basker_matrix<Int,Entry> *B)
860 for(j=0; j < B->ncol; j++)
863 B->col_ptr[i+1] = A->col_ptr[j+1] - A->col_ptr[j];
867 for(j=0; j < B->ncol; j++)
869 B->col_ptr[j+1] = B->col_ptr[j+1] + B->col_ptr[j];
874 for(ii = 0 ; ii < B->ncol; ii++)
876 ko = B->col_ptr(p[ii]);
877 for(k = A->col_ptr[ii]; k < A->col_ptr[ii+1]; k++)
879 B->row_index[ko] = A->row_index[k];
880 B->val[ko] = A->val[ko];
887 template <
class Int,
class Entry>
888 int Basker <Int, Entry>::permute_row(Int *p, basker_matrix<Int,Entry> *B)
891 for(k=0; k < A->nnz; k++)
893 B->row_idx[k] = p[A->row_idx[k]];
898 template <
class Int,
class Entry>
899 int Basker <Int, Entry>::sort_factors()
906 for(i = 0 ; i < L->ncol; i++)
911 for(j = L->col_ptr[i]+1; j < (L->col_ptr[i+1]); j++)
913 if(L->row_idx[j] < val)
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;
929 for(i = 0 ; i < U->ncol; i++)
931 p = U->col_ptr[i+1]-1;
934 for(j = U->col_ptr[i]; j < (U->col_ptr[i+1]-1); j++)
936 if(U->row_idx[j] > val)
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;
953 template <
class Int,
class Entry>
954 Entry* Basker <Int, Entry>::entry_realloc(Entry *old , Int old_size, Int new_size)
956 Entry *new_entry =
new Entry[new_size];
957 for(Int i = 0; i < old_size; i++)
960 new_entry[i] = old[i];
967 template <
class Int,
class Entry>
968 Int* Basker <Int, Entry>::int_realloc(Int *old, Int old_size, Int new_size)
970 Int *new_int =
new Int[new_size];
971 for(Int i =0; i < old_size; i++)
Definition: basker.cpp:35