Parent Directory | Revision Log

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

*Fri Oct 29 20:54:12 2004 UTC*
(18 years, 5 months ago)
by *aw0a*

File MIME type: text/x-chdr

File size: 42413 byte(s)

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 |

Name | Value |
---|---|

svn:executable |
* |

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

Powered by ViewVC 1.1.22 |