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
53template <typename Entry, typename Int>
54static 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
128template <typename Entry, typename Int>
129Int 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