?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Adding remaining functions.
[JKLU.git] / src / main / java / edu / ufl / cise / klu / tdouble / Dklu_analyze_given.java
blob334866686aa8c7b5a12d8d305a8fee07b3899a79
1 /**
2 * KLU: a sparse LU factorization algorithm.
3 * Copyright (C) 2004-2009, Timothy A. Davis.
4 * Copyright (C) 2011, Richard W. Lincoln.
5 * http://www.cise.ufl.edu/research/sparse/klu
6 *
7 * -------------------------------------------------------------------------
8 *
9 * KLU is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 * KLU is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with this Module; if not, write to the Free Software
21 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
22 *
23 */
24
25 package edu.ufl.cise.klu.tdouble;
26
27 import edu.ufl.cise.klu.common.KLU_common;
28 import edu.ufl.cise.klu.common.KLU_symbolic;
29
30 /**
31 * Given an input permutation P and Q, create the Symbolic object. BTF can
32 * be done to modify the user's P and Q (does not perform the max transversal;
33 * just finds the strongly-connected components).
34 */
35 public class Dklu_analyze_given extends Dklu_internal
36 {
37
38 /**
39 * Allocate Symbolic object, and check input matrix. Not user callable.
40 *
41 * @param n
42 * @param Ap
43 * @param Ai
44 * @param Common
45 * @return
46 */
47 protected static KLU_symbolic klu_alloc_symbolic(int n, int[] Ap, int[] Ai,
48 KLU_common Common)
49 {
50 KLU_symbolic Symbolic ;
51 int[] P, Q, R ;
52 double[] Lnz ;
53 int nz, i, j, p, pend ;
54
55 if (Common == null)
56 {
57 return (null) ;
58 }
59 Common.status = KLU_common.KLU_OK ;
60
61 /* A is n-by-n, with n > 0. Ap [0] = 0 and nz = Ap [n] >= 0 required.
62 * Ap [j] <= Ap [j+1] must hold for all j = 0 to n-1. Row indices in Ai
63 * must be in the range 0 to n-1, and no duplicate entries can be present.
64 * The list of row indices in each column of A need not be sorted.
65 */
66
67 if (n <= 0 || Ap == null || Ai == null)
68 {
69 /* Ap and Ai must be present, and n must be > 0 */
70 Common.status = KLU_common.KLU_INVALID ;
71 return (null) ;
72 }
73
74 nz = Ap [n] ;
75 if (Ap [0] != 0 || nz < 0)
76 {
77 /* nz must be >= 0 and Ap [0] must equal zero */
78 Common.status = KLU_common.KLU_INVALID ;
79 return (null) ;
80 }
81
82 for (j = 0 ; j < n ; j++)
83 {
84 if (Ap [j] > Ap [j+1])
85 {
86 /* column pointers must be non-decreasing */
87 Common.status = KLU_common.KLU_INVALID ;
88 return (null) ;
89 }
90 }
91 P = Dklu_memory.klu_malloc (n, sizeof (Int), Common) ;
92 if (Common.status < KLU_common.KLU_OK)
93 {
94 /* out of memory */
95 Common.status = KLU_common.KLU_OUT_OF_MEMORY ;
96 return (null) ;
97 }
98 for (i = 0 ; i < n ; i++)
99 {
100 P [i] = EMPTY ;
101 }
102 for (j = 0 ; j < n ; j++)
103 {
104 pend = Ap [j+1] ;
105 for (p = Ap [j] ; p < pend ; p++)
106 {
107 i = Ai [p] ;
108 if (i < 0 || i >= n || P [i] == j)
109 {
110 /* row index out of range, or duplicate entry */
111 Dklu_memory.klu_free (P, n, sizeof (Int), Common) ;
112 Common.status = KLU_common.KLU_INVALID ;
113 return (null) ;
114 }
115 /* flag row i as appearing in column j */
116 P [i] = j ;
117 }
118 }
119
120 /* ---------------------------------------------------------------------- */
121 /* allocate the Symbolic object */
122 /* ---------------------------------------------------------------------- */
123
124 Symbolic = Dklu_memory.klu_malloc (sizeof (KLU_symbolic), 1, Common) ;
125 if (Common.status < KLU_common.KLU_OK)
126 {
127 /* out of memory */
128 Dklu_memory.klu_free (P, n, sizeof (Int), Common) ;
129 Common.status = KLU_common.KLU_OUT_OF_MEMORY ;
130 return (null) ;
131 }
132
133 Q = Dklu_memory.klu_malloc (n, sizeof (Int), Common) ;
134 R = Dklu_memory.klu_malloc (n+1, sizeof (Int), Common) ;
135 Lnz = Dklu_memory.klu_malloc (n, sizeof (double), Common) ;
136
137 Symbolic.n = n ;
138 Symbolic.nz = nz ;
139 Symbolic.P = P ;
140 Symbolic.Q = Q ;
141 Symbolic.R = R ;
142 Symbolic.Lnz = Lnz ;
143
144 if (Common.status < KLU_common.KLU_OK)
145 {
146 /* out of memory */
147 Dklu_free_symbolic.klu_free_symbolic (Symbolic, Common) ;
148 Common.status = KLU_common.KLU_OUT_OF_MEMORY ;
149 return (null) ;
150 }
151
152 return (Symbolic) ;
153 }
154
155 /**
156 *
157 * @param n A is n-by-n
158 * @param Ap size n+1, column pointers
159 * @param Ai size nz, row indices
160 * @param Puser size n, user's row permutation (may be null)
161 * @param Quser size n, user's column permutation (may be null)
162 * @param Common
163 * @return
164 */
165 public static KLU_symbolic klu_analyze_given(int n, int[] Ap, int[] Ai,
166 int[] Puser, int[] Quser, KLU_common Common)
167 {
168 KLU_symbolic Symbolic ;
169 double[] Lnz ;
170 int nblocks, nz, block, maxblock, nzoff, p, pend, do_btf, k ;
171 int[] P, Q, R ;
172
173 /* ---------------------------------------------------------------------- */
174 /* determine if input matrix is valid, and get # of nonzeros */
175 /* ---------------------------------------------------------------------- */
176
177 Symbolic = Dklu_alloc_symbolic.klu_alloc_symbolic (n, Ap, Ai, Common) ;
178 if (Symbolic == null)
179 {
180 return (null) ;
181 }
182 P = Symbolic.P ;
183 Q = Symbolic.Q ;
184 R = Symbolic.R ;
185 Lnz = Symbolic.Lnz ;
186 nz = Symbolic.nz ;
187
188 /* ---------------------------------------------------------------------- */
189 /* Q = Quser, or identity if Quser is null */
190 /* ---------------------------------------------------------------------- */
191
192 if (Quser == null)
193 {
194 for (k = 0 ; k < n ; k++)
195 {
196 Q [k] = k ;
197 }
198 }
199 else
200 {
201 for (k = 0 ; k < n ; k++)
202 {
203 Q [k] = Quser [k] ;
204 }
205 }
206
207 /* ---------------------------------------------------------------------- */
208 /* get the control parameters for BTF and ordering method */
209 /* ---------------------------------------------------------------------- */
210
211 do_btf = Common.btf ;
212 do_btf = (do_btf) ? TRUE : FALSE ;
213 Symbolic.ordering = 2 ;
214 Symbolic.do_btf = do_btf ;
215
216 /* ---------------------------------------------------------------------- */
217 /* find the block triangular form, if requested */
218 /* ---------------------------------------------------------------------- */
219
220 if (do_btf == 1)
221 {
222
223 /* ------------------------------------------------------------------ */
224 /* get workspace for BTF_strongcomp */
225 /* ------------------------------------------------------------------ */
226
227 int[] Pinv, Work, Bi ;
228 int k1, k2, nk, oldcol ;
229
230 Work = Dklu_memory.klu_malloc (4*n, sizeof (Int), Common) ;
231 Pinv = Dklu_memory.klu_malloc (n, sizeof (Int), Common) ;
232 if (Puser != null)
233 {
234 Bi = Dklu_memory.klu_malloc (nz+1, sizeof (Int), Common) ;
235 }
236 else
237 {
238 Bi = Ai ;
239 }
240
241 if (Common.status < KLU_common.KLU_OK)
242 {
243 /* out of memory */
244 Dklu_memory.klu_free (Work, 4*n, sizeof (Int), Common) ;
245 Dklu_memory.klu_free (Pinv, n, sizeof (Int), Common) ;
246 if (Puser != null)
247 {
248 Dklu_memory.klu_free (Bi, nz+1, sizeof (Int), Common) ;
249 }
250 Dklu_free_symbolic.klu_free_symbolic (Symbolic, Common) ;
251 Common.status = KLU_common.KLU_OUT_OF_MEMORY ;
252 return (null) ;
253 }
254
255 /* ------------------------------------------------------------------ */
256 /* B = Puser * A */
257 /* ------------------------------------------------------------------ */
258
259 if (Puser != null)
260 {
261 for (k = 0 ; k < n ; k++)
262 {
263 Pinv [Puser [k]] = k ;
264 }
265 for (p = 0 ; p < nz ; p++)
266 {
267 Bi [p] = Pinv [Ai [p]] ;
268 }
269 }
270
271 /* ------------------------------------------------------------------ */
272 /* find the strongly-connected components */
273 /* ------------------------------------------------------------------ */
274
275 /* modifies Q, and determines P and R */
276 nblocks = Dbtf_strongcomp.btf_strongcomp (n, Ap, Bi, Q, P, R, Work) ;
277
278 /* ------------------------------------------------------------------ */
279 /* P = P * Puser */
280 /* ------------------------------------------------------------------ */
281
282 if (Puser != null)
283 {
284 for (k = 0 ; k < n ; k++)
285 {
286 Work [k] = Puser [P [k]] ;
287 }
288 for (k = 0 ; k < n ; k++)
289 {
290 P [k] = Work [k] ;
291 }
292 }
293
294 /* ------------------------------------------------------------------ */
295 /* Pinv = inverse of P */
296 /* ------------------------------------------------------------------ */
297
298 for (k = 0 ; k < n ; k++)
299 {
300 Pinv [P [k]] = k ;
301 }
302
303 /* ------------------------------------------------------------------ */
304 /* analyze each block */
305 /* ------------------------------------------------------------------ */
306
307 nzoff = 0 ; /* nz in off-diagonal part */
308 maxblock = 1 ; /* size of the largest block */
309
310 for (block = 0 ; block < nblocks ; block++)
311 {
312
313 /* -------------------------------------------------------------- */
314 /* the block is from rows/columns k1 to k2-1 */
315 /* -------------------------------------------------------------- */
316
317 k1 = R [block] ;
318 k2 = R [block+1] ;
319 nk = k2 - k1 ;
320 PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk)) ;
321 maxblock = MAX (maxblock, nk) ;
322
323 /* -------------------------------------------------------------- */
324 /* scan the kth block, C */
325 /* -------------------------------------------------------------- */
326
327 for (k = k1 ; k < k2 ; k++)
328 {
329 oldcol = Q [k] ;
330 pend = Ap [oldcol+1] ;
331 for (p = Ap [oldcol] ; p < pend ; p++)
332 {
333 if (Pinv [Ai [p]] < k1)
334 {
335 nzoff++ ;
336 }
337 }
338 }
339
340 /* fill-in not estimated */
341 Lnz [block] = EMPTY ;
342 }
343
344 /* ------------------------------------------------------------------ */
345 /* free all workspace */
346 /* ------------------------------------------------------------------ */
347
348 Dklu_memory.klu_free (Work, 4*n, sizeof (Int), Common) ;
349 Dklu_memory.klu_free (Pinv, n, sizeof (Int), Common) ;
350 if (Puser != null)
351 {
352 Dklu_memory.klu_free (Bi, nz+1, sizeof (Int), Common) ;
353 }
354
355 }
356 else
357 {
358
359 /* ------------------------------------------------------------------ */
360 /* BTF not requested */
361 /* ------------------------------------------------------------------ */
362
363 nzoff = 0 ;
364 nblocks = 1 ;
365 maxblock = n ;
366 R [0] = 0 ;
367 R [1] = n ;
368 Lnz [0] = EMPTY ;
369
370 /* ------------------------------------------------------------------ */
371 /* P = Puser, or identity if Puser is null */
372 /* ------------------------------------------------------------------ */
373
374 for (k = 0 ; k < n ; k++)
375 {
376 P [k] = (Puser == null) ? k : Puser [k] ;
377 }
378 }
379
380 /* ---------------------------------------------------------------------- */
381 /* return the symbolic object */
382 /* ---------------------------------------------------------------------- */
383
384 Symbolic.nblocks = nblocks ;
385 Symbolic.maxblock = maxblock ;
386 Symbolic.lnz = EMPTY ;
387 Symbolic.unz = EMPTY ;
388 Symbolic.nzoff = nzoff ;
389
390 return (Symbolic) ;
391 }
392
393 }
A sparse LU factorization algorithm in Java