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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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