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