Parent Directory | Revision Log

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

*Sat Apr 21 15:11:45 2007 UTC*
(17 years, 2 months ago)
by *jpye*

File MIME type: text/x-chdr

File size: 39874 byte(s)

File MIME type: text/x-chdr

File size: 39874 byte(s)

Added full jacobian preconditioner (ongoing). Some tweaks to dsgsat & related models.

1 | /* ASCEND modelling environment |

2 | Copyright (C) 1990 Karl Michael Westerberg |

3 | Copyright (C) 1993 Joseph Zaher |

4 | Copyright (C) 1994 Boyd Safrit, Joseph Zaher, Benjamin Andrew Allan |

5 | Copyright (C) 1995 Benjamin Andrew Allan, Kirk Andre Abbott |

6 | Copyright (C) 2006-2007 Carnegie Mellon University |

7 | |

8 | This program is free software; you can redistribute it and/or modify |

9 | it under the terms of the GNU General Public License as published by |

10 | the Free Software Foundation; either version 2, or (at your option) |

11 | any later version. |

12 | |

13 | This program is distributed in the hope that it will be useful, |

14 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

15 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |

16 | GNU General Public License for more details. |

17 | |

18 | You should have received a copy of the GNU General Public License |

19 | along with this program; if not, write to the Free Software |

20 | Foundation, Inc., 59 Temple Place - Suite 330, |

21 | Boston, MA 02111-1307, USA. |

22 | *//** @file |

23 | * linsolqr: Ascend Linear Equation Solver. |

24 | * |

25 | * A linear system consists of a coefficient matrix (A) |

26 | * and possibly several right-hand-sides (rhs). The |

27 | * solution vector x sought can be that of either |

28 | *<pre> |

29 | * T |

30 | * A x = rhs or A x = rhs |

31 | *</pre> |

32 | * depending on a specification inherent with rhs which |

33 | * dictates whether or not A should be transposed. If a |

34 | * rhs specifies transpose, then the vector itself is |

35 | * expected to be indexed by the original column numbers |

36 | * of the coefficient matrix and the solution vector shall |

37 | * be indexed by original row numbers. Otherwise, rhs |

38 | * is expected to be indexed by original row numbers while |

39 | * the solution can be referenced using original column |

40 | * numbers. The coefficient matrix and each rhs will be |

41 | * preserved throughout solving, except that the |

42 | * coefficient matrix may be permuted during reordering |

43 | * or (depending on the method) solving. |

44 | * |

45 | * The method of linear solution is an option new to |

46 | * linsolqr. Also, some of the data structures have been |

47 | * expanded as required. The factoring methods are |

48 | * enumerated and summarized below. |

49 | * |

50 | * Linsol is a general sparse solver. It uses mtx, though any |

51 | * proper sparse matrix package can be substituted. By |

52 | * proper, of course, we mean one with at least the |

53 | * functionality of mtx. If you find something with all |

54 | * the features of mtx that is faster than mtx, let us |

55 | * know, please. If you do any comparisons of linsol with |

56 | * other sparse codes (on sequential processors) also |

57 | * please let us know how you find its performance. |

58 | * You probably won't find algorithms for banded or |

59 | * symmetric matrices here. Why? |

60 | * The efficiency of most banded algorithms comes from |

61 | * the way they access memory with three assumptions: |

62 | * that all possible coefficient locations are known before |

63 | * a solution is attempted, that they are accessible |

64 | * in some ORDERED fashion, and that the matrix will not |

65 | * be reordered (except perhaps in very limited ways) |

66 | * during solution. Linsol is predicated on just |

67 | * the opposite assumptions at the moment. |

68 | * These assumptions would relieve the solver of |

69 | * memory allocation and a great deal of range checking. |

70 | * While we would do our flops efficiently with linsol/mtx |

71 | * on a banded matrix, the cost of accessing linear |

72 | * coefficients in an UNORDERED fashion is such that a good |

73 | * band, variable band, or symmetric skyline code making the |

74 | * ORDERED assumption will always beat linsol in speed, |

75 | * though perhaps not in robustness. |

76 | * |

77 | * Requires: #include <string.h> |

78 | * #include "utilities/ascConfig.h" |

79 | * #include "mtx.h" |

80 | * |

81 | * Each method is described here, including the unimplemented ones, |

82 | * however vaguely. By convention all methods assign unpivotable variables |

83 | * the value of 0.0 in the linear solution vector. This makes sense because |

84 | * the linear system usually arises from a nonlinear scheme wherein |

85 | * stepping in the unpivoted variables is left for a future iteration. |

86 | * |

87 | * The currently recommended methods are ranki_*2, in particular ba2. |

88 | * |

89 | * ranki_kw: |

90 | * Reordering method: SPK1, Stadtherr & Wood (1984), and like. |

91 | * The region used is confined to the main diagonal. |

92 | * The square region is reordered to lower triangular with spikes. |

93 | * Each spike is higher than or equal to all spikes to the left of it |

94 | * if the matrix is partitioned. |

95 | * Solution method: RANKI, Stadtherr & Wood |

96 | * Ref: Computers and Chemical Engineering (8)1 pp. 9-33 (1984) |

97 | * The square region is U/L factored with spikes dragged in to replace |

98 | * unacceptable pivots. Fill is in spike columns and |

99 | * dragged rows where it usually does take place. |

100 | * Numerically dependent rows and columns do not have dependencies |

101 | * automatically calculated. This is a change from linsol, so be |

102 | * careful in making comparisons. |

103 | * Authors: Karl Westerberg, Joe Zaher |

104 | * with some clean ups and restructuring by Ben Allan |

105 | * |

106 | * ranki_jz: |

107 | * As ranki_kw, except where ranki_kw looks only across the current row |

108 | * for pivots if the diagonal element is too small, ranki_jz also looks |

109 | * down the current column. This introduces additional spikes to the |

110 | * current block, but sometimes the overall effect is reduced fill. The |

111 | * main motivation is to broaden pivot choices on those problems |

112 | * where restricting the choice to the current row results in immediate |

113 | * explosions of fill due to a full subdiagonal and a moderately small |

114 | * pivot. |

115 | * |

116 | * ranki_kw2: |

117 | * ranki_jz2: |

118 | * As the ranki_XX factorizations, except uses drop tolerance and |

119 | * uses much faster bilevel matrices (master + slave matrix of mtx) |

120 | * to segregate the contents of U from L. Also, dependency information |

121 | * is not calculated automatically. |

122 | * Author: Ben Allan |

123 | * |

124 | * ranki_ba2: |

125 | * As the ranki_kw2 factorization, except pivot and row being eliminated |

126 | * are stored as linked lists/full org vectors (dual access) during |

127 | * elimination which removes the last significant quadratic effect from |

128 | * the ranki factorization. |

129 | * |

130 | * plain_qr: |

131 | * Reordering method: Transpose SPK1 (tspk1), LORA |

132 | * The region used may be rectangular, but will be modified |

133 | * if needed so that the upper left corner is on the diagonal. |

134 | * The region in sys->coef is not permuted in the solution process. |

135 | * Solution method: |

136 | * The rectangular region is QR factored with standard column pivoting |

137 | * modified by the pivot_tolerance. That is, at each step the column |

138 | * to factor in next is the leftmost col with norm alpha that satisfies |

139 | * alpha >= ptol * maxalpha. |

140 | * ptol = 1 --> Stewart's QRDC style pivoting. |

141 | * ptol = 0 --> no pivoting (Sometimes a dumb idea). |

142 | * plain_qr is to QRDC as ranki_kw is to LU factorization with strict |

143 | * partial pivoting. |

144 | * This implementation is best used on matrices where nrows about = ncols |

145 | * and data are not highly correlated. |

146 | * There are better QR methods (not available here) for vector processors, |

147 | * distributed processors, data fitting and so forth. |

148 | * Some care has been taken that the implementation here will scale up |

149 | * reasonably to larger matrices. |

150 | * Detailed references are given in the plain_qr section of the .c file. |

151 | * Author: Ben Allan |

152 | * |

153 | * cond_qr: (broken) |

154 | * Reordering method: Transpose SPK1 (tspk1) |

155 | * The region used is confined to the main diagonal. |

156 | * The square region is permuted in the solution process. |

157 | * Solution method: |

158 | * The square region is QR factored with pivoting by a local minimum |

159 | * fill criteria balanced against choosing pivots based on the |

160 | * incremental inverse R condition number. |

161 | * Ref:Incremental Condition Calculation and Column Selection, |

162 | * UMIACS TR 90-87, G.W.Stewart, July 1990) |

163 | * (Allan, Safrit, Westerberg, 1995) |

164 | * Authors: Ben Allan |

165 | * |

166 | * opt_qr: |

167 | * As cond_qr, except region may be underspecified (MxN, M < N) |

168 | * The region's upper left corner should be (or will be forced to be) |

169 | * on the main diagonal. The columns pivoted are expected to be a |

170 | * good basis. Solution is not a least squares solution of the |

171 | * transpose NxM problem; for this we need rowwise Householder |

172 | * transforms and spk1 reordering. |

173 | * Authors: Ben Allan |

174 | * |

175 | * ls_qr: |

176 | * As cond_qr, except region may be overspecified (MxN, M > N) |

177 | * Solves the system in a linear least squares sense. |

178 | * Authors: Ben Allan |

179 | * |

180 | * band_lu: |

181 | * symmetric_lu: |

182 | * (cholesky?) |

183 | * The best codes we can find already implemented to complement and/or |

184 | * verify our own. More detailed semantics unknown at this time. |

185 | * |

186 | * The parameters pivot_tolerance, condition_tolerance and drop_tolerance |

187 | * have semantics varying with the method. |

188 | * See the header of linsolqr_set_pivot_tolerance below for details. |

189 | * </pre> |

190 | *//* |

191 | by Karl Westerberg |

192 | Created: 2/6/90 |

193 | Last in CVS: $Revision: 1.13 $ $Date: 1997/07/25 17:23:49 $ $Author: ballan $ |

194 | QR options by Ben Allan |

195 | C version based on sqr.pas v1.5, Boyd Safrit (1994) |

196 | |

197 | Dates: 06/90 - original version (KW) |

198 | 04/91 - removed output assignment and partitioning |

199 | (belong in structural analysis) (JZ) |

200 | 08/92 - added transpose ability (JZ) |

201 | 01/94 - broke out linsol_invert() and linsol_solve() (JZ) |

202 | 10/94 - reorganized for multiple methods and QR (BA) |

203 | moved drag operator to mtx (BA,JZ) |

204 | 11/94 - added linsolqr_get_factors/inverse for debugging.BA |

205 | 09/95 - added factor_class to better manage families. BA |

206 | 09/95 - added ranki_*2 methods (BA) |

207 | 09/95 - added gauss (broken) (KA) |

208 | 11/96 - added ranki_ba2 method (BA) |

209 | |

210 | */ |

211 | |

212 | #ifndef ASC_LINSOLQR_H |

213 | #define ASC_LINSOLQR_H |

214 | |

215 | #include "mtx.h" |

216 | |

217 | /** @addtogroup linear Linear |

218 | @{ |

219 | */ |

220 | |

221 | #define LINSOLMTX_DEBUG FALSE |

222 | /**< Debug mode. */ |

223 | #define LINQR_DROP_TOLERANCE 1.0e-16 |

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

225 | |

226 | struct rhs_list; |

227 | |

228 | typedef struct linsolqr_header *linsolqr_system_t; |

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

230 | |

231 | enum factor_class { |

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

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

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

235 | }; |

236 | |

237 | enum reorder_method { |

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

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

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

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

242 | /* future work: |

243 | invspk1, spk1 then diagonally inverted |

244 | invtspk1, tspk1 then diagonally inverted |

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

246 | */ |

247 | }; |

248 | |

249 | enum factor_method { |

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

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

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

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

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

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

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

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

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

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

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

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

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

263 | }; |

264 | |

265 | ASC_DLLSPEC int g_linsolqr_timing; |

266 | /**< |

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

268 | * generates a fill and timing message. |

269 | */ |

270 | |

271 | /* KAA_DEBUG */ |

272 | /* unimplemented function |

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

274 | */ |

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

276 | |

277 | |

278 | /* Functions for managing interfaces */ |

279 | |

280 | extern char *linsolqr_rmethods(void); |

281 | /**< |

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

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

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

285 | * |

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

287 | */ |

288 | |

289 | ASC_DLLSPEC char *linsolqr_fmethods(void); |

290 | /**< |

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

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

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

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

295 | */ |

296 | |

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

298 | /**< |

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

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

301 | */ |

302 | |

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

304 | /**< |

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

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

307 | */ |

308 | |

309 | ASC_DLLSPEC enum factor_class linsolqr_fmethod_to_fclass(enum factor_method fm); |

310 | /**< |

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

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

313 | */ |

314 | |

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

316 | /**< |

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

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

319 | * Do not free the name. |

320 | */ |

321 | |

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

323 | /**< |

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

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

326 | * Do not free the name. |

327 | */ |

328 | |

329 | ASC_DLLSPEC const char *linsolqr_rmethod_description(const enum reorder_method meth); |

330 | /**< |

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

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

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

334 | */ |

335 | |

336 | ASC_DLLSPEC const char *linsolqr_fmethod_description(const enum factor_method meth); |

337 | /**< |

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

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

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

341 | */ |

342 | |

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

344 | |

345 | ASC_DLLSPEC linsolqr_system_t linsolqr_create(void); |

346 | /**< |

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

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

349 | */ |

350 | |

351 | ASC_DLLSPEC linsolqr_system_t linsolqr_create_default(void); |

352 | /**< |

353 | Creates a linsolqr system as per linsolqr_default but sets default |

354 | 'recommended' preferences for reordering and factoring methods, so the |

355 | ignorant user (me) doesn't have to worry about it every time they use this. |

356 | */ |

357 | |

358 | ASC_DLLSPEC void linsolqr_destroy(linsolqr_system_t sys); |

359 | /**< |

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

361 | * not destroyed by this call. |

362 | */ |

363 | |

364 | ASC_DLLSPEC void linsolqr_set_matrix(linsolqr_system_t sys, mtx_matrix_t mtx); |

365 | /**< |

366 | * Sets the coefficient matrix to mtx. |

367 | */ |

368 | |

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

370 | /**< |

371 | * Sets the reg to region. |

372 | */ |

373 | |

374 | ASC_DLLSPEC mtx_matrix_t linsolqr_get_matrix(linsolqr_system_t sys); |

375 | /**< |

376 | * Returns the coefficient matrix. |

377 | */ |

378 | |

379 | ASC_DLLSPEC void linsolqr_add_rhs(linsolqr_system_t sys, |

380 | real64 *rhs, |

381 | boolean transpose); |

382 | /**< |

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

384 | * system. |

385 | * |

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

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

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

389 | * |

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

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

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

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

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

395 | * <pre> |

396 | * For the mathematically impaired: |

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

398 | * |

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

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

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

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

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

404 | * |

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

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

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

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

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

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

411 | * mtx_transpose in order to improve factorization. |

412 | * |

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

414 | * </pre> |

415 | */ |

416 | |

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

418 | /**< |

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

420 | * removed from the system. |

421 | */ |

422 | |

423 | extern int32 linsolqr_number_of_rhs(linsolqr_system_t sys); |

424 | /**< |

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

426 | */ |

427 | |

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

429 | /**< |

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

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

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

433 | */ |

434 | |

435 | ASC_DLLSPEC void linsolqr_matrix_was_changed(linsolqr_system_t sys); |

436 | /**< |

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

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

439 | * made. |

440 | */ |

441 | |

442 | ASC_DLLSPEC void linsolqr_rhs_was_changed(linsolqr_system_t sys, |

443 | real64 *rhs); |

444 | /**< |

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

446 | * called whenever the rhs is modified. |

447 | */ |

448 | |

449 | extern void linsolqr_set_pivot_zero(linsolqr_system_t sys, |

450 | real64 pivot_zero); |

451 | /**< |

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

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

454 | * automatically call linsolqr_matrix_was_changed(). |

455 | */ |

456 | extern real64 linsolqr_pivot_zero(linsolqr_system_t sys); |

457 | /**< |

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

459 | * this value are regarded as zero. |

460 | */ |

461 | |

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

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

464 | extern real64 linsolqr_pivot_tolerance(linsolqr_system_t sys); |

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

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

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

468 | extern real64 linsolqr_condition_tolerance(linsolqr_system_t sys); |

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

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

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

472 | extern real64 linsolqr_drop_tolerance(linsolqr_system_t sys); |

473 | /**< |

474 | * <pre> |

475 | * linsolqr_set_pivot_tolerance(sys,ptol) |

476 | * ptol = linsolqr_pivot_tolerance(sys) |

477 | * |

478 | * linsolqr_set_condition_tolerance(sys,ctol) |

479 | * ctol = linsolqr_condition_tolerance(sys) |

480 | * |

481 | * linsolqr_set_drop_tolerance(sys,dtol) |

482 | * dtol = linsolqr_drop_tolerance(sys) |

483 | * |

484 | * linsolqr_system_t sys; |

485 | * real64 ptol, ctol,dtol; |

486 | * |

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

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

489 | * |

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

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

492 | * |

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

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

495 | * ranki_ba2: |

496 | * ranki_kw2: |

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

498 | * dtol in magnitude are dropped. |

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

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

501 | * |

502 | * cond_qr and variants: |

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

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

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

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

507 | * selected column. |

508 | * |

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

510 | * linsolqr_matrix_was_changed(). |

511 | * </pre> |

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

513 | */ |

514 | |

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

516 | |

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

518 | /**< |

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

520 | * The system should be previously prepped. |

521 | */ |

522 | |

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

524 | /**< |

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

526 | * The system should be previously prepped. |

527 | */ |

528 | |

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

530 | /**< |

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

532 | * The system should be previously prepped. |

533 | */ |

534 | |

535 | ASC_DLLSPEC int32 linsolqr_rank(linsolqr_system_t sys); |

536 | /**< |

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

538 | * factored. |

539 | */ |

540 | |

541 | extern real64 linsolqr_smallest_pivot(linsolqr_system_t sys); |

542 | /**< |

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

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

545 | * MAXDOUBLE is returned. |

546 | */ |

547 | |

548 | ASC_DLLSPEC int linsolqr_prep(linsolqr_system_t sys, enum factor_class fclass); |

549 | /**< |

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

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

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

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

554 | * |

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

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

557 | * the matrix before you call prep. |

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

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

560 | * |

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

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

563 | * prep again with the new class. |

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

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

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

567 | * |

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

569 | */ |

570 | |

571 | ASC_DLLSPEC int linsolqr_reorder(linsolqr_system_t sys, |

572 | mtx_region_t *region, |

573 | enum reorder_method rmeth); |

574 | /**< |

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

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

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

578 | * diagonal. |

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

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

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

582 | * |

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

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

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

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

587 | * not a particularly cheap search. |

588 | * |

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

590 | * takes care of any structural analysis necessary |

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

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

593 | * on structurally identical matrices. |

594 | * |

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

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

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

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

599 | * these methods. |

600 | * |

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

602 | * <pre> |

603 | * Brief notes on the reorderings available. |

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

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

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

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

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

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

610 | * </pre> |

611 | */ |

612 | |

613 | ASC_DLLSPEC int linsolqr_factor(linsolqr_system_t sys, |

614 | enum factor_method fmethod); |

615 | /**< |

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

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

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

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

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

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

622 | * with the class set when the prep was done, otherwise the call fails. |

623 | * |

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

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

626 | * status = linsolqr_factor(sys,linsolqr_fmethod(sys)); |

627 | * |

628 | * linsolqr_reorder must be called before factoring the matrix -- JP -- check? |

629 | * |

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

631 | */ |

632 | |

633 | extern int linsolqr_get_pivot_sets(linsolqr_system_t sys, |

634 | unsigned *org_rowpivots, |

635 | unsigned *org_colpivots); |

636 | /**< |

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

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

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

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

641 | * The system must be previously factored. |

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

643 | * There is no association of rows with columns here. |

644 | * |

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

646 | * changed. |

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

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

649 | */ |

650 | |

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

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

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

654 | /**< |

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

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

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

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

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

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

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

662 | * |

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

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

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

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

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

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

669 | */ |

670 | |

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

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

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

674 | /**< |

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

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

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

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

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

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

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

682 | * |

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

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

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

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

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

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

689 | */ |

690 | |

691 | extern int32 linsolqr_org_row_to_org_col(linsolqr_system_t sys, |

692 | int32 org_row); |

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

694 | extern int32 linsolqr_org_col_to_org_row(linsolqr_system_t sys, |

695 | int32 org_col); |

696 | /**< |

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

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

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

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

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

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

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

704 | */ |

705 | |

706 | ASC_DLLSPEC void linsolqr_calc_row_dependencies(linsolqr_system_t sys); |

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

708 | |

709 | ASC_DLLSPEC void linsolqr_calc_col_dependencies(linsolqr_system_t sys); |

710 | /**< |

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

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

713 | * called before either linsolqr_org_row/col_dependency |

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

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

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

717 | * calculated. |

718 | * |

719 | * The calculation method is fairly standard. |

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

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

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

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

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

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

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

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

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

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

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

731 | * work of a regular solve. |

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

733 | */ |

734 | |

735 | ASC_DLLSPEC mtx_sparse_t *linsolqr_row_dependence_coefs( |

736 | linsolqr_system_t sys, int32 orgrow); |

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

738 | |

739 | ASC_DLLSPEC mtx_sparse_t *linsolqr_col_dependence_coefs( |

740 | linsolqr_system_t sys, int32 orgcol |

741 | ); |

742 | /**< |

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

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

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

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

747 | * return is NULL. |

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

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

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

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

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

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

754 | * tation to the the caller. |

755 | * |

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

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

758 | * |

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

760 | */ |

761 | |

762 | extern real64 linsolqr_org_row_dependency(linsolqr_system_t sys, |

763 | int32 dep, int32 ind); |

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

765 | extern real64 linsolqr_org_col_dependency(linsolqr_system_t sys, |

766 | int32 dep, int32 ind); |

767 | /**< |

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

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

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

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

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

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

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

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

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

777 | * |

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

779 | */ |

780 | |

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

782 | /**< |

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

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

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

786 | * or not (see linsolqr_add_rhs.) |

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

788 | * in linsolqr_factor. |

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

790 | */ |

791 | |

792 | extern real64 linsolqr_var_value(linsolqr_system_t sys, |

793 | real64 *rhs, int32 var); |

794 | /**< |

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

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

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

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

799 | * original column number. |

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

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

802 | */ |

803 | |

804 | ASC_DLLSPEC boolean linsolqr_copy_solution(linsolqr_system_t sys, |

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

806 | /**< |

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

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

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

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

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

812 | * if the rhs is a transpose. |

813 | * |

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

815 | * |

816 | */ |

817 | |

818 | extern real64 linsolqr_eqn_residual(linsolqr_system_t sys, |

819 | real64 *rhs, int32 eqn); |

820 | /**< |

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

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

823 | * <pre> |

824 | * T |

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

826 | * </pre> |

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

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

829 | * be an original row number. |

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

831 | * value is MAXDOUBLE. |

832 | */ |

833 | |

834 | extern boolean linsolqr_calc_residual(linsolqr_system_t sys, |

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

836 | /**< |

837 | * Returns the residuals using the solution associated |

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

839 | * <pre> |

840 | * T |

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

842 | * </pre> |

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

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

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

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

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

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

849 | * |

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

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

852 | */ |

853 | |

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

855 | |

856 | extern size_t linsolqr_size(linsolqr_system_t sys); |

857 | /**< |

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

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

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

861 | * must do accounting for those. |

862 | */ |

863 | |

864 | extern void linsolqr_free_reused_mem(void); |

865 | /**< |

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

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

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

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

870 | */ |

871 | |

872 | /*----------------------------------------------------------------------------- |

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

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

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

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

877 | * functionality to linsol. |

878 | */ |

879 | |

880 | ASC_DLLSPEC mtx_matrix_t linsolqr_get_factors(linsolqr_system_t sys); |

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

882 | extern mtx_matrix_t linsolqr_get_inverse(linsolqr_system_t sys); |

883 | /**< |

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

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

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

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

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

889 | * not do so yourself. |

890 | * |

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

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

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

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

895 | * your hands on the dense vectors. |

896 | * |

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

898 | * and linsolqr_get_inverse()? |

899 | */ |

900 | |

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

902 | /**< |

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

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

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

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

907 | */ |

908 | |

909 | extern int linsolqr_setup_ngslv(linsolqr_system_t sys, |

910 | real64 *rhs, |

911 | mtx_range_t *un_p_rng, |

912 | real64 *tmpvec); |

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

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

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

916 | |

917 | /** @} */ |

918 | |

919 | #endif /* ASC_LINSOLQR_H */ |

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

Powered by ViewVC 1.1.22 |