/[ascend]/trunk/ascend4/solver/linsolqr.h
ViewVC logotype

Contents of /trunk/ascend4/solver/linsolqr.h

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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