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