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()
110 template <
class Int,
class Entry>
111 int Basker<Int,Entry>:: basker_dfs
127 Int start, end, done, *store ;
132 bool has_elements =
true;
149 ASSERT (color[j] == 1) ;
157 for ( i1 = start ; i1 < end ; i1++ )
171 pattern[--*top] = j ;
175 has_elements =
false;
184 std::cout <<
"Out of DFS: " << j << std::endl;
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)
205 A->col_ptr = col_ptr;
206 A->row_idx = row_idx;
218 L->col_ptr =
new Int[ncol+1]();
220 L->row_idx =
new Int[L->nnz]();
222 L->val =
new Entry[L->nnz]();
231 U->col_ptr =
new Int[ncol+1]();
233 U->row_idx =
new Int[U->nnz]();
235 U->val =
new Entry[U->nnz]();
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))
249 tptr =
new Int[(ncol)+(4*nrow)]();
251 X =
new Entry[2*nrow]();
253 pinv =
new Int[ncol+1]();
256 if( (tptr == NULL) || (X == NULL) || (pinv == NULL) )
266 Int *color, *pattern, *stack;
267 Int top, top1, maxindex, t;
268 Int lnnz, unnz, xnnz, lcnt, ucnt;
269 Int cu_ltop, cu_utop;
272 Entry pivot, value, xj;
291 for(k = 0 ; k < ncol; k++)
297 for (k = 0; k < ncol; k++)
301 cout <<
"k = " << k << endl;
313 ASSERT (top == ncol);
315 for(i = 0; i < nrow; i++)
317 ASSERT(X[i] == (Entry)0);
319 for(i = 0; i < ncol; i++)
321 ASSERT(color[i] == 0);
325 for( i = col_ptr[k]; i < col_ptr[k+1]; i++)
333 basker_dfs(nrow, j, L->row_idx, L->col_ptr, color, pattern, &top, pinv, stack);
342 cout << ncol << endl;
343 cout << xnnz << endl;
348 for(pp = 0; pp < xnnz; pp++)
357 p2 = L->col_ptr[t+1];
358 for(p = L->col_ptr[t]+1; p < p2; p++)
360 X[L->row_idx[p]] -= L->val[p] * xj;
368 for(i = top; i < nrow; i++)
375 absv = BASKER_ScalarTraits<Entry>::approxABS(value);
380 if( BASKER_ScalarTraits<Entry>::gt(absv , maxv))
389 ucnt = nrow - top - lcnt + 1;
391 if(maxindex == ncol || pivot == ((Entry)0))
393 cout <<
"Matrix is singular at index: " << maxindex <<
" pivot: " << pivot << endl;
402 cout <<
"Permuting pivot: " << k <<
" for row: " << maxindex << endl;
406 if(lnnz + lcnt >= L->nnz)
409 newsize = L->nnz * 1.1 + 2*nrow + 1;
411 cout <<
"Out of memory -- Reallocating. Old Size: " << L->nnz <<
" New Size: " << newsize << endl;
414 L->row_idx = int_realloc(L->row_idx , L->nnz, newsize);
417 cout <<
"WARNING: Cannot Realloc Memory" << endl;
422 L->val = entry_realloc(L->val, L->nnz, newsize);
425 cout <<
"WARNING: Cannot Realloc Memory" << endl;
433 if(unnz + ucnt >= U->nnz)
435 newsize = U->nnz*1.1 + 2*nrow + 1;
437 cout <<
"Out of memory -- Reallocating. Old Size: " << U->nnz <<
" New Size: " << newsize << endl;
440 U->row_idx = int_realloc(U->row_idx, U->nnz, newsize);
443 cout <<
"WARNING: Cannot Realloc Memory" << endl;
449 U->val = entry_realloc(U->val, U->nnz, newsize);
452 cout <<
"WARNING: Cannot Realloc Memory" << endl;
460 L->row_idx[lnnz] = maxindex;
464 Entry last_v_temp = 0;
466 for(i = top; i < nrow; i++)
474 if(X[j] != ((Entry)0))
481 cout <<
"BASKER: Insufficent memory for U" << endl;
488 U->row_idx[unnz] = pinv[j];
504 cout <<
"BASKER: Insufficent memory for L" << endl;
509 L->row_idx[lnnz] = j;
511 L->val[lnnz] = BASKER_ScalarTraits<Entry>::divide(X[j],pivot);
523 U->row_idx[unnz] = k;
524 U->val[unnz] = last_v_temp;
530 L->col_ptr[k] = cu_ltop;
531 L->col_ptr[k+1] = lnnz;
534 U->col_ptr[k] = cu_utop;
535 U->col_ptr[k+1] = unnz;
542 for(k = 0; k < lnnz; k++)
544 printf(
"L[%d]=%g" , k , L->val[k]);
547 for(k = 0; k < lnnz; k++)
549 printf(
"Li[%d]=%d", k, L->row_idx[k]);
552 for(k = 0; k < nrow; k++)
554 printf(
"p[%d]=%d", k, pinv[k]);
559 for(k = 0; k < ncol; k++)
561 printf(
"Up[%d]=%d", k, U->col_ptr[k]);
565 for(k = 0; k < unnz; k++)
567 printf(
"U[%d]=%g" , k , U->val[k]);
570 for(k = 0; k < unnz; k++)
572 printf(
"Ui[%d]=%d", k, U->row_idx[k]);
579 for( i = 0; i < ncol; i++)
581 for(k = L->col_ptr[i]; k < L->col_ptr[i+1]; k++)
591 cout <<
"After Permuting" << endl;
592 for(k = 0; k < lnnz; k++)
594 printf(
"Li[%d]=%d", k, L->row_idx[k]);
605 template <
class Int,
class Entry>
606 int Basker<Int, Entry>::returnL(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
615 *col_ptr =
new Int[L->nrow+1];
617 *row_idx =
new Int[L->nnz];
619 *val =
new Entry[L->nnz];
621 if( (*col_ptr == NULL) || (*row_idx == NULL) || (*val == NULL) )
626 for(i = 0; i < L->nrow+1; i++)
628 (*col_ptr)[i] = L->col_ptr[i];
631 for(i = 0; i < L->nnz; i++)
633 (*row_idx)[i] = pinv[L->row_idx[i]];
634 (*val)[i] = L->val[i];
640 template <
class Int,
class Entry>
641 int Basker<Int, Entry>::returnU(Int *dim, Int *nnz, Int **col_ptr, Int **row_idx, Entry **val)
648 *col_ptr =
new Int[U->nrow+1];
650 *row_idx =
new Int[U->nnz];
652 *val =
new Entry[U->nnz];
654 if( (*col_ptr == NULL) || (*row_idx == NULL) || (*val == NULL) )
659 for(i = 0; i < U->nrow+1; i++)
661 (*col_ptr)[i] = U->col_ptr[i];
663 for(i = 0; i < U->nnz; i++)
665 (*row_idx)[i] = U->row_idx[i];
666 (*val)[i] = U->val[i];
671 template <
class Int,
class Entry>
672 int Basker<Int, Entry>::returnP(Int** p)
676 *p =
new Int[A->nrow];
683 for(i = 0; i < A->nrow; i++)
690 template <
class Int,
class Entry>
691 void Basker<Int, Entry>::free_factor()
710 template <
class Int,
class Entry>
711 void Basker<Int, Entry>::free_perm_matrix()
718 template <
class Int,
class Entry>
719 int Basker<Int, Entry>::solveMultiple(Int nrhs, Entry *b, Entry *x)
722 for(i = 0; i < nrhs; i++)
725 int result = solve(&(b[k]), &(x[k]));
728 cout <<
"Error in Solving \n";
736 template <
class Int,
class Entry>
737 int Basker<Int, Entry>::solve(Entry* b, Entry* x)
745 Entry* temp =
new Entry[A->nrow]();
748 for(i = 0 ; i < A->ncol; i++)
754 result = low_tri_solve_csc(L->nrow, L->col_ptr, L->row_idx, L->val, temp, x);
757 result = up_tri_solve_csc(U->nrow, U->col_ptr, U->row_idx, U->val, x, temp);
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)
771 for(i = 0; i < n ; i++)
774 ASSERT(val[col_ptr[i]] != (Entry)0);
776 if(val[col_ptr[i]] == (Entry) 0)
781 x[i] = BASKER_ScalarTraits<Entry>::divide(b[i], val[col_ptr[i]]);
783 for(j = col_ptr[i]+1; j < (col_ptr[i+1]); j++)
785 b[pinv[row_idx[j]]] = b[pinv[row_idx[j]]] - (val[j]*x[i]);
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)
796 for(i = n; i > 1 ; i--)
800 ASSERT(val[col_ptr[i]-1] != (Entry)0);
802 if(val[col_ptr[i]-1] == (Entry) 0)
804 cout <<
"Dig(" << i <<
") = " << val[col_ptr[i]-1] << endl;
809 x[ii] = BASKER_ScalarTraits<Entry>::divide(b[ii],val[col_ptr[i]-1]);
811 for(j = (col_ptr[i]-2); j >= (col_ptr[ii]); j--)
813 b[row_idx[j]] = b[row_idx[j]] - (val[j]*x[ii]);
817 x[0] = BASKER_ScalarTraits<Entry>::divide(b[0],val[col_ptr[1]-1]);
821 template <
class Int,
class Entry>
822 int Basker<Int, Entry>::preorder(Int *row_perm, Int *col_perm)
825 basker_matrix <Int, Entry> *B;
826 B =
new basker_matrix<Int, Entry>;
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));
834 if( (B->col_ptr == NULL) || (B->row_idx == NULL) || (B->val == NULL) )
840 (void) permute_column(col_perm, B);
841 (void) permute_row(row_perm, B);
845 A->col_ptr = B->col_ptr;
846 A->row_idx = B->row_idx;
854 template <
class Int,
class Entry>
855 int Basker <Int, Entry>::permute_column(Int *p, basker_matrix<Int,Entry> *B)
861 for(j=0; j < B->ncol; j++)
864 B->col_ptr[i+1] = A->col_ptr[j+1] - A->col_ptr[j];
868 for(j=0; j < B->ncol; j++)
870 B->col_ptr[j+1] = B->col_ptr[j+1] + B->col_ptr[j];
875 for(ii = 0 ; ii < B->ncol; ii++)
877 ko = B->col_ptr(p[ii]);
878 for(k = A->col_ptr[ii]; k < A->col_ptr[ii+1]; k++)
880 B->row_index[ko] = A->row_index[k];
881 B->val[ko] = A->val[ko];
888 template <
class Int,
class Entry>
889 int Basker <Int, Entry>::permute_row(Int *p, basker_matrix<Int,Entry> *B)
892 for(k=0; k < A->nnz; k++)
894 B->row_idx[k] = p[A->row_idx[k]];
899 template <
class Int,
class Entry>
900 int Basker <Int, Entry>::sort_factors()
907 for(i = 0 ; i < L->ncol; i++)
912 for(j = L->col_ptr[i]+1; j < (L->col_ptr[i+1]); j++)
914 if(L->row_idx[j] < val)
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;
930 for(i = 0 ; i < U->ncol; i++)
932 p = U->col_ptr[i+1]-1;
935 for(j = U->col_ptr[i]; j < (U->col_ptr[i+1]-1); j++)
937 if(U->row_idx[j] > val)
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;
954 template <
class Int,
class Entry>
955 Entry* Basker <Int, Entry>::entry_realloc(Entry *old , Int old_size, Int new_size)
957 Entry *new_entry =
new Entry[new_size];
958 for(Int i = 0; i < old_size; i++)
961 new_entry[i] = old[i];
968 template <
class Int,
class Entry>
969 Int* Basker <Int, Entry>::int_realloc(Int *old, Int old_size, Int new_size)
971 Int *new_int =
new Int[new_size];
972 for(Int i =0; i < old_size; i++)
Definition: basker.cpp:35