Amesos2 - Direct Sparse Solver Interfaces Version of the Day
klu2_ext.hpp
1/* ========================================================================== */
2/* === klu include file ===================================================== */
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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
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
38/* Include file for user programs that call klu_* routines */
39
40#ifndef _TKLU_H
41#define _TKLU_H
42
43
44#include "trilinos_amd.h"
45#include "trilinos_colamd.h"
46#include "trilinos_btf_decl.h"
47
48/* -------------------------------------------------------------------------- */
49/* Symbolic object - contains the pre-ordering computed by klu_analyze */
50/* -------------------------------------------------------------------------- */
51
52/* TODO : Entry is not needed in symbolic, numeric and common. Remove TODO */
53template <typename Entry, typename Int> struct klu_symbolic
54{
55 /* A (P,Q) is in upper block triangular form. The kth block goes from
56 * row/col index R [k] to R [k+1]-1. The estimated number of nonzeros
57 * in the L factor of the kth block is Lnz [k].
58 */
59
60 /* only computed if the AMD ordering is chosen: */
61 double symmetry ; /* symmetry of largest block */
62 double est_flops ; /* est. factorization flop count */
63 double lnz, unz ; /* estimated nz in L and U, including diagonals */
64 double *Lnz ; /* size n, but only Lnz [0..nblocks-1] is used */
65
66 /* computed for all orderings: */
67 Int
68 n, /* input matrix A is n-by-n */
69 nz, /* # entries in input matrix */
70 *P, /* size n */
71 *Q, /* size n */
72 *R, /* size n+1, but only R [0..nblocks] is used */
73 nzoff, /* nz in off-diagonal blocks */
74 nblocks, /* number of blocks */
75 maxblock, /* size of largest block */
76 ordering, /* ordering used (AMD, COLAMD, or GIVEN) */
77 do_btf ; /* whether or not BTF preordering was requested */
78
79 /* only computed if BTF preordering requested */
80 Int structural_rank ; /* 0 to n-1 if the matrix is structurally rank
81 * deficient. -1 if not computed. n if the matrix has
82 * full structural rank */
83
84} ;
85
86/* -------------------------------------------------------------------------- */
87/* Numeric object - contains the factors computed by klu_factor */
88/* -------------------------------------------------------------------------- */
89template <typename Entry, typename Int> struct klu_numeric
90{
91 /* LU factors of each block, the pivot row permutation, and the
92 * entries in the off-diagonal blocks */
93
94 Int n ; /* A is n-by-n */
95 Int nblocks ; /* number of diagonal blocks */
96 Int lnz ; /* actual nz in L, including diagonal */
97 Int unz ; /* actual nz in U, including diagonal */
98 Int max_lnz_block ; /* max actual nz in L in any one block, incl. diag */
99 Int max_unz_block ; /* max actual nz in U in any one block, incl. diag */
100 Int *Pnum ; /* size n. final pivot permutation */
101 Int *Pinv ; /* size n. inverse of final pivot permutation */
102
103 /* LU factors of each block */
104 Int *Lip ; /* size n. pointers into LUbx[block] for L */
105 Int *Uip ; /* size n. pointers into LUbx[block] for U */
106 Int *Llen ; /* size n. Llen [k] = # of entries in kth column of L */
107 Int *Ulen ; /* size n. Ulen [k] = # of entries in kth column of U */
108 void **LUbx ; /* L and U indices and entries (excl. diagonal of U) */
109 size_t *LUsize ; /* size of each LUbx [block], in sizeof (Unit) */
110 void *Udiag ; /* diagonal of U */
111
112 /* scale factors; can be NULL if no scaling */
113 double *Rs ; /* size n. Rs [i] is scale factor for row i */
114
115 /* permanent workspace for factorization and solve */
116 size_t worksize ; /* size (in bytes) of Work */
117 void *Work ; /* workspace */
118 void *Xwork ; /* alias into Numeric->Work */
119 Int *Iwork ; /* alias into Numeric->Work */
120
121 /* off-diagonal entries in a conventional compressed-column sparse matrix */
122 Int *Offp ; /* size n+1, column pointers */
123 Int *Offi ; /* size nzoff, row indices */
124 void *Offx ; /* size nzoff, numerical values */
125 Int nzoff ;
126
127} ;
128
129/* -------------------------------------------------------------------------- */
130/* KLU control parameters and statistics */
131/* -------------------------------------------------------------------------- */
132
133/* Common->status values */
134#define KLU_OK 0
135#define KLU_SINGULAR (1) /* status > 0 is a warning, not an error */
136#define KLU_OUT_OF_MEMORY (-2)
137#define KLU_INVALID (-3)
138#define KLU_TOO_LARGE (-4) /* integer overflow has occured */
139
140template <typename Entry, typename Int> struct klu_common
141{
142
143 /* --------------------------------------------------------------------- */
144 /* parameters */
145 /* --------------------------------------------------------------------- */
146
147 double tol ; /* pivot tolerance for diagonal preference */
148 double memgrow ; /* realloc memory growth size for LU factors */
149 double initmem_amd ; /* init. memory size with AMD: c*nnz(L) + n */
150 double initmem ; /* init. memory size: c*nnz(A) + n */
151 double maxwork ; /* maxwork for BTF, <= 0 if no limit */
152
153 Int btf ; /* use BTF pre-ordering, or not */
154 Int ordering ; /* 0: AMD, 1: COLAMD, 2: user P and Q,
155 * 3: user function */
156 Int scale ; /* row scaling: -1: none (and no error check),
157 * 0: none, 1: sum, 2: max */
158
159 /* memory management routines */
160 void *(*malloc_memory) (size_t) ; /* pointer to malloc */
161 void *(*realloc_memory) (void *, size_t) ; /* pointer to realloc */
162 void (*free_memory) (void *) ; /* pointer to free */
163 void *(*calloc_memory) (size_t, size_t) ; /* pointer to calloc */
164
165 /* pointer to user ordering function */
166 int (*user_order) (Int, Int *, Int *, Int *, struct klu_common<Entry, Int> *) ;
167
168 /* pointer to user data, passed unchanged as the last parameter to the
169 * user ordering function (optional, the user function need not use this
170 * information). */
171 void *user_data ;
172
173 Int halt_if_singular ; /* how to handle a singular matrix:
174 * FALSE: keep going. Return a Numeric object with a zero U(k,k). A
175 * divide-by-zero may occur when computing L(:,k). The Numeric object
176 * can be passed to klu_solve (a divide-by-zero will occur). It can
177 * also be safely passed to klu_refactor.
178 * TRUE: stop quickly. klu_factor will free the partially-constructed
179 * Numeric object. klu_refactor will not free it, but will leave the
180 * numerical values only partially defined. This is the default. */
181
182 /* ---------------------------------------------------------------------- */
183 /* statistics */
184 /* ---------------------------------------------------------------------- */
185
186 Int status ; /* KLU_OK if OK, < 0 if error */
187 Int nrealloc ; /* # of reallocations of L and U */
188
189 Int structural_rank ; /* 0 to n-1 if the matrix is structurally rank
190 * deficient (as determined by maxtrans). -1 if not computed. n if the
191 * matrix has full structural rank. This is computed by klu_analyze
192 * if a BTF preordering is requested. */
193
194 Int numerical_rank ; /* First k for which a zero U(k,k) was found,
195 * if the matrix was singular (in the range 0 to n-1). n if the matrix
196 * has full rank. This is not a true rank-estimation. It just reports
197 * where the first zero pivot was found. -1 if not computed.
198 * Computed by klu_factor and klu_refactor. */
199
200 Int singular_col ; /* n if the matrix is not singular. If in the
201 * range 0 to n-1, this is the column index of the original matrix A that
202 * corresponds to the column of U that contains a zero diagonal entry.
203 * -1 if not computed. Computed by klu_factor and klu_refactor. */
204
205 Int noffdiag ; /* # of off-diagonal pivots, -1 if not computed */
206
207 double flops ; /* actual factorization flop count, from klu_flops */
208 double rcond ; /* crude reciprocal condition est., from klu_rcond */
209 double condest ; /* accurate condition est., from klu_condest */
210 double rgrowth ; /* reciprocal pivot rgrowth, from klu_rgrowth */
211 double work ; /* actual work done in BTF, in klu_analyze */
212
213 size_t memusage ; /* current memory usage, in bytes */
214 size_t mempeak ; /* peak memory usage, in bytes */
215
216} ;
217
218/* -------------------------------------------------------------------------- */
219/* klu_defaults: sets default control parameters */
220/* -------------------------------------------------------------------------- */
221
222template <typename Entry, typename Int>
223Int klu_defaults
224(
225 klu_common<Entry, Int> *Common
226) ;
227
228/* -------------------------------------------------------------------------- */
229/* klu_analyze: orders and analyzes a matrix */
230/* -------------------------------------------------------------------------- */
231
232/* Order the matrix with BTF (or not), then order each block with AMD, COLAMD,
233 * a natural ordering, or with a user-provided ordering function */
234
235template <typename Entry, typename Int>
236klu_symbolic<Entry, Int> *klu_analyze
237(
238 /* inputs, not modified */
239 Int n, /* A is n-by-n */
240 Int Ap [ ], /* size n+1, column pointers */
241 Int Ai [ ], /* size nz, row indices */
242 klu_common<Entry, Int> *Common
243) ;
244
245/* -------------------------------------------------------------------------- */
246/* klu_analyze_given: analyzes a matrix using given P and Q */
247/* -------------------------------------------------------------------------- */
248
249/* Order the matrix with BTF (or not), then use natural or given ordering
250 * P and Q on the blocks. P and Q are interpretted as identity
251 * if NULL. */
252
253template <typename Entry, typename Int>
254klu_symbolic<Entry, Int> *klu_analyze_given
255(
256 /* inputs, not modified */
257 Int n, /* A is n-by-n */
258 Int Ap [ ], /* size n+1, column pointers */
259 Int Ai [ ], /* size nz, row indices */
260 Int P [ ], /* size n, user's row permutation (may be NULL) */
261 Int Q [ ], /* size n, user's column permutation (may be NULL) */
262 klu_common<Entry, Int> *Common
263) ;
264
265/* -------------------------------------------------------------------------- */
266/* klu_factor: factors a matrix using the klu_analyze results */
267/* -------------------------------------------------------------------------- */
268
269template <typename Entry, typename Int>
270klu_numeric<Entry, Int> *klu_factor /* returns KLU_OK if OK, < 0 if error */
271(
272 /* inputs, not modified */
273 Int Ap [ ], /* size n+1, column pointers */
274 Int Ai [ ], /* size nz, row indices */
275 Entry Ax [ ], /* size nz, numerical values */
276 klu_symbolic<Entry, Int> *Symbolic,
277 klu_common<Entry, Int> *Common
278) ;
279
280/* -------------------------------------------------------------------------- */
281/* klu_solve: solves Ax=b using the Symbolic and Numeric objects */
282/* -------------------------------------------------------------------------- */
283
284template <typename Entry, typename Int>
285Int klu_solve
286(
287 /* inputs, not modified */
288 klu_symbolic<Entry, Int> *Symbolic,
289 klu_numeric<Entry, Int> *Numeric,
290 Int ldim, /* leading dimension of B */
291 Int nrhs, /* number of right-hand-sides */
292
293 /* right-hand-side on input, overwritten with solution to Ax=b on output */
294 Entry B [ ], /* size ldim*nrhs */
295 klu_common<Entry, Int> *Common
296) ;
297
298
299/* -------------------------------------------------------------------------- */
300/* klu_tsolve: solves A'x=b using the Symbolic and Numeric objects */
301/* -------------------------------------------------------------------------- */
302
303template <typename Entry, typename Int>
304Int klu_tsolve
305(
306 /* inputs, not modified */
307 klu_symbolic<Entry, Int> *Symbolic,
308 klu_numeric<Entry, Int> *Numeric,
309 Int ldim, /* leading dimension of B */
310 Int nrhs, /* number of right-hand-sides */
311
312 /* right-hand-side on input, overwritten with solution to Ax=b on output */
313 Entry B [ ], /* size ldim*nrhs */
314 klu_common<Entry, Int> *Common
315) ;
316
317/* -------------------------------------------------------------------------- */
318/* klu_refactor: refactorizes matrix with same ordering as klu_factor */
319/* -------------------------------------------------------------------------- */
320
321template <typename Entry, typename Int>
322Int klu_refactor /* return TRUE if successful, FALSE otherwise */
323(
324 /* inputs, not modified */
325 Int Ap [ ], /* size n+1, column pointers */
326 Int Ai [ ], /* size nz, row indices */
327 Entry Ax [ ], /* size nz, numerical values */
328 klu_symbolic<Entry, Int> *Symbolic,
329 /* input, and numerical values modified on output */
330 klu_numeric<Entry, Int> *Numeric,
331 klu_common<Entry, Int> *Common
332) ;
333
334/* -------------------------------------------------------------------------- */
335/* klu_free_symbolic: destroys the Symbolic object */
336/* -------------------------------------------------------------------------- */
337
338template <typename Entry, typename Int>
339Int klu_free_symbolic
340(
341 klu_symbolic<Entry, Int> **Symbolic,
342 klu_common<Entry, Int> *Common
343) ;
344
345/* -------------------------------------------------------------------------- */
346/* klu_free_numeric: destroys the Numeric object */
347/* -------------------------------------------------------------------------- */
348
349/* Note that klu_free_numeric and klu_z_free_numeric are identical; each can
350 * free both kinds of Numeric objects (real and complex) */
351
352template <typename Entry, typename Int>
353Int klu_free_numeric
354(
355 klu_numeric<Entry, Int> **Numeric,
356 klu_common<Entry, Int> *Common
357) ;
358
359/* -------------------------------------------------------------------------- */
360/* klu_sort: sorts the columns of the LU factorization */
361/* -------------------------------------------------------------------------- */
362
363/* this is not needed except for the MATLAB interface */
364
365template <typename Entry, typename Int>
366Int klu_sort
367(
368 /* inputs, not modified */
369 klu_symbolic<Entry, Int> *Symbolic,
370 /* input/output */
371 klu_numeric<Entry, Int> *Numeric,
372 klu_common<Entry, Int> *Common
373) ;
374
375/* -------------------------------------------------------------------------- */
376/* klu_flops: determines # of flops performed in numeric factorzation */
377/* -------------------------------------------------------------------------- */
378
379template <typename Entry, typename Int>
380Int klu_flops
381(
382 /* inputs, not modified */
383 klu_symbolic<Entry, Int> *Symbolic,
384 klu_numeric<Entry, Int> *Numeric,
385 /* input/output */
386 klu_common<Entry, Int> *Common
387) ;
388
389/* -------------------------------------------------------------------------- */
390/* klu_rgrowth : compute the reciprocal pivot growth */
391/* -------------------------------------------------------------------------- */
392
393/* Pivot growth is computed after the input matrix is permuted, scaled, and
394 * off-diagonal entries pruned. This is because the LU factorization of each
395 * block takes as input the scaled diagonal blocks of the BTF form. The
396 * reciprocal pivot growth in column j of an LU factorization of a matrix C
397 * is the largest entry in C divided by the largest entry in U; then the overall
398 * reciprocal pivot growth is the smallest such value for all columns j. Note
399 * that the off-diagonal entries are not scaled, since they do not take part in
400 * the LU factorization of the diagonal blocks.
401 *
402 * In MATLAB notation:
403 *
404 * rgrowth = min (max (abs ((R \ A(p,q)) - F)) ./ max (abs (U))) */
405
406template <typename Entry, typename Int>
407Int klu_rgrowth
408(
409 Int Ap [ ],
410 Int Ai [ ],
411 Entry Ax [ ],
412 klu_symbolic<Entry, Int> *Symbolic,
413 klu_numeric<Entry, Int> *Numeric,
414 klu_common<Entry, Int> *Common /* Common->rgrowth = reciprocal pivot growth */
415) ;
416
417/* -------------------------------------------------------------------------- */
418/* klu_condest */
419/* -------------------------------------------------------------------------- */
420
421/* Computes a reasonably accurate estimate of the 1-norm condition number, using
422 * Hager's method, as modified by Higham and Tisseur (same method as used in
423 * MATLAB's condest */
424
425template <typename Entry, typename Int>
426Int klu_condest
427(
428 Int Ap [ ], /* size n+1, column pointers, not modified */
429 Entry Ax [ ], /* size nz = Ap[n], numerical values, not modified*/
430 klu_symbolic<Entry, Int> *Symbolic, /* symbolic analysis, not modified */
431 klu_numeric<Entry, Int> *Numeric, /* numeric factorization, not modified */
432 klu_common<Entry, Int> *Common /* result returned in Common->condest */
433) ;
434
435/* -------------------------------------------------------------------------- */
436/* klu_rcond: compute min(abs(diag(U))) / max(abs(diag(U))) */
437/* -------------------------------------------------------------------------- */
438
439template <typename Entry, typename Int>
440Int klu_rcond
441(
442 klu_symbolic<Entry, Int> *Symbolic, /* input, not modified */
443 klu_numeric<Entry, Int> *Numeric, /* input, not modified */
444 klu_common<Entry, Int> *Common /* result in Common->rcond */
445) ;
446
447/* -------------------------------------------------------------------------- */
448/* klu_scale */
449/* -------------------------------------------------------------------------- */
450
451template <typename Entry, typename Int>
452Int klu_scale /* return TRUE if successful, FALSE otherwise */
453(
454 /* inputs, not modified */
455 Int scale, /* <0: none, no error check; 0: none, 1: sum, 2: max */
456 Int n,
457 Int Ap [ ], /* size n+1, column pointers */
458 Int Ai [ ], /* size nz, row indices */
459 Entry Ax [ ],
460 /* outputs, not defined on input */
461 double Rs [ ],
462 /* workspace, not defined on input or output */
463 Int W [ ], /* size n, can be NULL */
464 klu_common<Entry, Int> *Common
465) ;
466
467/* -------------------------------------------------------------------------- */
468/* klu_extract */
469/* -------------------------------------------------------------------------- */
470
471template <typename Entry, typename Int>
472Int klu_extract /* returns TRUE if successful, FALSE otherwise */
473(
474 /* inputs: */
475 klu_numeric<Entry, Int> *Numeric,
476 klu_symbolic<Entry, Int> *Symbolic,
477
478 /* outputs, either allocated on input, or ignored otherwise */
479
480 /* L */
481 Int *Lp, /* size n+1 */
482 Int *Li, /* size Numeric->lnz */
483 Entry *Lx, /* size Numeric->lnz */
484
485 /* U */
486 Int *Up, /* size n+1 */
487 Int *Ui, /* size Numeric->unz */
488 Entry *Ux, /* size Numeric->unz */
489
490 /* F */
491 Int *Fp, /* size n+1 */
492 Int *Fi, /* size Numeric->nzoff */
493 Entry *Fx, /* size Numeric->nzoff */
494
495 /* P, row permutation */
496 Int *P, /* size n */
497
498 /* Q, column permutation */
499 Int *Q, /* size n */
500
501 /* Rs, scale factors */
502 Entry *Rs, /* size n */
503
504 /* R, block boundaries */
505 Int *R, /* size Symbolic->nblocks+1 (nblocks is at most n) */
506
507 klu_common<Entry, Int> *Common
508) ;
509
510
511/* -------------------------------------------------------------------------- */
512/* KLU memory management routines */
513/* -------------------------------------------------------------------------- */
514
515template <typename Entry, typename Int>
516void *klu_malloc /* returns pointer to the newly malloc'd block */
517(
518 /* ---- input ---- */
519 size_t n, /* number of items */
520 size_t size, /* size of each item */
521 /* --------------- */
522 klu_common<Entry, Int> *Common
523) ;
524
525template <typename Entry, typename Int>
526void *klu_free /* always returns NULL */
527(
528 /* ---- in/out --- */
529 void *p, /* block of memory to free */
530 size_t n, /* number of items */
531 size_t size, /* size of each item */
532 /* --------------- */
533 klu_common<Entry, Int> *Common
534) ;
535
536template <typename Entry, typename Int>
537void *klu_realloc /* returns pointer to reallocated block */
538(
539 /* ---- input ---- */
540 size_t nnew, /* requested # of items in reallocated block */
541 size_t nold, /* current size of block, in # of items */
542 size_t size, /* size of each item */
543 /* ---- in/out --- */
544 void *p, /* block of memory to realloc */
545 /* --------------- */
546 klu_common<Entry, Int> *Common
547) ;
548
549/* ========================================================================== */
550/* === KLU version ========================================================== */
551/* ========================================================================== */
552
553/* All versions of KLU include these definitions.
554 * As an example, to test if the version you are using is 1.2 or later:
555 *
556 * if (KLU_VERSION >= KLU_VERSION_CODE (1,2)) ...
557 *
558 * This also works during compile-time:
559 *
560 * #if (KLU >= KLU_VERSION_CODE (1,2))
561 * printf ("This is version 1.2 or later\n") ;
562 * #else
563 * printf ("This is an early version\n") ;
564 * #endif
565 */
566
567#define KLU_DATE "Mar 24, 2009"
568#define KLU_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
569#define KLU_MAIN_VERSION 1
570#define KLU_SUB_VERSION 1
571#define KLU_SUBSUB_VERSION 0
572#define KLU_VERSION KLU_VERSION_CODE(KLU_MAIN_VERSION,KLU_SUB_VERSION)
573
574#endif
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.