Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
klu2_sort.hpp
1 /* ========================================================================== */
2 /* === KLU_sort ============================================================= */
3 /* ========================================================================== */
4 // @HEADER
5 // ***********************************************************************
6 //
7 // KLU2: A Direct Linear Solver package
8 // Copyright 2011 Sandia Corporation
9 //
10 // Under terms of Contract DE-AC04-94AL85000, with Sandia Corporation, the
11 // U.S. Government retains certain rights in this software.
12 //
13 // This library is free software; you can redistribute it and/or modify
14 // it under the terms of the GNU Lesser General Public License as
15 // published by the Free Software Foundation; either version 2.1 of the
16 // License, or (at your option) any later version.
17 //
18 // This library is distributed in the hope that it will be useful, but
19 // WITHOUT ANY WARRANTY; without even the implied warranty of
20 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21 // Lesser General Public License for more details.
22 //
23 // You should have received a copy of the GNU Lesser General Public
24 // License along with this library; if not, write to the Free Software
25 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
26 // USA
27 // Questions? Contact Mike A. Heroux (maherou@sandia.gov)
28 //
29 // KLU2 is derived work from KLU, licensed under LGPL, and copyrighted by
30 // University of Florida. The Authors of KLU are Timothy A. Davis and
31 // Eka Palamadai. See Doc/KLU_README.txt for the licensing and copyright
32 // information for KLU.
33 //
34 // ***********************************************************************
35 // @HEADER
36 
37 /* sorts the columns of L and U so that the row indices appear in strictly
38  * increasing order.
39  */
40 
41 #ifndef KLU2_SORT_HPP
42 #define KLU2_SORT_HPP
43 
44 #include "klu2_internal.h"
45 #include "klu2_memory.hpp"
46 
47 /* ========================================================================== */
48 /* === sort ================================================================= */
49 /* ========================================================================== */
50 
51 /* Sort L or U using a double-transpose */
52 
53 template <typename Entry, typename Int>
54 static void sort (Int n, Int *Xip, Int *Xlen, Unit *LU, Int *Tp, Int *Tj,
55  Entry *Tx, Int *W)
56 {
57  Int *Xi ;
58  Entry *Xx ;
59  Int p, i, j, len, nz, tp, xlen, pend ;
60 
61  ASSERT (KLU_valid_LU (n, FALSE, Xip, Xlen, LU)) ;
62 
63  /* count the number of entries in each row of L or U */
64  for (i = 0 ; i < n ; i++)
65  {
66  W [i] = 0 ;
67  }
68  for (j = 0 ; j < n ; j++)
69  {
70  GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len) ;
71  for (p = 0 ; p < len ; p++)
72  {
73  W [Xi [p]]++ ;
74  }
75  }
76 
77  /* construct the row pointers for T */
78  nz = 0 ;
79  for (i = 0 ; i < n ; i++)
80  {
81  Tp [i] = nz ;
82  nz += W [i] ;
83  }
84  Tp [n] = nz ;
85  for (i = 0 ; i < n ; i++)
86  {
87  W [i] = Tp [i] ;
88  }
89 
90  /* transpose the matrix into Tp, Ti, Tx */
91  for (j = 0 ; j < n ; j++)
92  {
93  GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len) ;
94  for (p = 0 ; p < len ; p++)
95  {
96  tp = W [Xi [p]]++ ;
97  Tj [tp] = j ;
98  Tx [tp] = Xx [p] ;
99  }
100  }
101 
102  /* transpose the matrix back into Xip, Xlen, Xi, Xx */
103  for (j = 0 ; j < n ; j++)
104  {
105  W [j] = 0 ;
106  }
107  for (i = 0 ; i < n ; i++)
108  {
109  pend = Tp [i+1] ;
110  for (p = Tp [i] ; p < pend ; p++)
111  {
112  j = Tj [p] ;
113  GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len) ;
114  xlen = W [j]++ ;
115  Xi [xlen] = i ;
116  Xx [xlen] = Tx [p] ;
117  }
118  }
119 
120  ASSERT (KLU_valid_LU (n, FALSE, Xip, Xlen, LU)) ;
121 }
122 
123 
124 /* ========================================================================== */
125 /* === KLU_sort ============================================================= */
126 /* ========================================================================== */
127 
128 template <typename Entry, typename Int>
129 Int KLU_sort
130 (
131  KLU_symbolic<Entry, Int> *Symbolic,
132  KLU_numeric<Entry, Int> *Numeric,
133  KLU_common<Entry, Int> *Common
134 )
135 {
136  Int *R, *W, *Tp, *Ti, *Lip, *Uip, *Llen, *Ulen ;
137  Entry *Tx ;
138  Unit **LUbx ;
139  Int n, nk, nz, block, nblocks, maxblock, k1 ;
140  size_t m1 ;
141 
142  if (Common == NULL)
143  {
144  return (FALSE) ;
145  }
146  Common->status = KLU_OK ;
147 
148  n = Symbolic->n ;
149  R = Symbolic->R ;
150  nblocks = Symbolic->nblocks ;
151  maxblock = Symbolic->maxblock ;
152 
153  Lip = Numeric->Lip ;
154  Llen = Numeric->Llen ;
155  Uip = Numeric->Uip ;
156  Ulen = Numeric->Ulen ;
157  LUbx = (Unit **) Numeric->LUbx ;
158 
159  m1 = ((size_t) maxblock) + 1 ;
160 
161  /* allocate workspace */
162  nz = MAX (Numeric->max_lnz_block, Numeric->max_unz_block) ;
163  W = (Int *) KLU_malloc (maxblock, sizeof (Int), Common) ;
164  Tp = (Int *) KLU_malloc (m1, sizeof (Int), Common) ;
165  Ti = (Int *) KLU_malloc (nz, sizeof (Int), Common) ;
166  Tx = (Entry *) KLU_malloc (nz, sizeof (Entry), Common) ;
167 
168  PRINTF (("\n======================= Start sort:\n")) ;
169 
170  if (Common->status == KLU_OK)
171  {
172  /* sort each block of L and U */
173  for (block = 0 ; block < nblocks ; block++)
174  {
175  k1 = R [block] ;
176  nk = R [block+1] - k1 ;
177  if (nk > 1)
178  {
179  PRINTF (("\n-------------------block: %d nk %d\n", block, nk)) ;
180  sort (nk, Lip + k1, Llen + k1, LUbx [block], Tp, Ti, Tx, W) ;
181  sort (nk, Uip + k1, Ulen + k1, LUbx [block], Tp, Ti, Tx, W) ;
182  }
183  }
184  }
185 
186  PRINTF (("\n======================= sort done.\n")) ;
187 
188  /* free workspace */
189  KLU_free (W, maxblock, sizeof (Int), Common) ;
190  KLU_free (Tp, m1, sizeof (Int), Common) ;
191  KLU_free (Ti, nz, sizeof (Int), Common) ;
192  KLU_free (Tx, nz, sizeof (Entry), Common) ;
193  return (Common->status == KLU_OK) ;
194 }
195 
196 #endif