Amesos2 - Direct Sparse Solver Interfaces Version of the Day
klu2_analyze_given.hpp
1/* ========================================================================== */
2/* === klu_analyze_given ==================================================== */
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/* Given an input permutation P and Q, create the Symbolic object. BTF can
38 * be done to modify the user's P and Q (does not perform the max transversal;
39 * just finds the strongly-connected components). */
40
41#ifndef KLU2_ANALYZE_GIVEN_HPP
42#define KLU2_ANALYZE_GIVEN_HPP
43
44#include "klu2_internal.h"
45#include "klu2_memory.hpp"
46
47/* ========================================================================== */
48/* === klu_alloc_symbolic =================================================== */
49/* ========================================================================== */
50
51/* Allocate Symbolic object, and check input matrix. Not user callable. */
52
53template <typename Entry, typename Int>
54KLU_symbolic<Entry, Int> *KLU_alloc_symbolic
55(
56 Int n,
57 Int *Ap,
58 Int *Ai,
59 KLU_common<Entry, Int> *Common
60)
61{
62 KLU_symbolic<Entry, Int> *Symbolic ;
63 Int *P, *Q, *R ;
64 double *Lnz ;
65 Int nz, i, j, p, pend ;
66
67 if (Common == NULL)
68 {
69 return (NULL) ;
70 }
71 Common->status = KLU_OK ;
72
73 /* A is n-by-n, with n > 0. Ap [0] = 0 and nz = Ap [n] >= 0 required.
74 * Ap [j] <= Ap [j+1] must hold for all j = 0 to n-1. Row indices in Ai
75 * must be in the range 0 to n-1, and no duplicate entries can be present.
76 * The list of row indices in each column of A need not be sorted.
77 */
78
79 if (n <= 0 || Ap == NULL || Ai == NULL)
80 {
81 /* Ap and Ai must be present, and n must be > 0 */
82 Common->status = KLU_INVALID ;
83 return (NULL) ;
84 }
85
86 nz = Ap [n] ;
87 if (Ap [0] != 0 || nz < 0)
88 {
89 /* nz must be >= 0 and Ap [0] must equal zero */
90 Common->status = KLU_INVALID ;
91 return (NULL) ;
92 }
93
94 for (j = 0 ; j < n ; j++)
95 {
96 if (Ap [j] > Ap [j+1])
97 {
98 /* column pointers must be non-decreasing */
99 Common->status = KLU_INVALID ;
100 return (NULL) ;
101 }
102 }
103 P = (Int *) KLU_malloc (n, sizeof (Int), Common) ;
104 if (Common->status < KLU_OK)
105 {
106 /* out of memory */
107 Common->status = KLU_OUT_OF_MEMORY ;
108 return (NULL) ;
109 }
110 for (i = 0 ; i < n ; i++)
111 {
112 P [i] = AMESOS2_KLU2_EMPTY ;
113 }
114 for (j = 0 ; j < n ; j++)
115 {
116 pend = Ap [j+1] ;
117 for (p = Ap [j] ; p < pend ; p++)
118 {
119 i = Ai [p] ;
120 if (i < 0 || i >= n || P [i] == j)
121 {
122 /* row index out of range, or duplicate entry */
123 KLU_free (P, n, sizeof (Int), Common) ;
124 Common->status = KLU_INVALID ;
125 return (NULL) ;
126 }
127 /* flag row i as appearing in column j */
128 P [i] = j ;
129 }
130 }
131
132 /* ---------------------------------------------------------------------- */
133 /* allocate the Symbolic object */
134 /* ---------------------------------------------------------------------- */
135
136 Symbolic = (KLU_symbolic<Entry, Int> *) KLU_malloc (sizeof (KLU_symbolic<Entry, Int>), 1, Common) ;
137 if (Common->status < KLU_OK)
138 {
139 /* out of memory */
140 KLU_free (P, n, sizeof (Int), Common) ;
141 Common->status = KLU_OUT_OF_MEMORY ;
142 return (NULL) ;
143 }
144
145 Q = (Int *) KLU_malloc (n, sizeof (Int), Common) ;
146 R = (Int *) KLU_malloc (n+1, sizeof (Int), Common) ;
147 Lnz = (double *) KLU_malloc (n, sizeof (double), Common) ;
148
149 Symbolic->n = n ;
150 Symbolic->nz = nz ;
151 Symbolic->P = P ;
152 Symbolic->Q = Q ;
153 Symbolic->R = R ;
154 Symbolic->Lnz = Lnz ;
155
156 if (Common->status < KLU_OK)
157 {
158 /* out of memory */
159 KLU_free_symbolic (&Symbolic, Common) ;
160 Common->status = KLU_OUT_OF_MEMORY ;
161 return (NULL) ;
162 }
163
164 return (Symbolic) ;
165}
166
167
168/* ========================================================================== */
169/* === KLU_analyze_given ==================================================== */
170/* ========================================================================== */
171
172template <typename Entry, typename Int>
173KLU_symbolic<Entry, Int> *KLU_analyze_given /* returns NULL if error, or a valid
174 KLU_symbolic object if successful */
175(
176 /* inputs, not modified */
177 Int n, /* A is n-by-n */
178 Int Ap [ ], /* size n+1, column pointers */
179 Int Ai [ ], /* size nz, row indices */
180 Int Puser [ ], /* size n, user's row permutation (may be NULL) */
181 Int Quser [ ], /* size n, user's column permutation (may be NULL) */
182 /* -------------------- */
183 KLU_common<Entry, Int> *Common
184)
185{
186 KLU_symbolic<Entry, Int> *Symbolic ;
187 double *Lnz ;
188 Int nblocks, nz, block, maxblock, *P, *Q, *R, nzoff, p, pend, do_btf, k ;
189
190 /* ---------------------------------------------------------------------- */
191 /* determine if input matrix is valid, and get # of nonzeros */
192 /* ---------------------------------------------------------------------- */
193
194 Symbolic = KLU_alloc_symbolic (n, Ap, Ai, Common) ;
195 if (Symbolic == NULL)
196 {
197 return (NULL) ;
198 }
199 P = Symbolic->P ;
200 Q = Symbolic->Q ;
201 R = Symbolic->R ;
202 Lnz = Symbolic->Lnz ;
203 nz = Symbolic->nz ;
204
205 /* ---------------------------------------------------------------------- */
206 /* Q = Quser, or identity if Quser is NULL */
207 /* ---------------------------------------------------------------------- */
208
209 if (Quser == (Int *) NULL)
210 {
211 for (k = 0 ; k < n ; k++)
212 {
213 Q [k] = k ;
214 }
215 }
216 else
217 {
218 for (k = 0 ; k < n ; k++)
219 {
220 Q [k] = Quser [k] ;
221 }
222 }
223
224 /* ---------------------------------------------------------------------- */
225 /* get the control parameters for BTF and ordering method */
226 /* ---------------------------------------------------------------------- */
227
228 do_btf = Common->btf ;
229 do_btf = (do_btf) ? TRUE : FALSE ;
230 Symbolic->ordering = 2 ;
231 Symbolic->do_btf = do_btf ;
232
233 /* ---------------------------------------------------------------------- */
234 /* find the block triangular form, if requested */
235 /* ---------------------------------------------------------------------- */
236
237 if (do_btf)
238 {
239
240 /* ------------------------------------------------------------------ */
241 /* get workspace for BTF_strongcomp */
242 /* ------------------------------------------------------------------ */
243
244 Int *Pinv, *Work, *Bi, k1, k2, nk, oldcol ;
245
246 Work = (Int *) KLU_malloc (4*n, sizeof (Int), Common) ;
247 Pinv = (Int *) KLU_malloc (n, sizeof (Int), Common) ;
248 if (Puser != (Int *) NULL)
249 {
250 Bi = (Int *) KLU_malloc (nz+1, sizeof (Int), Common) ;
251 }
252 else
253 {
254 Bi = Ai ;
255 }
256
257 if (Common->status < KLU_OK)
258 {
259 /* out of memory */
260 KLU_free (Work, 4*n, sizeof (Int), Common) ;
261 KLU_free (Pinv, n, sizeof (Int), Common) ;
262 if (Puser != (Int *) NULL)
263 {
264 KLU_free (Bi, nz+1, sizeof (Int), Common) ;
265 }
266 KLU_free_symbolic (&Symbolic, Common) ;
267 Common->status = KLU_OUT_OF_MEMORY ;
268 return (NULL) ;
269 }
270
271 /* ------------------------------------------------------------------ */
272 /* B = Puser * A */
273 /* ------------------------------------------------------------------ */
274
275 if (Puser != (Int *) NULL)
276 {
277 for (k = 0 ; k < n ; k++)
278 {
279 Pinv [Puser [k]] = k ;
280 }
281 for (p = 0 ; p < nz ; p++)
282 {
283 Bi [p] = Pinv [Ai [p]] ;
284 }
285 }
286
287 /* ------------------------------------------------------------------ */
288 /* find the strongly-connected components */
289 /* ------------------------------------------------------------------ */
290
291 /* TODO : Correct version of BTF */
292 /* modifies Q, and determines P and R */
293 /*nblocks = BTF_strongcomp (n, Ap, Bi, Q, P, R, Work) ;*/
294 nblocks = KLU_OrdinalTraits<Int>::btf_strongcomp (n, Ap, Bi, Q, P, R,
295 Work) ;
296
297 /* ------------------------------------------------------------------ */
298 /* P = P * Puser */
299 /* ------------------------------------------------------------------ */
300
301 if (Puser != (Int *) NULL)
302 {
303 for (k = 0 ; k < n ; k++)
304 {
305 Work [k] = Puser [P [k]] ;
306 }
307 for (k = 0 ; k < n ; k++)
308 {
309 P [k] = Work [k] ;
310 }
311 }
312
313 /* ------------------------------------------------------------------ */
314 /* Pinv = inverse of P */
315 /* ------------------------------------------------------------------ */
316
317 for (k = 0 ; k < n ; k++)
318 {
319 Pinv [P [k]] = k ;
320 }
321
322 /* ------------------------------------------------------------------ */
323 /* analyze each block */
324 /* ------------------------------------------------------------------ */
325
326 nzoff = 0 ; /* nz in off-diagonal part */
327 maxblock = 1 ; /* size of the largest block */
328
329 for (block = 0 ; block < nblocks ; block++)
330 {
331
332 /* -------------------------------------------------------------- */
333 /* the block is from rows/columns k1 to k2-1 */
334 /* -------------------------------------------------------------- */
335
336 k1 = R [block] ;
337 k2 = R [block+1] ;
338 nk = k2 - k1 ;
339 PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk)) ;
340 maxblock = MAX (maxblock, nk) ;
341
342 /* -------------------------------------------------------------- */
343 /* scan the kth block, C */
344 /* -------------------------------------------------------------- */
345
346 for (k = k1 ; k < k2 ; k++)
347 {
348 oldcol = Q [k] ;
349 pend = Ap [oldcol+1] ;
350 for (p = Ap [oldcol] ; p < pend ; p++)
351 {
352 if (Pinv [Ai [p]] < k1)
353 {
354 nzoff++ ;
355 }
356 }
357 }
358
359 /* fill-in not estimated */
360 Lnz [block] = AMESOS2_KLU2_EMPTY ;
361 }
362
363 /* ------------------------------------------------------------------ */
364 /* free all workspace */
365 /* ------------------------------------------------------------------ */
366
367 KLU_free (Work, 4*n, sizeof (Int), Common) ;
368 KLU_free (Pinv, n, sizeof (Int), Common) ;
369 if (Puser != (Int *) NULL)
370 {
371 KLU_free (Bi, nz+1, sizeof (Int), Common) ;
372 }
373
374 }
375 else
376 {
377
378 /* ------------------------------------------------------------------ */
379 /* BTF not requested */
380 /* ------------------------------------------------------------------ */
381
382 nzoff = 0 ;
383 nblocks = 1 ;
384 maxblock = n ;
385 R [0] = 0 ;
386 R [1] = n ;
387 Lnz [0] = AMESOS2_KLU2_EMPTY ;
388
389 /* ------------------------------------------------------------------ */
390 /* P = Puser, or identity if Puser is NULL */
391 /* ------------------------------------------------------------------ */
392
393 for (k = 0 ; k < n ; k++)
394 {
395 P [k] = (Puser == NULL) ? k : Puser [k] ;
396 }
397 }
398
399 /* ---------------------------------------------------------------------- */
400 /* return the symbolic object */
401 /* ---------------------------------------------------------------------- */
402
403 Symbolic->nblocks = nblocks ;
404 Symbolic->maxblock = maxblock ;
405 Symbolic->lnz = AMESOS2_KLU2_EMPTY ;
406 Symbolic->unz = AMESOS2_KLU2_EMPTY ;
407 Symbolic->nzoff = nzoff ;
408
409 return (Symbolic) ;
410}
411
412#endif /* KLU2_ANALYZE_GIVEN_HPP */