/[ascend]/trunk/base/generic/linear/linsolqr.h
ViewVC logotype

Contents of /trunk/base/generic/linear/linsolqr.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1395 - (show annotations) (download) (as text)
Sat Apr 21 15:11:45 2007 UTC (15 years, 2 months ago) by jpye
File MIME type: text/x-chdr
File size: 39874 byte(s)
Added full jacobian preconditioner (ongoing).
Some tweaks to dsgsat & related models.
1 /* ASCEND modelling environment
2 Copyright (C) 1990 Karl Michael Westerberg
3 Copyright (C) 1993 Joseph Zaher
4 Copyright (C) 1994 Boyd Safrit, Joseph Zaher, Benjamin Andrew Allan
5 Copyright (C) 1995 Benjamin Andrew Allan, Kirk Andre Abbott
6 Copyright (C) 2006-2007 Carnegie Mellon University
7
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2, or (at your option)
11 any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place - Suite 330,
21 Boston, MA 02111-1307, USA.
22 *//** @file
23 * linsolqr: Ascend Linear Equation Solver.
24 *
25 * A linear system consists of a coefficient matrix (A)
26 * and possibly several right-hand-sides (rhs). The
27 * solution vector x sought can be that of either
28 *<pre>
29 * T
30 * A x = rhs or A x = rhs
31 *</pre>
32 * depending on a specification inherent with rhs which
33 * dictates whether or not A should be transposed. If a
34 * rhs specifies transpose, then the vector itself is
35 * expected to be indexed by the original column numbers
36 * of the coefficient matrix and the solution vector shall
37 * be indexed by original row numbers. Otherwise, rhs
38 * is expected to be indexed by original row numbers while
39 * the solution can be referenced using original column
40 * numbers. The coefficient matrix and each rhs will be
41 * preserved throughout solving, except that the
42 * coefficient matrix may be permuted during reordering
43 * or (depending on the method) solving.
44 *
45 * The method of linear solution is an option new to
46 * linsolqr. Also, some of the data structures have been
47 * expanded as required. The factoring methods are
48 * enumerated and summarized below.
49 *
50 * Linsol is a general sparse solver. It uses mtx, though any
51 * proper sparse matrix package can be substituted. By
52 * proper, of course, we mean one with at least the
53 * functionality of mtx. If you find something with all
54 * the features of mtx that is faster than mtx, let us
55 * know, please. If you do any comparisons of linsol with
56 * other sparse codes (on sequential processors) also
57 * please let us know how you find its performance.
58 * You probably won't find algorithms for banded or
59 * symmetric matrices here. Why?
60 * The efficiency of most banded algorithms comes from
61 * the way they access memory with three assumptions:
62 * that all possible coefficient locations are known before
63 * a solution is attempted, that they are accessible
64 * in some ORDERED fashion, and that the matrix will not
65 * be reordered (except perhaps in very limited ways)
66 * during solution. Linsol is predicated on just
67 * the opposite assumptions at the moment.
68 * These assumptions would relieve the solver of
69 * memory allocation and a great deal of range checking.
70 * While we would do our flops efficiently with linsol/mtx
71 * on a banded matrix, the cost of accessing linear
72 * coefficients in an UNORDERED fashion is such that a good
73 * band, variable band, or symmetric skyline code making the
74 * ORDERED assumption will always beat linsol in speed,
75 * though perhaps not in robustness.
76 *
77 * Requires: #include <string.h>
78 * #include "utilities/ascConfig.h"
79 * #include "mtx.h"
80 *
81 * Each method is described here, including the unimplemented ones,
82 * however vaguely. By convention all methods assign unpivotable variables
83 * the value of 0.0 in the linear solution vector. This makes sense because
84 * the linear system usually arises from a nonlinear scheme wherein
85 * stepping in the unpivoted variables is left for a future iteration.
86 *
87 * The currently recommended methods are ranki_*2, in particular ba2.
88 *
89 * ranki_kw:
90 * Reordering method: SPK1, Stadtherr & Wood (1984), and like.
91 * The region used is confined to the main diagonal.
92 * The square region is reordered to lower triangular with spikes.
93 * Each spike is higher than or equal to all spikes to the left of it
94 * if the matrix is partitioned.
95 * Solution method: RANKI, Stadtherr & Wood
96 * Ref: Computers and Chemical Engineering (8)1 pp. 9-33 (1984)
97 * The square region is U/L factored with spikes dragged in to replace
98 * unacceptable pivots. Fill is in spike columns and
99 * dragged rows where it usually does take place.
100 * Numerically dependent rows and columns do not have dependencies
101 * automatically calculated. This is a change from linsol, so be
102 * careful in making comparisons.
103 * Authors: Karl Westerberg, Joe Zaher
104 * with some clean ups and restructuring by Ben Allan
105 *
106 * ranki_jz:
107 * As ranki_kw, except where ranki_kw looks only across the current row
108 * for pivots if the diagonal element is too small, ranki_jz also looks
109 * down the current column. This introduces additional spikes to the
110 * current block, but sometimes the overall effect is reduced fill. The
111 * main motivation is to broaden pivot choices on those problems
112 * where restricting the choice to the current row results in immediate
113 * explosions of fill due to a full subdiagonal and a moderately small
114 * pivot.
115 *
116 * ranki_kw2:
117 * ranki_jz2:
118 * As the ranki_XX factorizations, except uses drop tolerance and
119 * uses much faster bilevel matrices (master + slave matrix of mtx)
120 * to segregate the contents of U from L. Also, dependency information
121 * is not calculated automatically.
122 * Author: Ben Allan
123 *
124 * ranki_ba2:
125 * As the ranki_kw2 factorization, except pivot and row being eliminated
126 * are stored as linked lists/full org vectors (dual access) during
127 * elimination which removes the last significant quadratic effect from
128 * the ranki factorization.
129 *
130 * plain_qr:
131 * Reordering method: Transpose SPK1 (tspk1), LORA
132 * The region used may be rectangular, but will be modified
133 * if needed so that the upper left corner is on the diagonal.
134 * The region in sys->coef is not permuted in the solution process.
135 * Solution method:
136 * The rectangular region is QR factored with standard column pivoting
137 * modified by the pivot_tolerance. That is, at each step the column
138 * to factor in next is the leftmost col with norm alpha that satisfies
139 * alpha >= ptol * maxalpha.
140 * ptol = 1 --> Stewart's QRDC style pivoting.
141 * ptol = 0 --> no pivoting (Sometimes a dumb idea).
142 * plain_qr is to QRDC as ranki_kw is to LU factorization with strict
143 * partial pivoting.
144 * This implementation is best used on matrices where nrows about = ncols
145 * and data are not highly correlated.
146 * There are better QR methods (not available here) for vector processors,
147 * distributed processors, data fitting and so forth.
148 * Some care has been taken that the implementation here will scale up
149 * reasonably to larger matrices.
150 * Detailed references are given in the plain_qr section of the .c file.
151 * Author: Ben Allan
152 *
153 * cond_qr: (broken)
154 * Reordering method: Transpose SPK1 (tspk1)
155 * The region used is confined to the main diagonal.
156 * The square region is permuted in the solution process.
157 * Solution method:
158 * The square region is QR factored with pivoting by a local minimum
159 * fill criteria balanced against choosing pivots based on the
160 * incremental inverse R condition number.
161 * Ref:Incremental Condition Calculation and Column Selection,
162 * UMIACS TR 90-87, G.W.Stewart, July 1990)
163 * (Allan, Safrit, Westerberg, 1995)
164 * Authors: Ben Allan
165 *
166 * opt_qr:
167 * As cond_qr, except region may be underspecified (MxN, M < N)
168 * The region's upper left corner should be (or will be forced to be)
169 * on the main diagonal. The columns pivoted are expected to be a
170 * good basis. Solution is not a least squares solution of the
171 * transpose NxM problem; for this we need rowwise Householder
172 * transforms and spk1 reordering.
173 * Authors: Ben Allan
174 *
175 * ls_qr:
176 * As cond_qr, except region may be overspecified (MxN, M > N)
177 * Solves the system in a linear least squares sense.
178 * Authors: Ben Allan
179 *
180 * band_lu:
181 * symmetric_lu:
182 * (cholesky?)
183 * The best codes we can find already implemented to complement and/or
184 * verify our own. More detailed semantics unknown at this time.
185 *
186 * The parameters pivot_tolerance, condition_tolerance and drop_tolerance
187 * have semantics varying with the method.
188 * See the header of linsolqr_set_pivot_tolerance below for details.
189 * </pre>
190 *//*
191 by Karl Westerberg
192 Created: 2/6/90
193 Last in CVS: $Revision: 1.13 $ $Date: 1997/07/25 17:23:49 $ $Author: ballan $
194 QR options by Ben Allan
195 C version based on sqr.pas v1.5, Boyd Safrit (1994)
196
197 Dates: 06/90 - original version (KW)
198 04/91 - removed output assignment and partitioning
199 (belong in structural analysis) (JZ)
200 08/92 - added transpose ability (JZ)
201 01/94 - broke out linsol_invert() and linsol_solve() (JZ)
202 10/94 - reorganized for multiple methods and QR (BA)
203 moved drag operator to mtx (BA,JZ)
204 11/94 - added linsolqr_get_factors/inverse for debugging.BA
205 09/95 - added factor_class to better manage families. BA
206 09/95 - added ranki_*2 methods (BA)
207 09/95 - added gauss (broken) (KA)
208 11/96 - added ranki_ba2 method (BA)
209
210 */
211
212 #ifndef ASC_LINSOLQR_H
213 #define ASC_LINSOLQR_H
214
215 #include "mtx.h"
216
217 /** @addtogroup linear Linear
218 @{
219 */
220
221 #define LINSOLMTX_DEBUG FALSE
222 /**< Debug mode. */
223 #define LINQR_DROP_TOLERANCE 1.0e-16
224 /**< This is the default for drop tolerances in methods which use a drop tol */
225
226 struct rhs_list;
227
228 typedef struct linsolqr_header *linsolqr_system_t;
229 /**< linsolqr_system_t is the linear system handle. */
230
231 enum factor_class {
232 unknown_c, /**< error handling class */
233 ranki = 100, /**< all ranki (and gauss until broken out) methods */
234 s_qr = 200 /**< all sparse qr methods */
235 };
236
237 enum reorder_method {
238 unknown_r, /**< error handling method */
239 natural = 1000, /**< do nothing reorder */
240 spk1 = 2000, /**< Stadtherr's SPK1 reordering */
241 tspk1 = 3000 /**< transpose of Stadtherr's SPK1 reordering good for gauss */
242 /* future work:
243 invspk1, spk1 then diagonally inverted
244 invtspk1, tspk1 then diagonally inverted
245 widespk1, spk1 of an MxN region, N > M
246 */
247 };
248
249 enum factor_method {
250 unknown_f = 0, /**< error handling method */
251 ranki_kw = 1, /**< original linsol method */
252 ranki_jz = 2, /**< original linsol method with pseudo-complete pivoting */
253 ranki_kw2 = 3, /**< better data structure version of ranki_kw, w/drop tol */
254 ranki_jz2 = 4, /**< better data structure version of ranki_jz, w/drop tol */
255 ranki_ka = 12, /**< kirks hacked verion of the basic method. */
256 plain_qr = 5, /**< rectangular col pivoted qr variants */
257 cond_qr = 6, /**< stewart-based sparse QR method new Coke */
258 ranki_ba2 = 7, /**< proper linked list implementation of ranki (dragfree) */
259 opt_qr = 8, /**< coming soon */
260 ls_qr = 9, /**< anticipated */
261 gauss_ba2 = 10, /**< anticipated */
262 symmetric_lu = 11 /**< anticipated */
263 };
264
265 ASC_DLLSPEC int g_linsolqr_timing;
266 /**<
267 * If nonzero and certain internal defines are set, factorization
268 * generates a fill and timing message.
269 */
270
271 /* KAA_DEBUG */
272 /* unimplemented function
273 extern boolean linsolqr_col_is_a_spike (mtx_matrix_t mtx, int32 col);
274 */
275 /* Returns true if the column in question is a spike. */
276
277
278 /* Functions for managing interfaces */
279
280 extern char *linsolqr_rmethods(void);
281 /**<
282 * Returns a comma-separated list of the names of reordering methods
283 * implemented. If you implement a new method, update this
284 * function. Do not free the pointer returned.<br><br>
285 *
286 * Not all reorderings are appropriate to all factorizations.
287 */
288
289 ASC_DLLSPEC char *linsolqr_fmethods(void);
290 /**<
291 * Returns a comma-separated list of the names of factorization methods
292 * implemented. If you implement a new method, update this
293 * function. Do not free the pointer returned. The string returned
294 * will contain no whitespace of any sort, tabs, blanks, or \n.
295 */
296
297 extern enum reorder_method linsolqr_rmethod_to_enum(char *s);
298 /**<
299 * Returns the enum of a reorder method with the name s.
300 * If you implement a new method, update this function.
301 */
302
303 extern enum factor_method linsolqr_fmethod_to_enum(char *s);
304 /**<
305 * Returns the enum of a factor method with the name s.
306 * If you implement a new method, update this function.
307 */
308
309 ASC_DLLSPEC enum factor_class linsolqr_fmethod_to_fclass(enum factor_method fm);
310 /**<
311 * Returns the enum of the factor class containing the method given.
312 * If you implement a new method or class, update this function.
313 */
314
315 ASC_DLLSPEC char *linsolqr_enum_to_rmethod(enum reorder_method m);
316 /**<
317 * Returns the name of a reorder method with the enum m.
318 * If you implement a new method, update this function.
319 * Do not free the name.
320 */
321
322 extern char *linsolqr_enum_to_fmethod(enum factor_method m);
323 /**<
324 * Returns the name of a factor method with the enum m.
325 * If you implement a new method, update this function.
326 * Do not free the name.
327 */
328
329 ASC_DLLSPEC const char *linsolqr_rmethod_description(const enum reorder_method meth);
330 /**<
331 * Returns a string describing the method inquired on. Do not mess
332 * with the string: linsolqr owns it.
333 * If you implement a new method, update this function.
334 */
335
336 ASC_DLLSPEC const char *linsolqr_fmethod_description(const enum factor_method meth);
337 /**<
338 * Returns a string describing the method inquired on. Do not mess
339 * with the string: linsolqr owns it.
340 * If you implement a new method, update this function.
341 */
342
343 /* Functions for specifying problems and controlling them */
344
345 ASC_DLLSPEC linsolqr_system_t linsolqr_create(void);
346 /**<
347 * Creates a linear system and returns a pointer to it. Initially the
348 * system has no coefficient matrix and no rhs.
349 */
350
351 ASC_DLLSPEC linsolqr_system_t linsolqr_create_default(void);
352 /**<
353 Creates a linsolqr system as per linsolqr_default but sets default
354 'recommended' preferences for reordering and factoring methods, so the
355 ignorant user (me) doesn't have to worry about it every time they use this.
356 */
357
358 ASC_DLLSPEC void linsolqr_destroy(linsolqr_system_t sys);
359 /**<
360 * Destroys the linear system. The coefficient matrix and each rhs are
361 * not destroyed by this call.
362 */
363
364 ASC_DLLSPEC void linsolqr_set_matrix(linsolqr_system_t sys, mtx_matrix_t mtx);
365 /**<
366 * Sets the coefficient matrix to mtx.
367 */
368
369 ASC_DLLSPEC void linsolqr_set_region(linsolqr_system_t sys, mtx_region_t region);
370 /**<
371 * Sets the reg to region.
372 */
373
374 ASC_DLLSPEC mtx_matrix_t linsolqr_get_matrix(linsolqr_system_t sys);
375 /**<
376 * Returns the coefficient matrix.
377 */
378
379 ASC_DLLSPEC void linsolqr_add_rhs(linsolqr_system_t sys,
380 real64 *rhs,
381 boolean transpose);
382 /**<
383 * Adds the given rhs to the collection of rhs's already part of the
384 * system.
385 *
386 * You continue to be responsible for deallocating rhs. We do not take
387 * over management of that memory. Do not free the rhs you have given
388 * us without calling linsolqr_remove_rhs first.<br><br>
389 *
390 * Rhs should point to an array of reals indexed by original
391 * column number if the linear system is to be solved using the transpose
392 * of the matrix or by original row number if the matrix is not to be
393 * transposed. This is determined using the boolean transpose. The
394 * rhs should be refered to in the future by its pointer.
395 * <pre>
396 * For the mathematically impaired:
397 * M' denotes Transpose(M).
398 *
399 * transpose==FALSE --> rhs will be solved for x in A x = b.
400 * --> rhs b is an org row indexed column vector.
401 * --> x is an org col indexed column vector.
402 * --> yields dx = Newton direction in J.dx = -f
403 * if the b given is -f.
404 *
405 * transpose==TRUE --> rhs will be solved for y in A' y = d.
406 * --> rhs d is an org column indexed column vector.
407 * --> y is an org row indexed column vector.
408 * --> yields dx = Newton direction in J'.dx = -f.
409 * if the d given is -f and f,dx are properly indexed.
410 * This may be useful if A was created using
411 * mtx_transpose in order to improve factorization.
412 *
413 * Useful matrix identity: (M')^-1 == (M^-1)' for invertible M.
414 * </pre>
415 */
416
417 ASC_DLLSPEC void linsolqr_remove_rhs(linsolqr_system_t sys, real64 *rhs);
418 /**<
419 * Removes the given rhs from the system. The rhs is not destroyed, just
420 * removed from the system.
421 */
422
423 extern int32 linsolqr_number_of_rhs(linsolqr_system_t sys);
424 /**<
425 * Returns the number of rhs's currently part of the system.
426 */
427
428 ASC_DLLSPEC real64 *linsolqr_get_rhs(linsolqr_system_t sys, int n);
429 /**<
430 * Returns the n-th rhs, where rhs's are indexed in the order they were
431 * added using linsolqr_add_rhs() from 0 to (# rhs's)-1. NULL is returned
432 * if the index is out of range.
433 */
434
435 ASC_DLLSPEC void linsolqr_matrix_was_changed(linsolqr_system_t sys);
436 /**<
437 * Informs the solver that a numerical value of a non-zero was changed.
438 * This must be called whenever any numerical changes to the matrix are
439 * made.
440 */
441
442 ASC_DLLSPEC void linsolqr_rhs_was_changed(linsolqr_system_t sys,
443 real64 *rhs);
444 /**<
445 * Informs the solver that the given rhs has been modified. This must be
446 * called whenever the rhs is modified.
447 */
448
449 extern void linsolqr_set_pivot_zero(linsolqr_system_t sys,
450 real64 pivot_zero);
451 /**<
452 * Sets the pivot zero for the system. Pivots less than or equal to
453 * this value are regarded as zero. linsolqr_set_pivot_zero() will
454 * automatically call linsolqr_matrix_was_changed().
455 */
456 extern real64 linsolqr_pivot_zero(linsolqr_system_t sys);
457 /**<
458 * Gets the pivot zero for the system. Pivots less than or equal to
459 * this value are regarded as zero.
460 */
461
462 extern void linsolqr_set_pivot_tolerance(linsolqr_system_t sys, real64 ptol);
463 /**< See discussion under linsolqr_drop_tolerance(). */
464 extern real64 linsolqr_pivot_tolerance(linsolqr_system_t sys);
465 /**< See discussion under linsolqr_drop_tolerance(). */
466 extern void linsolqr_set_condition_tolerance(linsolqr_system_t sys, real64 ctol);
467 /**< See discussion under linsolqr_drop_tolerance(). */
468 extern real64 linsolqr_condition_tolerance(linsolqr_system_t sys);
469 /**< See discussion under linsolqr_drop_tolerance(). */
470 extern void linsolqr_set_drop_tolerance(linsolqr_system_t sys, real64 dtol);
471 /**< See discussion under linsolqr_drop_tolerance(). */
472 extern real64 linsolqr_drop_tolerance(linsolqr_system_t sys);
473 /**<
474 * <pre>
475 * linsolqr_set_pivot_tolerance(sys,ptol)
476 * ptol = linsolqr_pivot_tolerance(sys)
477 *
478 * linsolqr_set_condition_tolerance(sys,ctol)
479 * ctol = linsolqr_condition_tolerance(sys)
480 *
481 * linsolqr_set_drop_tolerance(sys,dtol)
482 * dtol = linsolqr_drop_tolerance(sys)
483 *
484 * linsolqr_system_t sys;
485 * real64 ptol, ctol,dtol;
486 *
487 * Sets/gets the pivot/condition/drop tolerance for the system. Semantics of
488 * tolerances vary with method. Not to be confused with pivot zero, epsilon.
489 *
490 * pivot_tolerance: Pivots less than this fraction of the maximum pivot
491 * value in the same row or column (depending on method) are disregarded.
492 *
493 * ranki_kw: pivot_tolerance applies to row. condition_tolerance ignored.
494 * ranki_jz: pivot_tolerance applies to row and col. condition ignored.
495 * ranki_ba2:
496 * ranki_kw2:
497 * ranki_jz2: as ranki_kw/ranki_jz except that matrix entries below
498 * dtol in magnitude are dropped.
499 * plain_qr: pivot_tolerance applies to cols. condition_tolerance used
500 * with a condition heuristic rather than actual condition.
501 *
502 * cond_qr and variants:
503 * Unpivoted columns are ordered by fill potential and then the first
504 * of these to meet the criterion CI*condition_tolerance <= min_CI
505 * is chosen as the next column to pivot with.
506 * Pivot_tolerance is applied to choose the pivot within the
507 * selected column.
508 *
509 * linsolqr_set_pivot/condition/drop_tolerance() will automatically call
510 * linsolqr_matrix_was_changed().
511 * </pre>
512 * @todo Separate documentation for these linsolqr functions?
513 */
514
515 /* Functions for analyzing and querying linear systems. */
516
517 extern enum factor_class linsolqr_fclass(linsolqr_system_t sys);
518 /**<
519 * Returns the most recently set factorization class of the system.
520 * The system should be previously prepped.
521 */
522
523 ASC_DLLSPEC enum factor_method linsolqr_fmethod(linsolqr_system_t sys);
524 /**<
525 * Returns the most recently used factorization method of the system.
526 * The system should be previously prepped.
527 */
528
529 ASC_DLLSPEC enum reorder_method linsolqr_rmethod(linsolqr_system_t sys);
530 /**<
531 * Returns the most recently set reorder method of the system.
532 * The system should be previously prepped.
533 */
534
535 ASC_DLLSPEC int32 linsolqr_rank(linsolqr_system_t sys);
536 /**<
537 * Returns the rank of the system. The system must be previously
538 * factored.
539 */
540
541 extern real64 linsolqr_smallest_pivot(linsolqr_system_t sys);
542 /**<
543 * Returns the smallest pivot accepted in solving the system. The
544 * system must be previously factored. If no pivoting was performed,
545 * MAXDOUBLE is returned.
546 */
547
548 ASC_DLLSPEC int linsolqr_prep(linsolqr_system_t sys, enum factor_class fclass);
549 /**<
550 * This function is analogous to slv_select_solver. It
551 * takes care of allocations. The prep call should be done (once) at
552 * the beginning of any series of linear solutions being performed on
553 * on the same linsolqr system with the same factorization class.<br><br>
554 *
555 * You must prep before doing a reordering, factorization or solution
556 * as prep sets up the appropriate internal settings. You should set
557 * the matrix before you call prep.
558 * Prep (or subsequent preps) do not affect the right hand sides of
559 * the system.<br><br>
560 *
561 * If you wish to change from 1 factoring method to another and they
562 * are not part of the class of compatible methods, you should call
563 * prep again with the new class.
564 * After prep is called, reorder should be called before factorization.
565 * As most of the sparse routines depend for performance on the
566 * prior reordering.<br><br>
567 *
568 * Return 0 if ok, or 1 otherwise.
569 */
570
571 ASC_DLLSPEC int linsolqr_reorder(linsolqr_system_t sys,
572 mtx_region_t *region,
573 enum reorder_method rmeth);
574 /**<
575 * The specified region of the coefficient matrix is reordered. This
576 * must be called before factoring the matrix. The specified region
577 * is assumed to contain only nonempty rows and columns and have a full
578 * diagonal.
579 * If the coefficient matrix in use is from a nonlinear system, the
580 * pattern in the coefficient matrix should be the structural one (as
581 * opposed to the numerically derived incidence which may be less.)<br><br>
582 *
583 * If you use the numerically derived incidence, you will need to reorder
584 * before every factorization. This is generally not cost effective.
585 * If region given is mtx_ENTIRE_MATRIX, a search will be done to find
586 * an appropriate bounding region in the coefficient mtx. This is
587 * not a particularly cheap search.
588 *
589 * This function is analogous to slv_presolve. It
590 * takes care of any structural analysis necessary
591 * to the linear solution. The reorder call should be done (once) at
592 * the beginning of any series of linear solutions being performed on
593 * on structurally identical matrices.
594 *
595 * The reordering done must be appropriate to the factorization class.
596 * You must reorder before doing a factorization or solution as reorder
597 * sets up the appropriate internal settings. Even were the internals
598 * method independent, the reordering is critical to the performance of
599 * these methods.
600 *
601 * @return 0 if ok, 1 if not.
602 * <pre>
603 * Brief notes on the reorderings available.
604 * SPK1: The region you give becomes the region for the problem,
605 * TSPK1: but for factorization purposes rows and columns that
606 * do not intersect the main diagonal within the region
607 * are considered structurally (therefore numerically) dependent.
608 * Natural: Blesses the system and does nothing.
609 * Again, the rows/cols not in the diagonal are dependent.
610 * </pre>
611 */
612
613 ASC_DLLSPEC int linsolqr_factor(linsolqr_system_t sys,
614 enum factor_method fmethod);
615 /**<
616 * Decompose the reordered region of a copy of the coefficient matrix
617 * into upper and lower triangular factors (if necessary) which can be
618 * inverted easily when applied to any rhs. Matrix must be factored in
619 * order to perform any of the operations below. The numerical rank and
620 * the smallest pivot encountered during pivoting are computed in the
621 * process. Factorization method used is that given if it is compatible
622 * with the class set when the prep was done, otherwise the call fails.
623 *
624 * If you are handed a linsolqr system and want to refactor it using the
625 * usual method, but don't know what that method is, call like:
626 * status = linsolqr_factor(sys,linsolqr_fmethod(sys));
627 *
628 * linsolqr_reorder must be called before factoring the matrix -- JP -- check?
629 *
630 * Return 0 if ok, 1 if not.
631 */
632
633 extern int linsolqr_get_pivot_sets(linsolqr_system_t sys,
634 unsigned *org_rowpivots,
635 unsigned *org_colpivots);
636 /**<
637 * Returns the set of original row numbers / original column numbers which
638 * have been pivoted. org_rowpivots and org_colpivots are assumed to be
639 * sets created by (or at least for) the set module with sufficient size
640 * before calling this function. They must also be previously NULLed.
641 * The system must be previously factored.
642 * The sets input should be the result of set_create(neqn),set_create(nvar).
643 * There is no association of rows with columns here.
644 *
645 * @return Status is 0 if not factored, 1 if factored. If 0, sets will not be
646 * changed.
647 * This bizarre piece of functionality should be done away with as soon
648 * as the equivalents below have been implemented.
649 */
650
651 ASC_DLLSPEC mtx_sparse_t *linsolqr_unpivoted_rows(linsolqr_system_t sys);
652 /**< See discussion under linsolqr_unpivoted_cols(). */
653 ASC_DLLSPEC mtx_sparse_t *linsolqr_unpivoted_cols(linsolqr_system_t sys);
654 /**<
655 * Returns the set of original row numbers / original column numbers which
656 * have NOT been pivoted and the rejected pivots in an mtx_sparse_t.
657 * Return is NULL if sys not factored or if no unpivoted rows/cols exist.
658 * E.g. singrows->idata[i] are the original rows that were not pivoted,
659 * singrows->data[i] is the pivot that was rejected,
660 * for i = 0 to singrows->len-1.
661 * If len is 0, the data and idata pointers may be NULL.
662 *
663 * The CALLER is responsible for deallocating the mtx_sparse_t returned;
664 * linsolqr wants nothing further to do with it.
665 * @bug Only deals with ranki based square systems at the moment.
666 * Then again, that's all we have at the moment (10/95).
667 * Returns NULL from non-ranki systems.
668 * @todo Separate documentation for these linsolqr functions?
669 */
670
671 extern mtx_sparse_t *linsolqr_pivoted_rows(linsolqr_system_t sys);
672 /**< See discussion under linsolqr_pivoted_cols(). */
673 extern mtx_sparse_t *linsolqr_pivoted_cols(linsolqr_system_t sys);
674 /**<
675 * Returns the set of original row numbers / original column numbers which
676 * have been pivoted and the pivots in an mtx_sparse_t.
677 * Return is NULL if sys not factored or if no pivoted rows/cols exist.
678 * E.g. pivrows->idata[i] are the original rows that were pivoted,
679 * pivrows->data[i] is the pivot,
680 * for i = 0 to pivrows->len-1.
681 * If len is 0, the data and idata pointers may be NULL.
682 *
683 * The CALLER is responsible for deallocating the mtx_sparse_t returned;
684 * linsolqr wants nothing further to do with it.
685 * @bug Only deals with ranki based square systems at the moment.
686 * Then again, that's all we have at the moment (10/95).
687 * Returns NULL from non-ranki systems.
688 * @todo Separate documentation for these linsolqr functions?
689 */
690
691 extern int32 linsolqr_org_row_to_org_col(linsolqr_system_t sys,
692 int32 org_row);
693 /**< See discussion under linsolqr_org_col_to_org_row(). */
694 extern int32 linsolqr_org_col_to_org_row(linsolqr_system_t sys,
695 int32 org_col);
696 /**<
697 * Pivoted original columns and pivoted original rows can be associated
698 * with one another via the pivot sequence. These functions returned the
699 * org_col/org_row associated with the given org_row/org_col. If the given
700 * org_row/org_col is not pivoted, a meaningless value is returned. The
701 * system must be previously factored. If not factored, these functions
702 * will return a value, but linsolqr may reorder making the value wrong.
703 * @todo Separate documentation for these linsolqr functions?
704 */
705
706 ASC_DLLSPEC void linsolqr_calc_row_dependencies(linsolqr_system_t sys);
707 /**< See discussion under linsolqr_calc_col_dependencies(). */
708
709 ASC_DLLSPEC void linsolqr_calc_col_dependencies(linsolqr_system_t sys);
710 /**<
711 * Given a factored system for which dependencies are not yet
712 * calculated, calculates row/column dependencies. This must be
713 * called before either linsolqr_org_row/col_dependency
714 * or linsolqr_row/col_dependence_coefs can be called.
715 * All rows/columns in the region specified to reorder
716 * and not part of the final factorization have their dependencies
717 * calculated.
718 *
719 * The calculation method is fairly standard.
720 * We will give the version for LU column dependency. Transpose it
721 * for row dependency. Similar things can be done for QR.
722 * Given a matrix region A, factor it and arrive at
723 * A11 | A12 where A11 = L U and columns of A12 are dependent,
724 * --------- that is A12 = A11 . Y. (The row A2 need not exist,
725 * A21 | eps but if it does, A22 must be small or we could have
726 * factored further.) Solve for Y. The Yi are the coefficients of the
727 * linear combination of columns in A11 which sum to A12.
728 * And, note that since the columns of A12 were treated during
729 * factoring as if they were ultimately going to be pivoted in,
730 * we only need to substitute on the other triangle factor, half the
731 * work of a regular solve.
732 * @todo Separate documentation for these linsolqr functions?
733 */
734
735 ASC_DLLSPEC mtx_sparse_t *linsolqr_row_dependence_coefs(
736 linsolqr_system_t sys, int32 orgrow);
737 /**< See discussion under linsolqr_col_dependence_coefs(). */
738
739 ASC_DLLSPEC mtx_sparse_t *linsolqr_col_dependence_coefs(
740 linsolqr_system_t sys, int32 orgcol
741 );
742 /**<
743 * Given an org row/col and a factored system, returns a mtx_sparse_t
744 * containing the org row/col indices and linear combination coefficients
745 * contributing to the given dependent row/col. If the orgrow/orgcol
746 * given is not dependent or dependencies have not been calculated,
747 * return is NULL.
748 * E.g. rowcoefs->idata[i] is the org row index of a contributing row,
749 * rowcoefs->data[i] is the dependency coefficient,
750 * for i = 0 to rowcoefs->len-1.
751 * Numeric dependency calculation is a numeric process with lots of room
752 * for interpretation of the results. Not everything that claims to
753 * contribute to a singularity really does so. We leave this interpre-
754 * tation to the the caller.
755 *
756 * The CALLER is responsible for deallocating the mtx_sparse_t returned;
757 * linsol wants nothing further to do with it.
758 *
759 * @todo Separate documentation for these linsolqr functions?
760 */
761
762 extern real64 linsolqr_org_row_dependency(linsolqr_system_t sys,
763 int32 dep, int32 ind);
764 /**< See discussion under linsolqr_org_col_dependency(). */
765 extern real64 linsolqr_org_col_dependency(linsolqr_system_t sys,
766 int32 dep, int32 ind);
767 /**<
768 * Any original row / column of the coefficient matrix which when submitted
769 * to the linear solver is not pivoted, is called dependent and can be
770 * written as a linear combination of the independent (pivoted) original
771 * rows / columns. These functions return the previously computed
772 * coefficients of the linear combination. The system must be previously
773 * factored and the independent row / column must be a member of the
774 * set of row / column pivots obtained by linsolqr_get_pivot_sets.
775 * This is a slow and clunky way to retrieve dependency info.
776 * This ought to be done away with when the above function is done.
777 *
778 * @todo Separate documentation for these linsolqr functions?
779 */
780
781 ASC_DLLSPEC int linsolqr_solve(linsolqr_system_t sys, real64 *rhs);
782 /**<
783 * Solves the system of linear equations (if necessary) utilizing the
784 * specified rhs along with the previously factored matrix. The rhs
785 * is automatically checked if the matrix factors need to be transposed
786 * or not (see linsolqr_add_rhs.)
787 * Solution method will be appropriate to the factorization method used
788 * in linsolqr_factor.
789 * @return 0 if ok, 1 if not.
790 */
791
792 extern real64 linsolqr_var_value(linsolqr_system_t sys,
793 real64 *rhs, int32 var);
794 /**<
795 * Returns the value of the variable in the solution vector associated
796 * with the given rhs and either the matrix or its transpose. The rhs
797 * must be solved for first. If rhs specifies transpose, then var is
798 * expected to be an original row number, otherwise it should be an
799 * original column number.
800 * If sys, rhs, and var are not proper, the value returned is 0.0
801 * and a warning is issued to stderr.
802 */
803
804 ASC_DLLSPEC boolean linsolqr_copy_solution(linsolqr_system_t sys,
805 real64 *rhs, real64 *vector);
806 /**<
807 * Once a sys has been factored and rhs solved, this
808 * fills in a copy of the solution vector associated with rhs. Caller
809 * must provide vector and vector must be of length at least as long
810 * as the order of the matrix that was solved. The vector entries
811 * will be in org_col order if the rhs is normal or in org_row order
812 * if the rhs is a transpose.
813 *
814 * @return TRUE if something is amiss, FALSE otherwise.
815 *
816 */
817
818 extern real64 linsolqr_eqn_residual(linsolqr_system_t sys,
819 real64 *rhs, int32 eqn);
820 /**<
821 * Returns the equation residual using the solution vector associated
822 * with the given rhs and either the matrix or its transpose.
823 * <pre>
824 * T
825 * residual = A x - b or residual = A x - b
826 * </pre>
827 * The rhs must be solved for first. If rhs specifies transpose, then
828 * eqn is expected to be an original column number, otherwise it should
829 * be an original row number.
830 * If the system and rhs and eqn are not properly solved, the return
831 * value is MAXDOUBLE.
832 */
833
834 extern boolean linsolqr_calc_residual(linsolqr_system_t sys,
835 real64 *rhs, real64 *vector);
836 /**<
837 * Returns the residuals using the solution associated
838 * with the given rhs and either the matrix or its transpose.
839 * <pre>
840 * T
841 * residual = A x - b or residual = A x - b
842 * </pre>
843 * The rhs must be solved for first. Caller
844 * must provide vector and vector must be of length at least as long
845 * as the order of the matrix that was solved. The vector entries
846 * will be in org_row order if the rhs is normal or in org_col order
847 * if the rhs is a transpose. Entries of the vector which do not
848 * correspond to rows/cols of factored system solved will not be modified.<br><br>
849 *
850 * @return If the system and rhs are not properly solved, or other is amiss
851 * return value is TRUE, else FALSE.
852 */
853
854 /* miscellaneous functions for C programmers wanting to know things. */
855
856 extern size_t linsolqr_size(linsolqr_system_t sys);
857 /**<
858 * Returns the amount of memory in use by a system and all its
859 * bits and pieces. User supplied RHS vectors and coefficient
860 * matrix are NOT included in the size calculation. The user
861 * must do accounting for those.
862 */
863
864 extern void linsolqr_free_reused_mem(void);
865 /**<
866 * Deallocates any memory that linsolqr may be squirrelling away for
867 * internal reuse. Calling this while any slv_system_t using linsolqr exists
868 * is likely to be fatal: handle with care.
869 * There isn't a way to query how many bytes this is.
870 */
871
872 /*-----------------------------------------------------------------------------
873 * The following calls exist to facilitate debugging of the linear
874 * solver when it is being tested on large systems. Do not use them
875 * in routine coding. If you need access to the factor/inverse matrices
876 * for computational purpose, you should probably consider adding that
877 * functionality to linsol.
878 */
879
880 ASC_DLLSPEC mtx_matrix_t linsolqr_get_factors(linsolqr_system_t sys);
881 /**< See discussion under linsolqr_get_inverse(). */
882 extern mtx_matrix_t linsolqr_get_inverse(linsolqr_system_t sys);
883 /**<
884 * Returns the handle of the factor (L\U, Q\R) or inverse matrix.
885 * The handle may be NULL, and should be checked before use.
886 * The matrix should not be tampered (other than to look at it)
887 * with if any other mathematical operations will take place with
888 * sys. Linsol expects to deallocate both factors and inverse; do
889 * not do so yourself.
890 *
891 * @NOTE All the factorization methods use some sort of dense vectors
892 * for storing pivots and whatnot. As with a FORTRAN factorization,
893 * there's not much you can do with a factorization without
894 * also understanding in detail the factorization routine and getting
895 * your hands on the dense vectors.
896 *
897 * @todo Separate documentation of linsolqr_get_factors()
898 * and linsolqr_get_inverse()?
899 */
900
901 extern mtx_region_t *linsolqr_get_region(linsolqr_system_t sys);
902 /**<
903 * Returns a pointer to the current linsolqr region.
904 * This is being created for use in the check numeric dependency
905 * routine and may be removed once the new "problem manager" window
906 * is created to take over this functionality (and others).
907 */
908
909 extern int linsolqr_setup_ngslv(linsolqr_system_t sys,
910 real64 *rhs,
911 mtx_range_t *un_p_rng,
912 real64 *tmpvec);
913 /**< Sets up the NGSlv solver. */
914 extern real64 *linsolqr_get_varvalue(linsolqr_system_t sys, int n);
915 /**< Returns the value of the nth variable in sys. */
916
917 /** @} */
918
919 #endif /* ASC_LINSOLQR_H */

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22