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