Parent Directory | Revision Log

Revision **1014** -
(**show annotations**)
(**download**)
(**as text**)

*Wed Jan 3 04:49:17 2007 UTC*
(16 years, 5 months ago)
by *johnpye*

File MIME type: text/x-chdr

File size: 39624 byte(s)

File MIME type: text/x-chdr

File size: 39624 byte(s)

Fixing @addtogroup comments

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 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 | * linsol II: 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 | /** @addtogroup linear |

216 | Linear solver routines |

217 | @{ */ |

218 | |

219 | #define LINSOLMTX_DEBUG FALSE |

220 | /**< Debug mode. */ |

221 | #define LINQR_DROP_TOLERANCE 1.0e-16 |

222 | /**< This is the default for drop tolerances in methods which use a drop tol */ |

223 | |

224 | typedef struct linsolqr_header *linsolqr_system_t; |

225 | /**< linsolqr_system_t is the linear system handle. */ |

226 | |

227 | enum factor_class { |

228 | unknown_c, /**< error handling class */ |

229 | ranki = 100, /**< all ranki (and gauss until broken out) methods */ |

230 | s_qr = 200 /**< all sparse qr methods */ |

231 | }; |

232 | |

233 | enum reorder_method { |

234 | unknown_r, /**< error handling method */ |

235 | natural = 1000, /**< do nothing reorder */ |

236 | spk1 = 2000, /**< Stadtherr's SPK1 reordering */ |

237 | tspk1 = 3000 /**< transpose of Stadtherr's SPK1 reordering good for gauss */ |

238 | /* future work: |

239 | invspk1, spk1 then diagonally inverted |

240 | invtspk1, tspk1 then diagonally inverted |

241 | widespk1, spk1 of an MxN region, N > M |

242 | */ |

243 | }; |

244 | |

245 | enum factor_method { |

246 | unknown_f = 0, /**< error handling method */ |

247 | ranki_kw = 1, /**< original linsol method */ |

248 | ranki_jz = 2, /**< original linsol method with pseudo-complete pivoting */ |

249 | ranki_kw2 = 3, /**< better data structure version of ranki_kw, w/drop tol */ |

250 | ranki_jz2 = 4, /**< better data structure version of ranki_jz, w/drop tol */ |

251 | ranki_ka = 12, /**< kirks hacked verion of the basic method. */ |

252 | plain_qr = 5, /**< rectangular col pivoted qr variants */ |

253 | cond_qr = 6, /**< stewart-based sparse QR method new Coke */ |

254 | ranki_ba2 = 7, /**< proper linked list implementation of ranki (dragfree) */ |

255 | opt_qr = 8, /**< coming soon */ |

256 | ls_qr = 9, /**< anticipated */ |

257 | gauss_ba2 = 10, /**< anticipated */ |

258 | symmetric_lu = 11 /**< anticipated */ |

259 | }; |

260 | |

261 | ASC_DLLSPEC(int ) g_linsolqr_timing; |

262 | /**< |

263 | * If nonzero and certain internal defines are set, factorization |

264 | * generates a fill and timing message. |

265 | */ |

266 | |

267 | /* KAA_DEBUG */ |

268 | /* unimplemented function |

269 | extern boolean linsolqr_col_is_a_spike (mtx_matrix_t mtx, int32 col); |

270 | */ |

271 | /* Returns true if the column in question is a spike. */ |

272 | |

273 | |

274 | /* Functions for managing interfaces */ |

275 | |

276 | extern char *linsolqr_rmethods(void); |

277 | /**< |

278 | * Returns a comma-separated list of the names of reordering methods |

279 | * implemented. If you implement a new method, update this |

280 | * function. Do not free the pointer returned.<br><br> |

281 | * |

282 | * Not all reorderings are appropriate to all factorizations. |

283 | */ |

284 | |

285 | ASC_DLLSPEC(char *) linsolqr_fmethods(void); |

286 | /**< |

287 | * Returns a comma-separated list of the names of factorization methods |

288 | * implemented. If you implement a new method, update this |

289 | * function. Do not free the pointer returned. The string returned |

290 | * will contain no whitespace of any sort, tabs, blanks, or \n. |

291 | */ |

292 | |

293 | extern enum reorder_method linsolqr_rmethod_to_enum(char *s); |

294 | /**< |

295 | * Returns the enum of a reorder method with the name s. |

296 | * If you implement a new method, update this function. |

297 | */ |

298 | |

299 | extern enum factor_method linsolqr_fmethod_to_enum(char *s); |

300 | /**< |

301 | * Returns the enum of a factor method with the name s. |

302 | * If you implement a new method, update this function. |

303 | */ |

304 | |

305 | extern enum factor_class linsolqr_fmethod_to_fclass(enum factor_method fm); |

306 | /**< |

307 | * Returns the enum of the factor class containing the method given. |

308 | * If you implement a new method or class, update this function. |

309 | */ |

310 | |

311 | ASC_DLLSPEC(char *) linsolqr_enum_to_rmethod(enum reorder_method m); |

312 | /**< |

313 | * <!-- s=linsolqr_enum_to_rmethod(m); --> |

314 | * |

315 | * Returns the name of a reorder method with the enum m. |

316 | * If you implement a new method, update this function. |

317 | * Do not free the name. |

318 | */ |

319 | |

320 | extern char *linsolqr_enum_to_fmethod(enum factor_method m); |

321 | /**< |

322 | * Returns the name of a factor method with the enum m. |

323 | * If you implement a new method, update this function. |

324 | * Do not free the name. |

325 | */ |

326 | |

327 | ASC_DLLSPEC(char *) linsolqr_rmethod_description(enum reorder_method meth); |

328 | /**< |

329 | * Returns a string describing the method inquired on. Do not mess |

330 | * with the string: linsolqr owns it. |

331 | * If you implement a new method, update this function. |

332 | */ |

333 | |

334 | ASC_DLLSPEC(char *) linsolqr_fmethod_description(enum factor_method meth); |

335 | /**< |

336 | * Returns a string describing the method inquired on. Do not mess |

337 | * with the string: linsolqr owns it. |

338 | * If you implement a new method, update this function. |

339 | */ |

340 | |

341 | /* Functions for specifying problems and controlling them */ |

342 | |

343 | extern linsolqr_system_t linsolqr_create(void); |

344 | /**< |

345 | * Creates a linear system and returns a pointer to it. Initially the |

346 | * system has no coefficient matrix and no rhs. |

347 | */ |

348 | |

349 | extern void linsolqr_destroy(linsolqr_system_t sys); |

350 | /**< |

351 | * Destroys the linear system. The coefficient matrix and each rhs are |

352 | * not destroyed by this call. |

353 | */ |

354 | |

355 | extern void linsolqr_set_matrix(linsolqr_system_t sys, mtx_matrix_t mtx); |

356 | /**< |

357 | * Sets the coefficient matrix to mtx. |

358 | */ |

359 | |

360 | ASC_DLLSPEC(void ) linsolqr_set_region(linsolqr_system_t sys, mtx_region_t region); |

361 | /**< |

362 | * Sets the reg to region. |

363 | */ |

364 | |

365 | ASC_DLLSPEC(mtx_matrix_t ) linsolqr_get_matrix(linsolqr_system_t sys); |

366 | /**< |

367 | * Returns the coefficient matrix. |

368 | */ |

369 | |

370 | ASC_DLLSPEC(void ) linsolqr_add_rhs(linsolqr_system_t sys, |

371 | real64 *rhs, |

372 | boolean transpose); |

373 | /**< |

374 | * Adds the given rhs to the collection of rhs's already part of the |

375 | * system. |

376 | * |

377 | * You continue to be responsible for deallocating rhs. We do not take |

378 | * over management of that memory. Do not free the rhs you have given |

379 | * us without calling linsolqr_remove_rhs first.<br><br> |

380 | * |

381 | * Rhs should point to an array of reals indexed by original |

382 | * column number if the linear system is to be solved using the transpose |

383 | * of the matrix or by original row number if the matrix is not to be |

384 | * transposed. This is determined using the boolean transpose. The |

385 | * rhs should be refered to in the future by its pointer. |

386 | * <pre> |

387 | * For the mathematically impaired: |

388 | * M' denotes Transpose(M). |

389 | * |

390 | * transpose==FALSE --> rhs will be solved for x in A x = b. |

391 | * --> rhs b is an org row indexed column vector. |

392 | * --> x is an org col indexed column vector. |

393 | * --> yields dx = Newton direction in J.dx = -f |

394 | * if the b given is -f. |

395 | * |

396 | * transpose==TRUE --> rhs will be solved for y in A' y = d. |

397 | * --> rhs d is an org column indexed column vector. |

398 | * --> y is an org row indexed column vector. |

399 | * --> yields dx = Newton direction in J'.dx = -f. |

400 | * if the d given is -f and f,dx are properly indexed. |

401 | * This may be useful if A was created using |

402 | * mtx_transpose in order to improve factorization. |

403 | * |

404 | * Useful matrix identity: (M')^-1 == (M^-1)' for invertible M. |

405 | * </pre> |

406 | */ |

407 | |

408 | ASC_DLLSPEC(void ) linsolqr_remove_rhs(linsolqr_system_t sys, real64 *rhs); |

409 | /**< |

410 | * Removes the given rhs from the system. The rhs is not destroyed, just |

411 | * removed from the system. |

412 | */ |

413 | |

414 | extern int32 linsolqr_number_of_rhs(linsolqr_system_t sys); |

415 | /**< |

416 | * Returns the number of rhs's currently part of the system. |

417 | */ |

418 | |

419 | ASC_DLLSPEC(real64 *) linsolqr_get_rhs(linsolqr_system_t sys, int n); |

420 | /**< |

421 | * Returns the n-th rhs, where rhs's are indexed in the order they were |

422 | * added using linsolqr_add_rhs() from 0 to (# rhs's)-1. NULL is returned |

423 | * if the index is out of range. |

424 | */ |

425 | |

426 | ASC_DLLSPEC(void ) linsolqr_matrix_was_changed(linsolqr_system_t sys); |

427 | /**< |

428 | * Informs the solver that a numerical value of a non-zero was changed. |

429 | * This must be called whenever any numerical changes to the matrix are |

430 | * made. |

431 | */ |

432 | |

433 | ASC_DLLSPEC(void ) linsolqr_rhs_was_changed(linsolqr_system_t sys, |

434 | real64 *rhs); |

435 | /**< |

436 | * Informs the solver that the given rhs has been modified. This must be |

437 | * called whenever the rhs is modified. |

438 | */ |

439 | |

440 | extern void linsolqr_set_pivot_zero(linsolqr_system_t sys, |

441 | real64 pivot_zero); |

442 | /**< |

443 | * Sets the pivot zero for the system. Pivots less than or equal to |

444 | * this value are regarded as zero. linsolqr_set_pivot_zero() will |

445 | * automatically call linsolqr_matrix_was_changed(). |

446 | */ |

447 | extern real64 linsolqr_pivot_zero(linsolqr_system_t sys); |

448 | /**< |

449 | * Gets the pivot zero for the system. Pivots less than or equal to |

450 | * this value are regarded as zero. |

451 | */ |

452 | |

453 | extern void linsolqr_set_pivot_tolerance(linsolqr_system_t sys, real64 ptol); |

454 | /**< See discussion under linsolqr_drop_tolerance(). */ |

455 | extern real64 linsolqr_pivot_tolerance(linsolqr_system_t sys); |

456 | /**< See discussion under linsolqr_drop_tolerance(). */ |

457 | extern void linsolqr_set_condition_tolerance(linsolqr_system_t sys, real64 ctol); |

458 | /**< See discussion under linsolqr_drop_tolerance(). */ |

459 | extern real64 linsolqr_condition_tolerance(linsolqr_system_t sys); |

460 | /**< See discussion under linsolqr_drop_tolerance(). */ |

461 | extern void linsolqr_set_drop_tolerance(linsolqr_system_t sys, real64 dtol); |

462 | /**< See discussion under linsolqr_drop_tolerance(). */ |

463 | extern real64 linsolqr_drop_tolerance(linsolqr_system_t sys); |

464 | /**< |

465 | * <pre> |

466 | * linsolqr_set_pivot_tolerance(sys,ptol) |

467 | * ptol = linsolqr_pivot_tolerance(sys) |

468 | * |

469 | * linsolqr_set_condition_tolerance(sys,ctol) |

470 | * ctol = linsolqr_condition_tolerance(sys) |

471 | * |

472 | * linsolqr_set_drop_tolerance(sys,dtol) |

473 | * dtol = linsolqr_drop_tolerance(sys) |

474 | * |

475 | * linsolqr_system_t sys; |

476 | * real64 ptol, ctol,dtol; |

477 | * |

478 | * Sets/gets the pivot/condition/drop tolerance for the system. Semantics of |

479 | * tolerances vary with method. Not to be confused with pivot zero, epsilon. |

480 | * |

481 | * pivot_tolerance: Pivots less than this fraction of the maximum pivot |

482 | * value in the same row or column (depending on method) are disregarded. |

483 | * |

484 | * ranki_kw: pivot_tolerance applies to row. condition_tolerance ignored. |

485 | * ranki_jz: pivot_tolerance applies to row and col. condition ignored. |

486 | * ranki_ba2: |

487 | * ranki_kw2: |

488 | * ranki_jz2: as ranki_kw/ranki_jz except that matrix entries below |

489 | * dtol in magnitude are dropped. |

490 | * plain_qr: pivot_tolerance applies to cols. condition_tolerance used |

491 | * with a condition heuristic rather than actual condition. |

492 | * |

493 | * cond_qr and variants: |

494 | * Unpivoted columns are ordered by fill potential and then the first |

495 | * of these to meet the criterion CI*condition_tolerance <= min_CI |

496 | * is chosen as the next column to pivot with. |

497 | * Pivot_tolerance is applied to choose the pivot within the |

498 | * selected column. |

499 | * |

500 | * linsolqr_set_pivot/condition/drop_tolerance() will automatically call |

501 | * linsolqr_matrix_was_changed(). |

502 | * </pre> |

503 | * @todo Separate documentation for these linsolqr functions? |

504 | */ |

505 | |

506 | /* Functions for analyzing and querying linear systems. */ |

507 | |

508 | extern enum factor_class linsolqr_fclass(linsolqr_system_t sys); |

509 | /**< |

510 | * Returns the most recently set factorization class of the system. |

511 | * The system should be previously prepped. |

512 | */ |

513 | |

514 | ASC_DLLSPEC(enum factor_method ) linsolqr_fmethod(linsolqr_system_t sys); |

515 | /**< |

516 | * Returns the most recently used factorization method of the system. |

517 | * The system should be previously prepped. |

518 | */ |

519 | |

520 | ASC_DLLSPEC(enum reorder_method ) linsolqr_rmethod(linsolqr_system_t sys); |

521 | /**< |

522 | * Returns the most recently set reorder method of the system. |

523 | * The system should be previously prepped. |

524 | */ |

525 | |

526 | extern int32 linsolqr_rank(linsolqr_system_t sys); |

527 | /**< |

528 | * Returns the rank of the system. The system must be previously |

529 | * factored. |

530 | */ |

531 | |

532 | extern real64 linsolqr_smallest_pivot(linsolqr_system_t sys); |

533 | /**< |

534 | * Returns the smallest pivot accepted in solving the system. The |

535 | * system must be previously factored. If no pivoting was performed, |

536 | * MAXDOUBLE is returned. |

537 | */ |

538 | |

539 | extern int linsolqr_prep(linsolqr_system_t sys, |

540 | enum factor_class fclass); |

541 | /**< |

542 | * This function is analogous to slv_select_solver. It |

543 | * takes care of allocations. The prep call should be done (once) at |

544 | * the beginning of any series of linear solutions being performed on |

545 | * on the same linsolqr system with the same factorization class.<br><br> |

546 | * |

547 | * You must prep before doing a reordering, factorization or solution |

548 | * as prep sets up the appropriate internal settings. You should set |

549 | * the matrix before you call prep. |

550 | * Prep (or subsequent preps) do not affect the right hand sides of |

551 | * the system.<br><br> |

552 | * |

553 | * If you wish to change from 1 factoring method to another and they |

554 | * are not part of the class of compatible methods, you should call |

555 | * prep again with the new class. |

556 | * After prep is called, reorder should be called before factorization. |

557 | * As most of the sparse routines depend for performance on the |

558 | * prior reordering.<br><br> |

559 | * |

560 | * Return 0 if ok, or 1 otherwise. |

561 | */ |

562 | |

563 | ASC_DLLSPEC(int ) linsolqr_reorder(linsolqr_system_t sys, |

564 | mtx_region_t *region, |

565 | enum reorder_method rmeth); |

566 | /**< |

567 | * The specified region of the coefficient matrix is reordered. This |

568 | * must be called before factoring the matrix. The specified region |

569 | * is assumed to contain only nonempty rows and columns and have a full |

570 | * diagonal. |

571 | * If the coefficient matrix in use is from a nonlinear system, the |

572 | * pattern in the coefficient matrix should be the structural one (as |

573 | * opposed to the numerically derived incidence which may be less.)<br><br> |

574 | * |

575 | * If you use the numerically derived incidence, you will need to reorder |

576 | * before every factorization. This is generally not cost effective. |

577 | * If region given is mtx_ENTIRE_MATRIX, a search will be done to find |

578 | * an appropriate bounding region in the coefficient mtx. This is |

579 | * not a particularly cheap search.<br><br> |

580 | * |

581 | * This function is analogous to slv_presolve. It |

582 | * takes care of any structural analysis necessary |

583 | * to the linear solution. The reorder call should be done (once) at |

584 | * the beginning of any series of linear solutions being performed on |

585 | * on structurally identical matrices.<br><br> |

586 | * |

587 | * The reordering done must be appropriate to the factorization class. |

588 | * You must reorder before doing a factorization or solution as reorder |

589 | * sets up the appropriate internal settings. Even were the internals |

590 | * method independent, the reordering is critical to the performance of |

591 | * these methods.<br><br> |

592 | * Return 0 if ok, 1 if not. |

593 | * <pre> |

594 | * Brief notes on the reorderings available. |

595 | * SPK1: The region you give becomes the region for the problem, |

596 | * TSPK1: but for factorization purposes rows and columns that |

597 | * do not intersect the main diagonal within the region |

598 | * are considered structurally (therefore numerically) dependent. |

599 | * Natural: Blesses the system and does nothing. |

600 | * Again, the rows/cols not in the diagonal are dependent. |

601 | * </pre> |

602 | */ |

603 | |

604 | ASC_DLLSPEC(int ) linsolqr_factor(linsolqr_system_t sys, |

605 | enum factor_method fmethod); |

606 | /**< |

607 | * Decompose the reordered region of a copy of the coefficient matrix |

608 | * into upper and lower triangular factors (if necessary) which can be |

609 | * inverted easily when applied to any rhs. Matrix must be factored in |

610 | * order to perform any of the operations below. The numerical rank and |

611 | * the smallest pivot encountered during pivoting are computed in the |

612 | * process. Factorization method used is that given if it is compatible |

613 | * with the class set when the prep was done, otherwise the call fails.<br><br> |

614 | * |

615 | * If you are handed a linsolqr system and want to refactor it using the |

616 | * usual method, but don't know what that method is, call like: |

617 | * status = linsolqr_factor(sys,linsolqr_fmethod(sys));<br><br> |

618 | * |

619 | * Return 0 if ok, 1 if not. |

620 | */ |

621 | |

622 | extern int linsolqr_get_pivot_sets(linsolqr_system_t sys, |

623 | unsigned *org_rowpivots, |

624 | unsigned *org_colpivots); |

625 | /**< |

626 | * Returns the set of original row numbers / original column numbers which |

627 | * have been pivoted. org_rowpivots and org_colpivots are assumed to be |

628 | * sets created by (or at least for) the set module with sufficient size |

629 | * before calling this function. They must also be previously NULLed. |

630 | * The system must be previously factored. |

631 | * The sets input should be the result of set_create(neqn),set_create(nvar). |

632 | * There is no association of rows with columns here.<br><br> |

633 | * |

634 | * @return Status is 0 if not factored, 1 if factored. If 0, sets will not be |

635 | * changed. |

636 | * This bizarre piece of functionality should be done away with as soon |

637 | * as the equivalents below have been implemented. |

638 | */ |

639 | |

640 | ASC_DLLSPEC(mtx_sparse_t *) linsolqr_unpivoted_rows(linsolqr_system_t sys); |

641 | /**< See discussion under linsolqr_unpivoted_cols(). */ |

642 | ASC_DLLSPEC(mtx_sparse_t *) linsolqr_unpivoted_cols(linsolqr_system_t sys); |

643 | /**< |

644 | * Returns the set of original row numbers / original column numbers which |

645 | * have NOT been pivoted and the rejected pivots in an mtx_sparse_t. |

646 | * Return is NULL if sys not factored or if no unpivoted rows/cols exist. |

647 | * E.g. singrows->idata[i] are the original rows that were not pivoted, |

648 | * singrows->data[i] is the pivot that was rejected, |

649 | * for i = 0 to singrows->len-1. |

650 | * If len is 0, the data and idata pointers may be NULL. |

651 | * |

652 | * The CALLER is responsible for deallocating the mtx_sparse_t returned; |

653 | * linsolqr wants nothing further to do with it. |

654 | * @bug Only deals with ranki based square systems at the moment. |

655 | * Then again, that's all we have at the moment (10/95). |

656 | * Returns NULL from non-ranki systems. |

657 | * @todo Separate documentation for these linsolqr functions? |

658 | */ |

659 | |

660 | extern mtx_sparse_t *linsolqr_pivoted_rows(linsolqr_system_t sys); |

661 | /**< See discussion under linsolqr_pivoted_cols(). */ |

662 | extern mtx_sparse_t *linsolqr_pivoted_cols(linsolqr_system_t sys); |

663 | /**< |

664 | * Returns the set of original row numbers / original column numbers which |

665 | * have been pivoted and the pivots in an mtx_sparse_t. |

666 | * Return is NULL if sys not factored or if no pivoted rows/cols exist. |

667 | * E.g. pivrows->idata[i] are the original rows that were pivoted, |

668 | * pivrows->data[i] is the pivot, |

669 | * for i = 0 to pivrows->len-1. |

670 | * If len is 0, the data and idata pointers may be NULL. |

671 | * |

672 | * The CALLER is responsible for deallocating the mtx_sparse_t returned; |

673 | * linsolqr wants nothing further to do with it. |

674 | * @bug Only deals with ranki based square systems at the moment. |

675 | * Then again, that's all we have at the moment (10/95). |

676 | * Returns NULL from non-ranki systems. |

677 | * @todo Separate documentation for these linsolqr functions? |

678 | */ |

679 | |

680 | extern int32 linsolqr_org_row_to_org_col(linsolqr_system_t sys, |

681 | int32 org_row); |

682 | /**< See discussion under linsolqr_org_col_to_org_row(). */ |

683 | extern int32 linsolqr_org_col_to_org_row(linsolqr_system_t sys, |

684 | int32 org_col); |

685 | /**< |

686 | * Pivoted original columns and pivoted original rows can be associated |

687 | * with one another via the pivot sequence. These functions returned the |

688 | * org_col/org_row associated with the given org_row/org_col. If the given |

689 | * org_row/org_col is not pivoted, a meaningless value is returned. The |

690 | * system must be previously factored. If not factored, these functions |

691 | * will return a value, but linsolqr may reorder making the value wrong. |

692 | * @todo Separate documentation for these linsolqr functions? |

693 | */ |

694 | |

695 | ASC_DLLSPEC(void) linsolqr_calc_row_dependencies(linsolqr_system_t sys); |

696 | /**< See discussion under linsolqr_calc_col_dependencies(). */ |

697 | |

698 | ASC_DLLSPEC(void) linsolqr_calc_col_dependencies(linsolqr_system_t sys); |

699 | /**< |

700 | * Given a factored system for which dependencies are not yet |

701 | * calculated, calculates row/column dependencies. This must be |

702 | * called before either linsolqr_org_row/col_dependency |

703 | * or linsolqr_row/col_dependence_coefs can be called. |

704 | * All rows/columns in the region specified to reorder |

705 | * and not part of the final factorization have their dependencies |

706 | * calculated. |

707 | * |

708 | * The calculation method is fairly standard. |

709 | * We will give the version for LU column dependency. Transpose it |

710 | * for row dependency. Similar things can be done for QR. |

711 | * Given a matrix region A, factor it and arrive at |

712 | * A11 | A12 where A11 = L U and columns of A12 are dependent, |

713 | * --------- that is A12 = A11 . Y. (The row A2 need not exist, |

714 | * A21 | eps but if it does, A22 must be small or we could have |

715 | * factored further.) Solve for Y. The Yi are the coefficients of the |

716 | * linear combination of columns in A11 which sum to A12. |

717 | * And, note that since the columns of A12 were treated during |

718 | * factoring as if they were ultimately going to be pivoted in, |

719 | * we only need to substitute on the other triangle factor, half the |

720 | * work of a regular solve. |

721 | * @todo Separate documentation for these linsolqr functions? |

722 | */ |

723 | |

724 | ASC_DLLSPEC(mtx_sparse_t *) linsolqr_row_dependence_coefs( |

725 | linsolqr_system_t sys, int32 orgrow); |

726 | /**< See discussion under linsolqr_col_dependence_coefs(). */ |

727 | |

728 | ASC_DLLSPEC(mtx_sparse_t *) linsolqr_col_dependence_coefs( |

729 | linsolqr_system_t sys, int32 orgcol |

730 | ); |

731 | /**< |

732 | * Given an org row/col and a factored system, returns a mtx_sparse_t |

733 | * containing the org row/col indices and linear combination coefficients |

734 | * contributing to the given dependent row/col. If the orgrow/orgcol |

735 | * given is not dependent or dependencies have not been calculated, |

736 | * return is NULL. |

737 | * E.g. rowcoefs->idata[i] is the org row index of a contributing row, |

738 | * rowcoefs->data[i] is the dependency coefficient, |

739 | * for i = 0 to rowcoefs->len-1. |

740 | * Numeric dependency calculation is a numeric process with lots of room |

741 | * for interpretation of the results. Not everything that claims to |

742 | * contribute to a singularity really does so. We leave this interpre- |

743 | * tation to the the caller. |

744 | * |

745 | * The CALLER is responsible for deallocating the mtx_sparse_t returned; |

746 | * linsol wants nothing further to do with it. |

747 | * |

748 | * @todo Separate documentation for these linsolqr functions? |

749 | */ |

750 | |

751 | extern real64 linsolqr_org_row_dependency(linsolqr_system_t sys, |

752 | int32 dep, int32 ind); |

753 | /**< See discussion under linsolqr_org_col_dependency(). */ |

754 | extern real64 linsolqr_org_col_dependency(linsolqr_system_t sys, |

755 | int32 dep, int32 ind); |

756 | /**< |

757 | * Any original row / column of the coefficient matrix which when submitted |

758 | * to the linear solver is not pivoted, is called dependent and can be |

759 | * written as a linear combination of the independent (pivoted) original |

760 | * rows / columns. These functions return the previously computed |

761 | * coefficients of the linear combination. The system must be previously |

762 | * factored and the independent row / column must be a member of the |

763 | * set of row / column pivots obtained by linsolqr_get_pivot_sets. |

764 | * This is a slow and clunky way to retrieve dependency info. |

765 | * This ought to be done away with when the above function is done. |

766 | * |

767 | * @todo Separate documentation for these linsolqr functions? |

768 | */ |

769 | |

770 | ASC_DLLSPEC(int ) linsolqr_solve(linsolqr_system_t sys, real64 *rhs); |

771 | /**< |

772 | * Solves the system of linear equations (if necessary) utilizing the |

773 | * specified rhs along with the previously factored matrix. The rhs |

774 | * is automatically checked if the matrix factors need to be transposed |

775 | * or not (see linsolqr_add_rhs.) |

776 | * Solution method will be appropriate to the factorization method used |

777 | * in linsolqr_factor. |

778 | * @return 0 if ok, 1 if not. |

779 | */ |

780 | |

781 | extern real64 linsolqr_var_value(linsolqr_system_t sys, |

782 | real64 *rhs, int32 var); |

783 | /**< |

784 | * Returns the value of the variable in the solution vector associated |

785 | * with the given rhs and either the matrix or its transpose. The rhs |

786 | * must be solved for first. If rhs specifies transpose, then var is |

787 | * expected to be an original row number, otherwise it should be an |

788 | * original column number. |

789 | * If sys, rhs, and var are not proper, the value returned is 0.0 |

790 | * and a warning is issued to stderr. |

791 | */ |

792 | |

793 | ASC_DLLSPEC(boolean ) linsolqr_copy_solution(linsolqr_system_t sys, |

794 | real64 *rhs, real64 *vector); |

795 | /**< |

796 | * Once a sys has been factored and rhs solved, this |

797 | * fills in a copy of the solution vector associated with rhs. Caller |

798 | * must provide vector and vector must be of length at least as long |

799 | * as the order of the matrix that was solved. The vector entries |

800 | * will be in org_col order if the rhs is normal or in org_row order |

801 | * if the rhs is a transpose. |

802 | * |

803 | * @return TRUE if something is amiss, FALSE otherwise. |

804 | * |

805 | */ |

806 | |

807 | extern real64 linsolqr_eqn_residual(linsolqr_system_t sys, |

808 | real64 *rhs, int32 eqn); |

809 | /**< |

810 | * Returns the equation residual using the solution vector associated |

811 | * with the given rhs and either the matrix or its transpose. |

812 | * <pre> |

813 | * T |

814 | * residual = A x - b or residual = A x - b |

815 | * </pre> |

816 | * The rhs must be solved for first. If rhs specifies transpose, then |

817 | * eqn is expected to be an original column number, otherwise it should |

818 | * be an original row number. |

819 | * If the system and rhs and eqn are not properly solved, the return |

820 | * value is MAXDOUBLE. |

821 | */ |

822 | |

823 | extern boolean linsolqr_calc_residual(linsolqr_system_t sys, |

824 | real64 *rhs, real64 *vector); |

825 | /**< |

826 | * Returns the residuals using the solution associated |

827 | * with the given rhs and either the matrix or its transpose. |

828 | * <pre> |

829 | * T |

830 | * residual = A x - b or residual = A x - b |

831 | * </pre> |

832 | * The rhs must be solved for first. Caller |

833 | * must provide vector and vector must be of length at least as long |

834 | * as the order of the matrix that was solved. The vector entries |

835 | * will be in org_row order if the rhs is normal or in org_col order |

836 | * if the rhs is a transpose. Entries of the vector which do not |

837 | * correspond to rows/cols of factored system solved will not be modified.<br><br> |

838 | * |

839 | * @return If the system and rhs are not properly solved, or other is amiss |

840 | * return value is TRUE, else FALSE. |

841 | */ |

842 | |

843 | /* miscellaneous functions for C programmers wanting to know things. */ |

844 | |

845 | extern size_t linsolqr_size(linsolqr_system_t sys); |

846 | /**< |

847 | * Returns the amount of memory in use by a system and all its |

848 | * bits and pieces. User supplied RHS vectors and coefficient |

849 | * matrix are NOT included in the size calculation. The user |

850 | * must do accounting for those. |

851 | */ |

852 | |

853 | extern void linsolqr_free_reused_mem(void); |

854 | /**< |

855 | * Deallocates any memory that linsolqr may be squirrelling away for |

856 | * internal reuse. Calling this while any slv_system_t using linsolqr exists |

857 | * is likely to be fatal: handle with care. |

858 | * There isn't a way to query how many bytes this is. |

859 | */ |

860 | |

861 | /*----------------------------------------------------------------------------- |

862 | * The following calls exist to facilitate debugging of the linear |

863 | * solver when it is being tested on large systems. Do not use them |

864 | * in routine coding. If you need access to the factor/inverse matrices |

865 | * for computational purpose, you should probably consider adding that |

866 | * functionality to linsol. |

867 | */ |

868 | |

869 | ASC_DLLSPEC(mtx_matrix_t ) linsolqr_get_factors(linsolqr_system_t sys); |

870 | /**< See discussion under linsolqr_get_inverse(). */ |

871 | extern mtx_matrix_t linsolqr_get_inverse(linsolqr_system_t sys); |

872 | /**< |

873 | * Returns the handle of the factor (L\U, Q\R) or inverse matrix. |

874 | * The handle may be NULL, and should be checked before use. |

875 | * The matrix should not be tampered (other than to look at it) |

876 | * with if any other mathematical operations will take place with |

877 | * sys. Linsol expects to deallocate both factors and inverse; do |

878 | * not do so yourself. |

879 | * |

880 | * @NOTE All the factorization methods use some sort of dense vectors |

881 | * for storing pivots and whatnot. As with a FORTRAN factorization, |

882 | * there's not much you can do with a factorization without |

883 | * also understanding in detail the factorization routine and getting |

884 | * your hands on the dense vectors. |

885 | * |

886 | * @todo Separate documentation of linsolqr_get_factors() |

887 | * and linsolqr_get_inverse()? |

888 | */ |

889 | |

890 | extern mtx_region_t *linsolqr_get_region(linsolqr_system_t sys); |

891 | /**< |

892 | * Returns a pointer to the current linsolqr region. |

893 | * This is being created for use in the check numeric dependency |

894 | * routine and may be removed once the new "problem manager" window |

895 | * is created to take over this functionality (and others). |

896 | */ |

897 | |

898 | extern int linsolqr_setup_ngslv(linsolqr_system_t sys, |

899 | real64 *rhs, |

900 | mtx_range_t *un_p_rng, |

901 | real64 *tmpvec); |

902 | /**< Sets up the NGSlv solver. */ |

903 | extern real64 *linsolqr_get_varvalue(linsolqr_system_t sys, int n); |

904 | /**< Returns the value of the nth variable in sys. */ |

905 | |

906 | /** @} */ |

907 | |

908 | #endif /* ASC_LINSOLQR_H */ |

john.pye@anu.edu.au | ViewVC Help |

Powered by ViewVC 1.1.22 |