Parent Directory | Revision Log

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

*Wed Jan 3 03:13:20 2007 UTC*
(15 years, 9 months ago)
by *johnpye*

File MIME type: text/x-chdr

File size: 23728 byte(s)

File MIME type: text/x-chdr

File size: 23728 byte(s)

Added idalinear.[ch], which will expose the linsolqr solver as a linear solver for IDA. Modified ida.c to support this new linear solver under the name 'ASCEND'.

1 | /* ASCEND modelling environment |

2 | Copyright (C) 1990 Karl Michael Westerberg |

3 | Copyright (C) 1993 Joseph Zaher |

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

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

6 | Copyright (C) 2006 Carnegie Mellon University |

7 | |

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

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

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

11 | any later version. |

12 | |

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

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

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

16 | GNU General Public License for more details. |

17 | |

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

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

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

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

22 | *//** @file |

23 | * mtx: Ascend Sparse Matrix Package. |

24 | * |

25 | * Description: This module allows the user to create and manipulate |

26 | * matrices. The following is list of what constitutes |

27 | * a "matrix": |

28 | * |

29 | * - A square nxn (where n is the order of the matrix) |

30 | * rectangular region of numbers, indexed by |

31 | * (row,col) pairs where 0<=row<n and 0<=col<n. |

32 | * |

33 | * - Row and column permutations which keep track of |

34 | * where a given row/column "came from" originally. |

35 | * |

36 | * The coefficient matrix (i.e. the nxn region of numbers) |

37 | * can be divided into two classes, non-zeros and zeros. |

38 | * Roughly speaking, a given element is a non-zero so long |

39 | * as its value has the POTENTIAL of not equaling zero. |

40 | * Thus, if the value of an element (r,c) is not equal to |

41 | * zero, then (r,c) must be classified as a non-zero: |

42 | * however the converse need not hold, unless either |

43 | * mtx_del_zr_in_row or mtx_del_zr_in_col has been |

44 | * recently called with r or c respectively. |

45 | * |

46 | * The mtx_matrix_t includes a data structure for block |

47 | * information storage, a block being a matrix subregion. |

48 | * This feature supports partitioning solvers. |

49 | * |

50 | * There are only 3 fundamental operations _on_ matrices |

51 | * vector scaling (mtx_mult_*) |

52 | * vector addition (mtx_add_*) |

53 | * permutation (mtx_swap_*) |

54 | * Matrix elements are maintained as |

55 | * relatively unordered linked lists. Mtx_next_in_* is |

56 | * most generally useful operator for doing sparse math |

57 | * with this matrix package. There are several operations |

58 | * that take matrices/vectors and return vectors/scalars. |

59 | * Sparse matrix-matrix computations are best coded with |

60 | * vector primitives and knowledge of the specific matrix, |

61 | * so they are not provided here. |

62 | * |

63 | * It's amazing for a package that only addresses the |

64 | * fundamentals, this still has over 100 calls and looks |

65 | * like an extremely well documented kitchen sink. |

66 | * |

67 | * The user may grep on extern to get a semiformatted list |

68 | * of the operator names in this file. Please follow the |

69 | * formatting if you add functions to this module. |

70 | * i.e. grep extern mtx2*h cooks up a reduced header. |

71 | * |

72 | ** vocabulary for this header file. |

73 | ** int32 a 32 bit signed integer (int usually, long on some archs.) |

74 | ** real64 a 64 bit real (double on most systems.) |

75 | ** By using these typedefs, all the code remains intact |

76 | ** across different systems/compilers: only the typedefs |

77 | ** have to change if there is a difference in data size. |

78 | ** |

79 | ** mtx, matrix: a peculiar sparse matrix structure which should not |

80 | ** be manipulated except by means of functions in this |

81 | ** header. Attempts to go around by another path are |

82 | ** 99.44% likely to fail, and 95% likely to fail |

83 | ** catastrophically. |

84 | ** |

85 | ** cur_row, |

86 | ** cur_col, or |

87 | ** row,col: an int32 denoting the ith(jth) row(column) as the |

88 | ** matrix is currently permuted (range 0..order-1). |

89 | ** org_row, org_col: |

90 | ** an int32 denoting the ith(jth) row(col) as the matrix |

91 | ** is ordered in the unpermuted state (range 0..order-1). |

92 | ** |

93 | ** vec: an array of real64, generally of size |

94 | ** mtx->order or larger. Indexing scheme is indeterminate. |

95 | ** org_vec: a vec indexed by the org_row or org_col index. |

96 | ** cur_vec: a vec indexed by the (current) row or col index. |

97 | ** |

98 | ** s,txxx: source or target xxx, where xxx may be a vec, row, col |

99 | ** org_row, or org_col. |

100 | ** |

101 | ** order: the working size of a matrix. |

102 | ** order N ==> cols/rows 0->N-1 exist. |

103 | ** capacity: the length, essentially, of the mtx edges. |

104 | ** This number is >= order, normally ==. |

105 | ** |

106 | ** coord: a pair of current indices (row,col). |

107 | ** range: a pair of current indices (low,high). |

108 | ** operations which take a range will do nothing if |

109 | ** low > high. |

110 | ** region: a pair of ranges (rowrange,colrange) |

111 | ** Coords, ranges and regions are passed by pointer to mtx functions. |

112 | ** |

113 | ** sparse: a compressed vector, with indeterminate indexing scheme. |

114 | ** a sparse has a capacity (maximum amount of data) |

115 | ** a length (the current amount of data) |

116 | ** an array of real64 (the data) and |

117 | ** an array of int32 (where the data correspond to.) |

118 | ** The array of int32 may be cur or org and this |

119 | ** is determined by the context in which the sparse is used. |

120 | ** |

121 | ** Sparses are generally passed by pointer. |

122 | *//* |

123 | 'mtx' by Karl Michael Westerberg, created 5/3/90 |

124 | 'mtx2' by Benjamin Andrew Allan |

125 | Last in CVS: $Revision: 1.11 $ $Date: 1998/07/05 20:46:27 $ $Author: ballan $ |

126 | |

127 | Dates: KW 06/90 - original version |

128 | JZ 01/94 - added output assignment and partitioning |

129 | operations allowing them to be performed |

130 | individually. |

131 | JZ 08/94 - reduced the internal practice of repeatedly |

132 | finding the existence of an element by re- |

133 | writing mtx_next_in_row and mtx_next_in_col |

134 | to use a static element pointer which can |

135 | advance and be rewound. Dramatic cpu time |

136 | savings were obtained by having the values |

137 | of subsequent elements returned by these |

138 | * functions to eliminate the need for calling |

139 | * mtx_value on the obtained coordinate. A new |

140 | * function mtx_add_value was added to replace |

141 | * mtx_set_value in applications where it is |

142 | * known by the user that no nonzero currently |

143 | * exists at the coordinate whose value is being |

144 | * set. |

145 | * BA 08/94 - Protoized header file to reduce coding errors. |

146 | * BA 11/94 - Added mtx_write_region as it is done elsewhere |

147 | * and seems a common lowlevel operation. |

148 | * Made error output file specifiable in __debug. |

149 | * Made functions explicitly void if needed. |

150 | * For the following, see headers for explanation. |

151 | * Added mtx_copy_complete and mtx_copy_wo_incidence. |

152 | * Added mtx_drag from JZ. |

153 | * Added mtx_write_region_plot from JZ for xgraphs. |

154 | * Added mtx_free_reused_mem. |

155 | * Added mtx_add_row_sparse and mtx_add_col_sparse. |

156 | * Added mtx_nonzeros_in_region. |

157 | * Added mtx_numbers_in_row. |

158 | * Added mtx_numbers_in_col. |

159 | * Added mtx_numbers_in_region. |

160 | * Annotated next_in operator descriptions. |

161 | * Renamed mtx_add_value to mtx_fill_value in |

162 | * anticipation of a set of fill commands and to |

163 | * keep the mtx_add namespace implying math only. |

164 | * BA 12/94 - Added a set of dense vector (double array) |

165 | * handling routines: |

166 | * mtx_org_row_vec \ |

167 | * mtx_org_col_vec | stuff an array which is cur or |

168 | * mtx_cur_row_vec | org indexed from a row/col. |

169 | * mtx_cur_col_vec / |

170 | * mtx_set_org_row_vec \ UNIMPLEMENTED |

171 | * mtx_set_org_col_vec | set a mtx row/col from |

172 | * mtx_set_cur_row_vec | cur or org indexed array. |

173 | * mtx_set_cur_col_vec / |

174 | * mtx_fill_org_row_vec \ |

175 | * mtx_fill_org_col_vec | set empty mtx row/col from |

176 | * mtx_fill_cur_row_vec | cur or org indexed array. |

177 | * mtx_fill_cur_col_vec / |

178 | * mtx_zr_org_vec_using_row \ |

179 | * mtx_zr_org_vec_using_col | zero org/cur vec entries |

180 | * mtx_zr_cur_vec_using_row | matching row/col incid. |

181 | * mtx_zr_cur_vec_using_col / |

182 | * mtx_org_vec_add_row \ |

183 | * mtx_org_vec_add_col | vec += mtx row/col is a |

184 | * mtx_cur_vec_add_row | cur/org indexed array. |

185 | * mtx_cur_vec_add_col / |

186 | * mtx_add_row_org_vec \ UNIMPLEMENTED |

187 | * mtx_add_col_org_vec | mtx row/col += vec. vec a |

188 | * mtx_add_row_cur_vec | cur/org indexed array. |

189 | * mtx_add_col_cur_vec / |

190 | * Added mtx_add_row/col_series for use in |

191 | * common elimination schemes. |

192 | * Added dot operators for dot products |

193 | * mtx_row_dot_full_org_vec \ |

194 | * mtx_col_dot_full_org_vec | row/col dot vec indexed |

195 | * mtx_row_dot_full_cur_vec | by org or cur row or col |

196 | * mtx_col_dot_full_cur_vec / |

197 | * |

198 | * Declared a sparse vector type that can be mapped |

199 | * directly to compressed row or column storage |

200 | * called mtx_sparse_t. |

201 | * |

202 | * Added mtx_write_region_human. |

203 | * Added mtx_read_region which is inverse of |

204 | * mtx_write_region. |

205 | * Revised mtx_row_max, mtx_col_max to return |

206 | * signed value of element with max abs(value). |

207 | * Consolidated mtx_copy* functions into one with |

208 | * a set of macros for the various flavors. |

209 | * |

210 | * 6/95 ba Added mtx_org_permute for 'natural' orderings. |

211 | * Added mtx_reverse_diagonal which does that. |

212 | * |

213 | * 8/95 ba Added slave matrices and attendant functions. |

214 | * mtx_create_slave adds a slave. |

215 | * mtx_chattel_size counts slave memory. |

216 | * |

217 | * Added map back functions with a drop tolerance. |

218 | * mtx_dropfill_cur_row_vec \ cur indexed vector |

219 | * mtx_dropfill_cur_col_vec / filled in mtx if big. |

220 | * Added a coordinate list struct for passing lists |

221 | * of elements around. |

222 | * |

223 | * 10/95 ba Added sparse vector fetch routines and destroy. |

224 | * mtx_org_row_sparse \ |

225 | * mtx_org_col_sparse | stuff a sparse which is cur or |

226 | * mtx_cur_row_sparse | org indexed from a row/col. |

227 | * mtx_cur_col_sparse / |

228 | * mtx_destroy_sparse > deallocate a sparse. |

229 | * 11/95 ba Added spec, but not code, for mtx_clear_rowlist |

230 | * in anticipation of nonlinear solvers recomputing |

231 | * the gradient of only the nonlinear rows. |

232 | * Added also: |

233 | * mtx_assemble, which does just that in much |

234 | * the finite element sense. |

235 | * mtx_fill_org_value which takes org coord. |

236 | * mtx_exception_recover which resets internals |

237 | * after floating point err has left things in |

238 | * an illdefined state. Also in this process |

239 | * unified several pieces of code that manage |

240 | * internal buffers. |

241 | * |

242 | * 12/95 ba mtx_steal_org_row_sparse \ clear incidence and |

243 | * mtx_steal_org_col_sparse | stuff a sparse which is |

244 | * mtx_steal_cur_row_sparse | cur or org indexed from |

245 | * mtx_steal_cur_col_sparse / a row/col. |

246 | * mtx_steal_org_row_vec \ clear incidence and |

247 | * mtx_steal_org_col_vec | stuff an array which is |

248 | * mtx_steal_cur_row_vec | cur or org indexed from |

249 | * mtx_steal_cur_col_vec / a row/col. |

250 | * mtx_get_pivot_col \ LU partial pivot selection |

251 | * mtx_get_pivot_row / with a sparsity tolerance. |

252 | * 1/96 ba mtx_fill_org_row_sparse \ |

253 | * mtx_fill_org_col_sparse | set empty row/col from |

254 | * mtx_fill_cur_row_sparse | cur or org indexed sparse. |

255 | * mtx_fill_cur_col_sparse / |

256 | * mtx_write_region_human \ Revised write human |

257 | * mtx_write_region_human_rows| so we can have row |

258 | * mtx_write_region_human_cols| or col oriented print. |

259 | * mtx_write_region_human_f / Supports old usage. |

260 | * |

261 | * 5/96 ba Split mtx.[ch] into several files to make life |

262 | * much easier for maintaining. All the mtx_*.c |

263 | * files are intended to stand together as a |

264 | * sparse matrix package. The split is for header |

265 | * digestibility and for ease of maintenance. |

266 | * The decomposition is as follows: |

267 | * |

268 | * mtx.h. |

269 | * |

270 | * This file, which includes the rest. |

271 | * This describes general mtx concepts. |

272 | * There is no mtx.c. |

273 | * |

274 | * mtx_basic.c mtx_basic.h Most memory management. |

275 | * mtx_perms.c mtx_perms.h Permuation management. |

276 | * mtx_query.c mtx_query.h Most queries, readonly math. |

277 | * mtx_linal.c mtx_linal.h Unjustifiable pieces. |

278 | * mtx_internal_use_only.[ch] Nobody's business. |

279 | * |

280 | * 4/97 ba Added mtx_transpose, mtx_isa_transpose to perms. |

281 | * Added mtx_duplicate_region to make a slave copy. |

282 | * |

283 | * 7/98 ba Added mtx_write_region_matlab for harwell format. |

284 | * |

285 | * |

286 | ** header conventions: (DON'T SKIP READING THIS SECTION!!) *yawn* |

287 | ** -$- in place of *** implies as described in compilation flags below. |

288 | ** ! ! in place of *** implies a warning which should be read before |

289 | ** ! ! using the function under discussion in certain ways. |

290 | ** |

291 | ** compilation flags: |

292 | ** MTX_DEBUG - if TRUE, check_matrix will be invoked all over the place |

293 | ** and -$- returns will be in effect. |

294 | ** check_sparse is also used where appropriate if TRUE. |

295 | ** - if FALSE, routines marked with -$- will assume that you |

296 | ** are sending in a good mtx_matrix_t. Check accordingly. |

297 | ** They will also assume any mtx_sparse_t * is good. |

298 | ** MTX_DEBUG is wired in file mtx.c, but shouldn't be. |

299 | ** MTX_DEBUG TRUE slows this code down incredibly and should not be |

300 | ** used when compiling production code. |

301 | ** |

302 | * right. let's get into it... |

303 | */ |

304 | |

305 | #ifndef ASC_MTX_H |

306 | #define ASC_MTX_H |

307 | |

308 | #include <stdio.h> |

309 | #include <string.h> |

310 | #include <utilities/ascConfig.h> |

311 | |

312 | /***********************************************************************\ |

313 | public mtx data structures |

314 | \***********************************************************************/ |

315 | typedef struct mtx_header *mtx_matrix_t; |

316 | /**< Handle to the matrix */ |

317 | |

318 | /* Used to index rows and columns of matrices */ |

319 | |

320 | /** |

321 | *** Refers to (row,col) element of a matrix |

322 | **/ |

323 | typedef struct mtx_coord_struct { |

324 | int32 row,col; |

325 | } mtx_coord_t; |

326 | |

327 | /** |

328 | *** provide a list type for giving functions a list of coordinates. |

329 | **/ |

330 | struct mtx_coord_list { |

331 | mtx_coord_t coord; |

332 | struct mtx_coord_list *next; |

333 | }; |

334 | |

335 | typedef struct mtx_range_struct { |

336 | int32 low,high; |

337 | } mtx_range_t; |

338 | /**< |

339 | *** Range refered to is low..high, inclusive. If low>high, range is empty |

340 | **/ |

341 | |

342 | typedef struct mtx_region_struct { |

343 | mtx_range_t row,col; |

344 | } mtx_region_t; |

345 | /**< |

346 | *** Rectangular region of a matrix |

347 | **/ |

348 | |

349 | typedef struct mtx_block_struct { |

350 | int32 nblocks; |

351 | mtx_region_t *block; /**< array of regions denoting blocks */ |

352 | } mtx_block_t; |

353 | /**< |

354 | *** Block structure of a matrix. |

355 | **/ |

356 | |

357 | typedef struct mtx_sparse_vector_struct { |

358 | real64 *data; |

359 | int32 *idata; |

360 | int32 cap; |

361 | int32 len; |

362 | } mtx_sparse_t; |

363 | /**< |

364 | *** A structure for holding a compressed row or column vector. |

365 | *** The semantics imposed for this structure are very limited |

366 | *** so that one may be used (or reused) in many contexts. |

367 | *** An mtx_sparse_t doesn't know where its data came from, or by |

368 | *** what scheme the data is indexed.<br><br> |

369 | *** |

370 | *** The value of data[i] goes with mtx location somehow indexed idata[i]. |

371 | *** |

372 | *** The values of idata[i] do not have to be in any particular |

373 | *** order, and the vector knows nothing about whether the indices |

374 | *** stored in idata are row indices, column indices, org_row indices |

375 | *** or org_column indices. This information is determined by the context |

376 | *** of the mtx_sparse_t.<br><br> |

377 | *** |

378 | *** cap is the size of data and idata. |

379 | *** len is the length of data and idata (starting from location 0) |

380 | *** which contains valid data.<br><br> |

381 | *** |

382 | *** A value of -1 for len indicates that pointers data and idata are invalid, |

383 | *** as does a value of 0 for cap.<br><br> |

384 | *** |

385 | *** NOTE: The mtx_sparse_t is not a pointer to a struct, it IS the struct. |

386 | **/ |

387 | |

388 | typedef struct mtx_block_perm_structure *mtx_block_perm_t; |

389 | /**< |

390 | *** The mtx_block_perm_structure contains a copy of the information |

391 | *** needed to apply a previously derived and saved ordering to the |

392 | *** matrix in question. Precisely what this information is is nobody's |

393 | *** business. We reserve the right to change at a moment's notice what |

394 | *** the mtx_block_perm_structure actually is.<br><br> |

395 | *** |

396 | *** The utility of the mtx_block_perm_structure is as follows: |

397 | *** Say in a sequence of nonlinear solves, one wishes to solve |

398 | *** many matrices of the same nonzero pattern, but different |

399 | *** nonzero values. One wants to reorder a specified block |

400 | *** of the matrix one time and reuse that initial ordering at each |

401 | *** factorization. Now suppose one also wants to use a linear |

402 | *** solution method that messes with the ordering during factorization. |

403 | *** One can do the following: |

404 | *** Create the mtx in question, output assign and partition it.<br><br> |

405 | *** <pre> |

406 | *** Call bp=mtx_create_block_perm(mtx). |

407 | *** This returns a copy of the information we need to store. |

408 | *** Working with individual blocks is then as follows: |

409 | *** 1- Reorder the block in question in the mtx in any way you like. |

410 | *** 2- Call mtx_update_block_perm(mtx,block_number,bp) |

411 | *** This updates the saved information for rows and columns |

412 | *** *** within *** the block specified by copying from the mtx. |

413 | *** All sorts of inexpensive sanity checks will be performed to |

414 | *** ensure that this is a wise thing to do. |

415 | *** 3- Factor, or otherwise do anything you want to the permutation |

416 | *** of the original mtx, so long as you do not mess with more than |

417 | *** one block. Changing the order (size) of the mtx is not permitted. |

418 | *** 4- Call mtx_restore_block_perm(mtx,block_number,bp) |

419 | *** This takes information stored in bp and resets the permutation |

420 | *** of the mtx within block_number to match it. |

421 | *** All sorts of inexpensive sanity checks will be performed to |

422 | *** ensure that this is a wise thing to do. |

423 | *** 5- Repeat steps 3 and 4 ad nauseum, or loop back through 1 and 2 |

424 | *** if moving from block to block. |

425 | *** 6- When done with the matrix, or at least the process in question, |

426 | *** don't forget to do an mtx_destroy_block_perm(bp). |

427 | *** </pre> |

428 | *** Other variations on this scheme are possible. See the headers for |

429 | *** details of the individual functions. |

430 | **/ |

431 | |

432 | #define mtx_ENTIRE_MATRIX NULL |

433 | /**< |

434 | *** Refers to the entire matrix in any function where a region is called for. |

435 | **/ |

436 | |

437 | #define mtx_ALL_ROWS NULL |

438 | #define mtx_ALL_COLS NULL |

439 | /**< |

440 | *** Refers to all rows/columns in any function where a row/column range |

441 | *** is called for. |

442 | **/ |

443 | |

444 | #define mtx_CREATE_SPARSE NULL |

445 | /**< |

446 | *** In calls that take an mtx_sparse_t as an argument AND return a |

447 | *** pointer to a sparse AND do not note otherwise in the header |

448 | *** mtx_CREATE_SPARSE can be passed the argument and the call will |

449 | *** allocate the sparse_t data structure it returns. In such cases, |

450 | *** it is the users responsibility to see that the sparse is |

451 | *** eventually destroyed. |

452 | **/ |

453 | |

454 | #define mtx_IGNORE_ZEROES (1) |

455 | #define mtx_SOFT_ZEROES (0) |

456 | /**< |

457 | *** Those functions which fill a sparse_t might optionally return |

458 | *** soft zeroes given mtx_SOFT_ZEROES. |

459 | **/ |

460 | |

461 | #define mtx_FIRST ((int32)(-1)) |

462 | #define mtx_LAST ((int32)(-1)) |

463 | #define mtx_NONE ((int32)(-1)) |

464 | #define mtx_ALL_BLOCKS ((int32)(-1)) |

465 | |

466 | #include "mtx_basic.h" |

467 | /** |

468 | *** All operations involving accounting/memory |

469 | *** management and CHANGES to matrix data structure. |

470 | *** This includes data file i/o and floating point |

471 | *** exception recovery handling. |

472 | **/ |

473 | |

474 | #include "mtx_perms.h" |

475 | #include "mtx_reorder.h" |

476 | /*** |

477 | *** Functions for permuting matrices, reordering, etc. |

478 | *** These probably ought to be combined. mtx_reorder is currently |

479 | *** thinly disguised linsol code that needs to be redone in the mtx idiom. |

480 | **/ |

481 | |

482 | #include "mtx_query.h" |

483 | /*** |

484 | *** READ_ONLY functions that don't make more sense |

485 | *** in mtx2_perms or mtx2_basic. Most vector data |

486 | *** and basic math operations fall in this file. |

487 | **/ |

488 | |

489 | #include "mtx_linal.h" |

490 | /*** |

491 | *** Functions which probably shouldn't exist at all, |

492 | *** but are simply too damned expensive to do without |

493 | *** knowing all the internals of the mtx. These tend |

494 | *** to be ultracritical path functions in linsol. |

495 | *** These functions generally have very, very |

496 | *** complex semantics and preconditions and assumptions. |

497 | *** You really DON'T want to add functions to this |

498 | *** file or use those that are in it. |

499 | **/ |

500 | |

501 | #endif /* ASC_MTX_H */ |

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

Powered by ViewVC 1.1.22 |