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 
29 template <typename Entry, typename Int>
30 KLU_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 
148 template <typename Entry, typename Int>
149 KLU_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 */
Ap
int Ap[]
Column offsets.
Definition: klu2_simple.cpp:52
Ai
int Ai[]
Row values.
Definition: klu2_simple.cpp:54
Generated by   doxygen 1.8.5