1 |
aw0a |
1 |
#ifdef __MTX_C_SEEN__ |
2 |
jds |
54 |
/* |
3 |
aw0a |
1 |
* mtx2: Ascend Sparse Matrix Package |
4 |
|
|
* by Benjamin Andrew Allan |
5 |
|
|
* Derived from mtx by Karl Michael Westerberg |
6 |
|
|
* Created: 5/3/90 |
7 |
|
|
* Version: $Revision: 1.9 $ |
8 |
|
|
* Version control file: $RCSfile: mtx_use_only.h,v $ |
9 |
|
|
* Date last modified: $Date: 2000/01/25 02:27:13 $ |
10 |
|
|
* Last modified by: $Author: ballan $ |
11 |
|
|
* |
12 |
|
|
* This file is part of the SLV solver. |
13 |
|
|
* |
14 |
|
|
* Copyright (C) 1996 Benjamin Andrew Allan |
15 |
|
|
* based (loosely) on mtx |
16 |
|
|
* Copyright (C) 1990 Karl Michael Westerberg |
17 |
|
|
* Copyright (C) 1993 Joseph Zaher |
18 |
|
|
* Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan |
19 |
|
|
* Copyright (C) 1995 Kirk Andre Abbott, Benjamin Andrew Allan |
20 |
|
|
* |
21 |
|
|
* The SLV solver is free software; you can redistribute |
22 |
|
|
* it and/or modify it under the terms of the GNU General Public License as |
23 |
|
|
* published by the Free Software Foundation; either version 2 of the |
24 |
|
|
* License, or (at your option) any later version. |
25 |
|
|
* |
26 |
|
|
* The SLV solver is distributed in hope that it will be |
27 |
|
|
* useful, but WITHOUT ANY WARRANTY; without even the implied warranty of |
28 |
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
29 |
|
|
* General Public License for more details. |
30 |
|
|
* |
31 |
|
|
* You should have received a copy of the GNU General Public License along with |
32 |
|
|
* the program; if not, write to the Free Software Foundation, Inc., 675 |
33 |
|
|
* Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING. |
34 |
|
|
* COPYING is found in ../compiler. |
35 |
|
|
*/ |
36 |
|
|
|
37 |
jds |
54 |
/** @file |
38 |
|
|
* mtx2: Ascend Sparse Matrix Package (Private). |
39 |
|
|
* |
40 |
|
|
* This file defines the private parts of an mtx and is only for mtx*.c |
41 |
|
|
* consumption. Any temptation to include this header in a linear or |
42 |
|
|
* nonlinear solver package is a symptom of extremely bad programming |
43 |
|
|
* and lack of proper task analysis. This header should be regarded as |
44 |
|
|
* highly unstable. We make ABSOLUTELY NO commitment to maintain |
45 |
|
|
* consistency between any two versions of this file.<br><br> |
46 |
|
|
* |
47 |
|
|
* Note to third parties: |
48 |
|
|
* mtx is PRODUCTION code in very long use at Carnegie Mellon University. |
49 |
|
|
* * As such, we maintain a very tight hold of the internals of our data |
50 |
|
|
* structure so that we can easily prove the code when apparent bugs arise. |
51 |
|
|
* 99.44% of "bugs" experienced in using mtx are the result of not |
52 |
|
|
* reading the public headers carefully.<br><br> |
53 |
|
|
* |
54 |
|
|
* The material in this file was never a part of any header until the |
55 |
|
|
* old mtx.c file got so big that we had to split it up to make it |
56 |
|
|
* manageable.<br><br> |
57 |
|
|
* |
58 |
|
|
* Note to future developers of the mtx module. If you change ANYTHING |
59 |
|
|
* in this file it is YOUR job to: |
60 |
|
|
* a) clear that change with all the other developers using this header |
61 |
|
|
* b) fix ALL the other mtx*.c files that depend on it. |
62 |
|
|
* If you are not willing to do that much work, why the hell are you |
63 |
|
|
* dabbling in sparse matrix math? go work on GUIs. |
64 |
|
|
* <pre> |
65 |
|
|
* requires: #include <stdio.h> |
66 |
|
|
* requires: #include "utilities/ascConfig.h" |
67 |
|
|
* requires: #include "mem.h" |
68 |
|
|
* requires: #include "mtx.h" |
69 |
|
|
* <pre> |
70 |
|
|
*/ |
71 |
aw0a |
1 |
|
72 |
|
|
#ifndef __MTX_INTERNAL_USE_ONLY_H__ |
73 |
|
|
#define __MTX_INTERNAL_USE_ONLY_H__ |
74 |
|
|
|
75 |
|
|
#ifndef FALSE |
76 |
jds |
54 |
/** These should have come in from ascConfig.h. |
77 |
|
|
* @todo Remove redefines of FALSE & TRUE to enforce pre-inclusion of ascCenfig.h? */ |
78 |
aw0a |
1 |
#define FALSE 0 |
79 |
|
|
#define TRUE 1 |
80 |
|
|
#endif |
81 |
|
|
|
82 |
|
|
#define MTX_DEBUG FALSE |
83 |
ben.allan |
33 |
/**< MTX_DEBUG is a no holds barred sanity checking flag for use when |
84 |
aw0a |
1 |
* nothing else is giving a clue why something is going wrong. It |
85 |
|
|
* slows down the code to a crawl. Do not under any conditions change |
86 |
|
|
* its value or undefine it except at this location. If you need some |
87 |
|
|
* other sort of debugging flag for debugging a particular function, |
88 |
|
|
* use some personal debugging flag. |
89 |
|
|
*/ |
90 |
|
|
|
91 |
|
|
#define EVEN FALSE |
92 |
|
|
#define ODD TRUE |
93 |
|
|
#define SWAPS_PRESERVE_ORDER TRUE |
94 |
ben.allan |
33 |
/**< |
95 |
aw0a |
1 |
*** Do row and column swaps preserve the ordering of non-zeros in rows |
96 |
|
|
*** and columns? Setting this to TRUE means swapping only entails the |
97 |
|
|
*** movement of integer row or column numbers and NOT the exchange of |
98 |
|
|
*** entire row or columns. |
99 |
|
|
**/ |
100 |
|
|
#define WIDTHMAGIC 2048 |
101 |
ben.allan |
33 |
/**< |
102 |
aw0a |
1 |
*** WIDTHMAGIC is the byte size to aim for in allocating groups of elements. |
103 |
jds |
54 |
**/ |
104 |
|
|
#define LENMAGIC 10 |
105 |
|
|
/**< |
106 |
aw0a |
1 |
*** LENMAGIC initial # of groups of elements, hence the smallest |
107 |
|
|
*** possible number of elements a matrix will ever have is LENM*WIDTHM/eltsize. |
108 |
|
|
**/ |
109 |
|
|
|
110 |
|
|
extern FILE *g_mtxerr; |
111 |
jds |
54 |
/**< |
112 |
aw0a |
1 |
*** Global file pointer to which errors are reported. Should never be |
113 |
|
|
*** NULL. Also useful when running ascend in gdb and you can't find |
114 |
|
|
*** any other file pointer to use. |
115 |
|
|
**/ |
116 |
|
|
|
117 |
jds |
54 |
/** just a struct to make resulting code more readable. */ |
118 |
aw0a |
1 |
struct next_element_t { |
119 |
|
|
struct element_t *row; |
120 |
|
|
struct element_t *col; |
121 |
|
|
}; |
122 |
|
|
|
123 |
jds |
54 |
/** |
124 |
aw0a |
1 |
*** This is the basic jacobian element of an mtx. |
125 |
|
|
*** It's size is 24 bytes on 4 byte pointer machines and |
126 |
|
|
*** 32 bytes on 8 byte pointer machines. |
127 |
|
|
*** The elements form a bidirectional singly linked list. |
128 |
|
|
*** The row and col indices in an element refer back to |
129 |
|
|
*** the header positions of the two lists that element is in. |
130 |
|
|
*** That is, each element knows its orgrow and orgcol. |
131 |
|
|
**/ |
132 |
jds |
54 |
struct element_t { |
133 |
|
|
real64 value; |
134 |
|
|
int32 row; |
135 |
|
|
int32 col; |
136 |
|
|
struct next_element_t next; |
137 |
|
|
}; |
138 |
aw0a |
1 |
|
139 |
jds |
54 |
/** |
140 |
aw0a |
1 |
*** Each matrix is really just a pair of arrays of pointers to |
141 |
|
|
*** elements. The index of a row or column in THESE arrays is |
142 |
|
|
*** what is referred to as an org index. A value of NULL in |
143 |
jds |
54 |
*** either array means that that row (or col) is empty.<br><br> |
144 |
|
|
*** |
145 |
aw0a |
1 |
*** When we insert elements in the matrix, we simply shove the |
146 |
|
|
*** element in at the head of the its row/column lists. |
147 |
|
|
*** When we delete an element in the matrix, we search in one |
148 |
|
|
*** direction and unlink the element, marking it "dead". Then a |
149 |
|
|
*** general pass in the other direction unlinks all the "dead" |
150 |
jds |
54 |
*** elements.<br><br> |
151 |
|
|
*** |
152 |
aw0a |
1 |
*** Special note: The -1th element of nz_header arrays is NOT allocated. |
153 |
|
|
**/ |
154 |
jds |
54 |
struct nz_headers_t { |
155 |
|
|
struct element_t **row; |
156 |
|
|
struct element_t **col; |
157 |
|
|
}; |
158 |
aw0a |
1 |
|
159 |
jds |
54 |
/**< |
160 |
aw0a |
1 |
*** We maintain, rather than rederiving, the information required to |
161 |
|
|
*** answer all possible permutation questions. |
162 |
|
|
*** This is a policy decision based on the fact that mtx is research |
163 |
jds |
54 |
*** code that needs maximal flexibility at reasonable speed.<br><br> |
164 |
|
|
*** |
165 |
aw0a |
1 |
*** The -1th element of org_to_cur and cur_to_org are defined because |
166 |
|
|
*** -1 is used all over mtx as an error return. It's easier to debug |
167 |
|
|
*** things without the memory access errors that would happen if |
168 |
jds |
54 |
*** -1 were not allocated or were part of memory in some other object.<br><br> |
169 |
aw0a |
1 |
*** |
170 |
jds |
54 |
*** Special note: The -1th element of nz_header arrays is NOT allocated.<br><br> |
171 |
|
|
*** |
172 |
|
|
*** Do not access the parity field of a slave matrix, refer to its master. |
173 |
|
|
*** Conduct all permuting operations on the master. |
174 |
aw0a |
1 |
**/ |
175 |
jds |
54 |
struct permutation_t { |
176 |
|
|
int32 *org_to_cur; /**< org_to_cur[-1] = -1 */ |
177 |
|
|
int32 *cur_to_org; /**< cur_to_org[-1] = -1 */ |
178 |
|
|
boolean parity; |
179 |
|
|
}; |
180 |
aw0a |
1 |
|
181 |
|
|
struct permutations_t { |
182 |
|
|
struct permutation_t row; |
183 |
|
|
struct permutation_t col; |
184 |
|
|
int32 transpose; |
185 |
|
|
}; |
186 |
|
|
|
187 |
jds |
54 |
/**< |
188 |
aw0a |
1 |
*** There is a list of blocks associated with a matrix. |
189 |
|
|
*** This is an artifact of POOR solver API design between |
190 |
|
|
*** Peter Piela and Karl Westerberg. The blockwise decomposition |
191 |
|
|
*** information properly belongs to a linear or nonlinear solver |
192 |
|
|
*** and not to the mtx. |
193 |
|
|
*** |
194 |
jds |
54 |
*** @todo We intend to fix this soon. |
195 |
aw0a |
1 |
**/ |
196 |
jds |
54 |
struct structural_data_t { |
197 |
|
|
int32 symbolic_rank; /**< Symbolic rank (< 0 if invalid) */ |
198 |
|
|
int32 nblocks; /**< # blocks in matrix */ |
199 |
|
|
mtx_region_t *block; /**< Pointer to array of blocks */ |
200 |
|
|
}; |
201 |
aw0a |
1 |
|
202 |
jds |
54 |
/**< |
203 |
aw0a |
1 |
*** capacity may be > order. |
204 |
|
|
*** A matrix of capacity 0 doesn't have a mem_store_t yet and elements |
205 |
|
|
*** cannot be queried about without a core dump. |
206 |
|
|
**/ |
207 |
jds |
54 |
struct mtx_header { |
208 |
|
|
int integrity; /**< Integrity integer */ |
209 |
|
|
int32 order; /**< Order of the matrix */ |
210 |
|
|
int32 capacity; /**< Capacity of all the arrays */ |
211 |
|
|
int32 nslaves; /**< number of slave matrices */ |
212 |
|
|
struct nz_headers_t hdr; /**< Non-zero headers of the matrix */ |
213 |
|
|
struct element_t *last_value; /**< value/set_value memory */ |
214 |
|
|
mem_store_t ms; /**< element cache memory */ |
215 |
|
|
struct permutations_t perm; /**< Permutation vectors */ |
216 |
|
|
struct structural_data_t *data; /**< Pointer to structural information */ |
217 |
|
|
mtx_matrix_t master; /**< the master of this mtx, if slave */ |
218 |
|
|
mtx_matrix_t *slaves; /**< array of slave matrices */ |
219 |
|
|
}; |
220 |
aw0a |
1 |
|
221 |
jds |
54 |
/**< |
222 |
aw0a |
1 |
*** If you want to save a permutation for restoration, you |
223 |
|
|
*** have to make a copy of that data, eh? Here's the place you |
224 |
|
|
*** put it. Note that the block list should be disappearing from |
225 |
|
|
*** from the structural data soon. |
226 |
|
|
**/ |
227 |
jds |
54 |
struct mtx_block_perm_structure { |
228 |
|
|
int integrity; |
229 |
|
|
int32 order; /**< Order of the matrix */ |
230 |
|
|
int32 capacity; /**< Capacity of all the arrays */ |
231 |
|
|
mtx_matrix_t mtx; /**< matrix of origin */ |
232 |
|
|
struct permutations_t perm; /**< Permutation vectors */ |
233 |
|
|
struct structural_data_t *data; /**< Pointers to structural information */ |
234 |
|
|
}; |
235 |
aw0a |
1 |
|
236 |
|
|
#define OK ((int)201539237) |
237 |
jds |
54 |
/**< Matrix integrity (ok) value. */ |
238 |
aw0a |
1 |
#define DESTROYED ((int)531503871) |
239 |
jds |
54 |
/**< matrix integrity (destroyed) value. */ |
240 |
aw0a |
1 |
|
241 |
jds |
54 |
#define ZERO ((int32)0) |
242 |
|
|
#define D_ZERO ((real64)0.0) |
243 |
|
|
#define D_ONE ((real64)1.0) |
244 |
ben.allan |
33 |
/**< useful constants if your C compiler is not too bright about ANSI */ |
245 |
aw0a |
1 |
|
246 |
|
|
#define ISSLAVE(m) ((m)->master!=NULL) |
247 |
jds |
54 |
/**< Returns 1 if m is a slave matrix, 0 if not. */ |
248 |
aw0a |
1 |
|
249 |
|
|
#define ordered3(a,b,c) ((a) <= (b) && (b) <= (c)) |
250 |
|
|
#define in_range(rng,ndx) ordered3((rng)->low,ndx,(rng)->high) |
251 |
|
|
#define legal(mtx,ndx) ordered3(ZERO,ndx,(mtx)->order-1) |
252 |
ben.allan |
33 |
/**< |
253 |
aw0a |
1 |
*** Boolean operators to compare a row or column |
254 |
|
|
*** index with some specified range or the maximum |
255 |
|
|
*** range of the matrix in which it is used. |
256 |
|
|
**/ |
257 |
|
|
|
258 |
|
|
#define fast_in_range(l,h,i) ( ordered3(l,i,h) ) |
259 |
|
|
#define not_in_range(l,h,i) ( (i)<(l) || (i)>(h) ) |
260 |
ben.allan |
33 |
/**< |
261 |
aw0a |
1 |
*** Boolean operators to compare 3 integers. |
262 |
|
|
*** l <= h must be TRUE or these will lie. In many cases, |
263 |
|
|
*** this condition can (or should) be met before in_range |
264 |
|
|
*** is called. Sometimes these are not faster since the lo,hi vals cost. |
265 |
|
|
*** In particular, queries like next_col do not profit while calls |
266 |
|
|
*** which must traverse an entire row/col do. |
267 |
|
|
*** Gains in cycle count on dec alphas+cc are about 10% per function, |
268 |
|
|
*** but the gains in time are more like 1%, so alpha pixie is lying a little. |
269 |
|
|
*** For compilers which are not as clever as Decs, (gcc, sun acc) the |
270 |
|
|
*** gains should be much more visible. (some do not realize rng->low |
271 |
|
|
*** is invariant even with -O.) |
272 |
|
|
*** Note that these are 'loose' comparisons if !(l<=h) |
273 |
|
|
**/ |
274 |
|
|
|
275 |
|
|
#define zero(ptr,nelts,type) \ |
276 |
|
|
mem_zero_byte_cast((ptr),0,(nelts)*sizeof(type)) |
277 |
ben.allan |
33 |
/**< |
278 |
aw0a |
1 |
*** Zeros a vector of specified length and type. |
279 |
|
|
*** It is inefficient to use, however, if you know the type |
280 |
|
|
*** is one of the basic types (int,double,ptr,char) |
281 |
|
|
**/ |
282 |
|
|
|
283 |
|
|
|
284 |
jds |
54 |
/* ************************************************************************ *\ |
285 |
aw0a |
1 |
Private check routines |
286 |
jds |
54 |
\* ************************************************************************ */ |
287 |
|
|
extern int super_check_matrix(mtx_matrix_t mtx); |
288 |
ben.allan |
33 |
/**< |
289 |
aw0a |
1 |
*** After somevery extensive checking, returns an error count. |
290 |
|
|
*** More or less assume MTX_DEBUG is TRUE, and that is the only |
291 |
|
|
*** condition under which this should be called. |
292 |
|
|
**/ |
293 |
|
|
|
294 |
jds |
54 |
/* ************************************************************************ *\ |
295 |
aw0a |
1 |
Element CREATE/find routines. Please try to confine use of these to |
296 |
|
|
mtx_basic.c as much as possible. |
297 |
jds |
54 |
Use of find should be avoided at all costs, and in particular |
298 |
aw0a |
1 |
absolutely noone outside mtx should put their fingers on elements. |
299 |
|
|
|
300 |
|
|
These functions are not exported to generic users because they are |
301 |
|
|
on the critical path and we cannot afford the sanity checking required. |
302 |
|
|
They should only be called in contexts where the arguments are |
303 |
|
|
guaranteed valid. |
304 |
jds |
54 |
\* ************************************************************************ */ |
305 |
aw0a |
1 |
|
306 |
jds |
54 |
struct element_t *mtx_find_element(mtx_matrix_t mtx, |
307 |
|
|
int32 org_row, |
308 |
|
|
int32 org_col); |
309 |
ben.allan |
33 |
/**< |
310 |
jds |
54 |
*** <!-- mtx_find_element(mtx,org_row,org_col) --> |
311 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
312 |
|
|
*** <!-- int32 org_row; --> |
313 |
|
|
*** <!-- int32 org_col; --> |
314 |
aw0a |
1 |
*** |
315 |
|
|
*** Searches for a given element of the matrix and returns a pointer to it |
316 |
|
|
*** if it exists, or NULL if it doesn't exist. |
317 |
|
|
*** It is *ASSUMED* that org_row |
318 |
|
|
*** and org_col are legal indices. May crash if they are not. |
319 |
|
|
**/ |
320 |
|
|
|
321 |
jds |
54 |
struct element_t *mtx_create_element(mtx_matrix_t mtx, |
322 |
|
|
int32 org_row, |
323 |
|
|
int32 org_col); |
324 |
|
|
/**< |
325 |
|
|
*** <!-- mtx_create_element(mtx,org_row,org_col); --> |
326 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
327 |
|
|
*** <!-- int32 org_row; --> |
328 |
|
|
*** <!-- int32 org_col; --> |
329 |
aw0a |
1 |
*** Creates the given element and returns a pointer to it. The value is |
330 |
|
|
*** initially zero. |
331 |
|
|
*** It is *ASSUMED* that org_row |
332 |
|
|
*** and org_col are legal indices. May crash if they are not. |
333 |
|
|
*** If mtx_DEBUG is TRUE, then we will whine if the element already |
334 |
|
|
*** exists, but go ahead and create it anyway. |
335 |
|
|
**/ |
336 |
|
|
|
337 |
jds |
54 |
struct element_t *mtx_create_element_value(mtx_matrix_t mtx, |
338 |
|
|
int32 org_row, |
339 |
|
|
int32 org_col, |
340 |
|
|
real64 val); |
341 |
|
|
/**< |
342 |
|
|
*** <!-- mtx_create_element_value(mtx,org_row,org_col,val); --> |
343 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
344 |
|
|
*** <!-- int32 org_row; --> |
345 |
|
|
*** <!-- int32 org_col; --> |
346 |
|
|
*** <!-- real64 val; --> |
347 |
aw0a |
1 |
*** Creates the given element and returns a pointer to it. The value is |
348 |
|
|
*** initialzed to val. |
349 |
|
|
*** It is *ASSUMED* that org_row |
350 |
|
|
*** and org_col are legal indices. May crash if they are not. |
351 |
|
|
*** If mtx_DEBUG is TRUE, then we will whine if the element already |
352 |
|
|
*** exists, but go ahead and create it anyway. |
353 |
|
|
**/ |
354 |
|
|
|
355 |
jds |
54 |
/* ************************************************************************ *\ |
356 |
aw0a |
1 |
Element list traversals. No linear algebra programmer with an ounce of |
357 |
|
|
intelligence would ever need to use these in critical path functions. |
358 |
jds |
54 |
\* ************************************************************************ */ |
359 |
|
|
extern struct element_t *mtx_next_col(register struct element_t *elt, |
360 |
|
|
mtx_range_t *rng, |
361 |
|
|
int32 *tocur); |
362 |
ben.allan |
33 |
/**< |
363 |
jds |
54 |
*** <!-- enext = struct element_t *mtx_next_col(elt,rng,tocur); --> |
364 |
|
|
*** <!-- struct element_t *elt, *enext; --> |
365 |
|
|
*** <!-- mtx_range_t *rng; --> |
366 |
|
|
*** <!-- int32 *tocur; --> |
367 |
|
|
*** |
368 |
aw0a |
1 |
*** Returns the next element after elt that is in the range |
369 |
|
|
*** rng according to the permutation vector tocur given. May return NULL. |
370 |
|
|
**/ |
371 |
|
|
|
372 |
jds |
54 |
extern struct element_t *mtx_next_row(register struct element_t *elt, |
373 |
|
|
mtx_range_t *rng, |
374 |
|
|
int32 *tocur); |
375 |
ben.allan |
33 |
/**< |
376 |
jds |
54 |
*** <!-- enext = struct element_t *mtx_next_row(elt,rng,tocur); --> |
377 |
|
|
*** <!-- struct element_t *elt, *enext; --> |
378 |
|
|
*** <!-- mtx_range_t *rng; --> |
379 |
|
|
*** <!-- int32 *tocur; --> |
380 |
|
|
*** |
381 |
aw0a |
1 |
*** Returns the next element after elt that is in the range |
382 |
|
|
*** rng according to the permutation vector tocur given. May return NULL. |
383 |
|
|
**/ |
384 |
|
|
|
385 |
jds |
54 |
/* ************************************************************************ *\ |
386 |
aw0a |
1 |
Permutation memory management. |
387 |
jds |
54 |
\* ************************************************************************ */ |
388 |
|
|
extern int32 *mtx_alloc_perm(int32 cap); |
389 |
|
|
/**< |
390 |
|
|
*** <!-- p = mtx_alloc_perm(cap); --> |
391 |
|
|
*** <!-- int32 cap, *p; --> |
392 |
aw0a |
1 |
*** Allocates a permutation vector. The user need |
393 |
|
|
*** not concern himself with the -1st element, which does exist. |
394 |
|
|
**/ |
395 |
|
|
|
396 |
jds |
54 |
extern void mtx_copy_perm(int32 *tarperm, int32 *srcperm, int32 cap); |
397 |
|
|
/**< |
398 |
|
|
*** <!-- mtx_copy_perm(tarperm,srcperm,cap) --> |
399 |
|
|
*** <!-- int32 *tarperm; --> |
400 |
|
|
*** <!-- int32 *srcperm; --> |
401 |
|
|
*** <!-- int32 cap; --> |
402 |
aw0a |
1 |
*** Copies srcperm to tarperm given the capacity of srcperm. |
403 |
|
|
*** If tarperm was obtained from alloc_perm(), the -1 has already been copied. |
404 |
|
|
**/ |
405 |
|
|
|
406 |
jds |
54 |
extern void mtx_free_perm(int32 *perm); |
407 |
|
|
/**< |
408 |
|
|
*** <!-- mtx_free_perm(perm); --> |
409 |
|
|
*** <!-- int32 *perm; --> |
410 |
|
|
*** Frees resources used by a permutation vector. |
411 |
aw0a |
1 |
**/ |
412 |
|
|
|
413 |
jds |
54 |
/* ************************************************************************ *\ |
414 |
aw0a |
1 |
It is advantageous in an interactive system to introduce reusable |
415 |
|
|
memory and monitor its integrity rather than to repeatedly allocate |
416 |
|
|
and zero it. The following code accomplishes this for mtx. |
417 |
|
|
A null_vector is an array of objects (size s, length n) with value 0. |
418 |
|
|
This sort of memory management is needed because there is always the chance |
419 |
|
|
that a floating point exception could cause premature return of an mtx |
420 |
|
|
client. This way we have a safe place to store pointers to the memory |
421 |
|
|
even if the user's algorithm loses them. |
422 |
jds |
54 |
\* ************************************************************************ */ |
423 |
aw0a |
1 |
|
424 |
|
|
struct reusable_data_vector { |
425 |
ben.allan |
33 |
void *arr; /**< pointer to array of objects size entrysize */ |
426 |
|
|
int capacity; /**< number of object slots in array */ |
427 |
|
|
size_t entry_size; /**< size of slots */ |
428 |
|
|
int last_line; /**< line most recently associated with this structure, |
429 |
aw0a |
1 |
should be 0 if the array is not in use. */ |
430 |
|
|
}; |
431 |
|
|
|
432 |
|
|
extern struct reusable_data_vector |
433 |
jds |
54 |
g_mtx_null_index_data, /**< bunch of int32 */ |
434 |
|
|
g_mtx_null_sum_data, /**< bunch of mtx_value_t */ |
435 |
|
|
g_mtx_null_mark_data, /**< bunch of char */ |
436 |
|
|
g_mtx_null_vector_data, /**< bunch of element pointers */ |
437 |
|
|
g_mtx_null_col_vector_data, /**< bunch of element pointers */ |
438 |
|
|
g_mtx_null_row_vector_data; /**< bunch of element pointers */ |
439 |
aw0a |
1 |
|
440 |
jds |
54 |
/* OLD GROUP COMMENTS */ |
441 |
|
|
/* |
442 |
aw0a |
1 |
*** vec = mtx_null_vector(nptrs); |
443 |
|
|
*** vec = mtx_null_col_vector(nptrs); |
444 |
|
|
*** vec = mtx_null_row_vector(nptrs); |
445 |
|
|
*** marks = mtx_null_mark(nchar); |
446 |
|
|
*** sums = mtx_null_sum(nnums); |
447 |
|
|
*** indexes = mtx_null_index(ninds); |
448 |
|
|
*** |
449 |
|
|
*** struct element_t **vec; |
450 |
|
|
*** char *marks; |
451 |
|
|
*** real64 *sums; |
452 |
|
|
*** int32 *indexes; |
453 |
|
|
*** int32 nptrs, nchar, nnums, ninds; |
454 |
|
|
*** |
455 |
|
|
*** Returns an array of chars, elt pointers, indexes or numbers all NULL/0. |
456 |
|
|
*** We need these a lot, but seldom simultaneously, and we know generally |
457 |
|
|
*** how to rezero them when done with them. |
458 |
|
|
*** These functions should not be |
459 |
|
|
*** called again until the vector is re-NULLED and out of use. |
460 |
|
|
*** If we detect a double call, we will whine loudly, renull |
461 |
|
|
*** the array ourselves, and give it to you again. |
462 |
|
|
*** To avoid whining, call the corresponding release functions |
463 |
|
|
*** each time you are done with one of these vectors. |
464 |
|
|
*** |
465 |
|
|
*** In the event of insufficient memory (alloc failed) we will |
466 |
|
|
*** return NULL. If we return NULL, you needn't call the release function. |
467 |
|
|
*** |
468 |
|
|
*** mtx_null_vector_release(); |
469 |
|
|
*** mtx_null_col_vector_release(); |
470 |
|
|
*** mtx_null_row_vector_release(); |
471 |
|
|
*** mtx_null_mark_release(); |
472 |
|
|
*** mtx_null_sum_release(); |
473 |
|
|
*** mtx_null_index_release(); |
474 |
|
|
*** |
475 |
|
|
*** These are a memory reuse promoter. |
476 |
|
|
*** Calling with cap==0 frees any memory in use. |
477 |
|
|
*** Clientlists -- PLEASE KEEP THIS UP TO DATE -- |
478 |
|
|
*** mtx_null_vector: |
479 |
|
|
*** expand_row,expand_col, mtx_assemble |
480 |
|
|
*** mtx_householder_transform |
481 |
|
|
*** mtx_null_row_vector: |
482 |
|
|
*** expand_row_series |
483 |
|
|
*** mtx_null_col_vector: |
484 |
|
|
*** expand_col_series |
485 |
|
|
*** mtx_null_mark: |
486 |
|
|
*** mtx_householder_transform |
487 |
|
|
*** mtx_null_sum: |
488 |
|
|
*** mtx_householder_transform |
489 |
|
|
*** mtx_null_index: |
490 |
|
|
*** mtx_householder_transform |
491 |
|
|
**/ |
492 |
jds |
54 |
#define mtx_null_vector(cap) \ |
493 |
|
|
((struct element_t **)mtx_null_vector_f(cap,__LINE__,__FILE__, \ |
494 |
aw0a |
1 |
&g_mtx_null_vector_data,"null_vector")) |
495 |
jds |
54 |
/**< |
496 |
|
|
* Returns an array of elt pointers all NULL/0. |
497 |
|
|
* This function should not be called again until the vector is |
498 |
|
|
* re-NULLED and out of use. If we detect a double call, we will |
499 |
|
|
* whine loudly, renull the array ourselves, and give it to you again. |
500 |
|
|
* To avoid whining, call mtx_null_vector_release() |
501 |
|
|
* each time you are done with the returned vector.<br><br> |
502 |
|
|
* |
503 |
|
|
* In the event of insufficient memory (alloc failed) we will |
504 |
|
|
* return NULL. If we return NULL, you needn't call the release function. |
505 |
|
|
* @param cap int32, the capacity of the matrix (0 to free memory). |
506 |
|
|
* @return No return value. |
507 |
|
|
* @see mtx_null_vector_f() |
508 |
|
|
*/ |
509 |
|
|
#define mtx_null_row_vector(cap) \ |
510 |
|
|
((struct element_t **)mtx_null_vector_f(cap,__LINE__,__FILE__, \ |
511 |
aw0a |
1 |
&g_mtx_null_row_vector_data,\ |
512 |
|
|
"null_row_vector")) |
513 |
jds |
54 |
/**< |
514 |
|
|
* Returns an array of elt pointers all NULL/0. |
515 |
|
|
* See mtx_null_vector() for more details. |
516 |
|
|
* @param cap int32, the capacity of the matrix (0 to free memory). |
517 |
|
|
* @return No return value. |
518 |
|
|
* @see mtx_null_vector_f() |
519 |
|
|
*/ |
520 |
|
|
#define mtx_null_col_vector(cap) \ |
521 |
|
|
((struct element_t **)mtx_null_vector_f(cap,__LINE__,__FILE__, \ |
522 |
aw0a |
1 |
&g_mtx_null_col_vector_data,\ |
523 |
|
|
"null_col_vector")) |
524 |
jds |
54 |
/**< |
525 |
|
|
* Returns an array of elt pointers all NULL/0. |
526 |
|
|
* See mtx_null_vector() for more details. |
527 |
|
|
* @param cap int32, the capacity of the matrix (0 to free memory). |
528 |
|
|
* @return No return value. |
529 |
|
|
* @see mtx_null_vector_f() |
530 |
|
|
*/ |
531 |
|
|
#define mtx_null_mark(cap) \ |
532 |
|
|
((char *)mtx_null_vector_f(cap,__LINE__,__FILE__, \ |
533 |
aw0a |
1 |
&g_mtx_null_mark_data,"null_mark")) |
534 |
jds |
54 |
/**< |
535 |
|
|
* Returns an array of chars all NULL/0. |
536 |
|
|
* This function should not be called again until the vector is |
537 |
|
|
* re-NULLED and out of use. If we detect a double call, we will |
538 |
|
|
* whine loudly, renull the array ourselves, and give it to you again. |
539 |
|
|
* To avoid whining, call mtx_null_mark_release() |
540 |
|
|
* each time you are done with the returned vector.<br><br> |
541 |
|
|
* |
542 |
|
|
* In the event of insufficient memory (alloc failed) we will |
543 |
|
|
* return NULL. If we return NULL, you needn't call the release function. |
544 |
|
|
* @param cap int32, the capacity of the array (0 to free memory). |
545 |
|
|
* @return No return value. |
546 |
|
|
* @see mtx_null_vector_f() |
547 |
|
|
*/ |
548 |
|
|
#define mtx_null_sum(cap) \ |
549 |
|
|
((real64 *)mtx_null_vector_f(cap,__LINE__,__FILE__, \ |
550 |
aw0a |
1 |
&g_mtx_null_sum_data,"null_sum")) |
551 |
jds |
54 |
/**< |
552 |
|
|
* Returns an array of real64 numbers all NULL/0. |
553 |
|
|
* This function should not be called again until the vector is |
554 |
|
|
* re-NULLED and out of use. If we detect a double call, we will |
555 |
|
|
* whine loudly, renull the array ourselves, and give it to you again. |
556 |
|
|
* To avoid whining, call mtx_null_sum_release() |
557 |
|
|
* each time you are done with the returned array.<br><br> |
558 |
|
|
* |
559 |
|
|
* In the event of insufficient memory (alloc failed) we will |
560 |
|
|
* return NULL. If we return NULL, you needn't call the release function. |
561 |
|
|
* @param cap int32, the capacity of the array (0 to free memory). |
562 |
|
|
* @return No return value. |
563 |
|
|
* @see mtx_null_vector_f() |
564 |
|
|
*/ |
565 |
|
|
#define mtx_null_index(cap) \ |
566 |
|
|
((int32 *)mtx_null_vector_f(cap,__LINE__,__FILE__, \ |
567 |
aw0a |
1 |
&g_mtx_null_index_data,"null_index")) |
568 |
|
|
|
569 |
jds |
54 |
/**< |
570 |
|
|
* Returns an array of int32 indexes all NULL/0. |
571 |
|
|
* This function should not be called again until the vector is |
572 |
|
|
* re-NULLED and out of use. If we detect a double call, we will |
573 |
|
|
* whine loudly, renull the array ourselves, and give it to you again. |
574 |
|
|
* To avoid whining, call mtx_null_index_release() |
575 |
|
|
* each time you are done with the returned array.<br><br> |
576 |
|
|
* |
577 |
|
|
* In the event of insufficient memory (alloc failed) we will |
578 |
|
|
* return NULL. If we return NULL, you needn't call the release function. |
579 |
|
|
* @param cap int32, the capacity of the array (0 to free memory). |
580 |
|
|
* @return No return value. |
581 |
|
|
* @see mtx_null_vector_f() |
582 |
|
|
*/ |
583 |
aw0a |
1 |
#define mtx_null_vector_release() \ |
584 |
|
|
mtx_null_vector_release_f(__LINE__,__FILE__, \ |
585 |
|
|
&g_mtx_null_vector_data,"null_vector") |
586 |
jds |
54 |
/**< |
587 |
|
|
* Marks a vector as not in use, or whines if it wasn't. |
588 |
|
|
* @param None. |
589 |
|
|
* @return No return value. |
590 |
|
|
* @see mtx_null_vector_release_f() |
591 |
|
|
*/ |
592 |
aw0a |
1 |
#define mtx_null_col_vector_release() \ |
593 |
|
|
mtx_null_vector_release_f(__LINE__,__FILE__, \ |
594 |
|
|
&g_mtx_null_col_vector_data,"null_col_vector") |
595 |
jds |
54 |
/**< |
596 |
|
|
* Marks a vector as not in use, or whines if it wasn't. |
597 |
|
|
* @param None. |
598 |
|
|
* @return No return value. |
599 |
|
|
* @see mtx_null_vector_release_f() |
600 |
|
|
*/ |
601 |
aw0a |
1 |
#define mtx_null_row_vector_release() \ |
602 |
|
|
mtx_null_vector_release_f(__LINE__,__FILE__, \ |
603 |
|
|
&g_mtx_null_row_vector_data,"null_row_vector") |
604 |
jds |
54 |
/**< |
605 |
|
|
* Marks a vector as not in use, or whines if it wasn't. |
606 |
|
|
* @param None. |
607 |
|
|
* @return No return value. |
608 |
|
|
* @see mtx_null_vector_release_f() |
609 |
|
|
*/ |
610 |
aw0a |
1 |
#define mtx_null_mark_release() \ |
611 |
|
|
mtx_null_vector_release_f(__LINE__,__FILE__, \ |
612 |
|
|
&g_mtx_null_mark_data,"null_mark") |
613 |
jds |
54 |
/**< |
614 |
|
|
* Marks a char array as not in use, or whines if it wasn't. |
615 |
|
|
* @param None. |
616 |
|
|
* @return No return value. |
617 |
|
|
* @see mtx_null_vector_release_f() |
618 |
|
|
*/ |
619 |
aw0a |
1 |
#define mtx_null_sum_release() \ |
620 |
|
|
mtx_null_vector_release_f(__LINE__,__FILE__, \ |
621 |
|
|
&g_mtx_null_sum_data,"null_sum") |
622 |
jds |
54 |
/**< |
623 |
|
|
* Marks a number array as not in use, or whines if it wasn't. |
624 |
|
|
* @param None. |
625 |
|
|
* @return No return value. |
626 |
|
|
* @see mtx_null_vector_release_f() |
627 |
|
|
*/ |
628 |
aw0a |
1 |
#define mtx_null_index_release() \ |
629 |
|
|
mtx_null_vector_release_f(__LINE__,__FILE__, \ |
630 |
|
|
&g_mtx_null_index_data,"null_index") |
631 |
jds |
54 |
/**< |
632 |
|
|
* Marks an index array as not in use, or whines if it wasn't. |
633 |
|
|
* @param None. |
634 |
|
|
* @return No return value. |
635 |
|
|
* @see mtx_null_vector_release_f() |
636 |
|
|
*/ |
637 |
aw0a |
1 |
|
638 |
jds |
54 |
extern void *mtx_null_vector_f(int32 cap, int line, CONST char *file, |
639 |
|
|
struct reusable_data_vector *ptr, char *fn); |
640 |
|
|
/**< |
641 |
|
|
*** <!-- v = mtx_null_vector_f(cap,line,file, ptr,fn); --> |
642 |
|
|
*** <!-- int32 cap; --> |
643 |
|
|
*** <!-- int line; --> |
644 |
|
|
*** <!-- CONST char *file; --> |
645 |
|
|
*** <!-- struct reusable_data_vector *ptr; --> |
646 |
|
|
*** <!-- char *fn; --> |
647 |
|
|
*** |
648 |
|
|
*** Implementation function for macros generating vectors of NULL |
649 |
|
|
*** elements. This includes: |
650 |
|
|
*** - mtx_null_vector() |
651 |
|
|
*** - mtx_null_col_vector() |
652 |
|
|
*** - mtx_null_row_vector() |
653 |
|
|
*** - mtx_null_mark() |
654 |
|
|
*** - mtx_null_sum() |
655 |
|
|
*** - mtx_null_index() |
656 |
|
|
*** |
657 |
|
|
*** Do not call this function directly - use the appropriate macro |
658 |
|
|
*** instead. |
659 |
|
|
*** Returns a pointer to cap*ptr->entry_size bytes, which must be cast. |
660 |
|
|
*** The memory pointed at is believed to be zero, and will be if the |
661 |
|
|
*** user is properly rezeroing the vector before it is released. |
662 |
|
|
*** If insufficient memory is available, this whines and returns NULL. |
663 |
|
|
*** Calling this with cap==0 causes the reused memory to be deallocated and |
664 |
|
|
*** returns NULL. |
665 |
aw0a |
1 |
**/ |
666 |
|
|
|
667 |
jds |
54 |
extern void mtx_null_vector_release_f(int line, |
668 |
|
|
CONST char *file, |
669 |
|
|
struct reusable_data_vector *ptr, |
670 |
|
|
char *fn); |
671 |
|
|
/**< |
672 |
|
|
*** <!-- mtx_null_vector_release_f(line,file,ptr,fn); --> |
673 |
|
|
*** <!-- int line; --> |
674 |
|
|
*** <!-- CONST char *file, --> |
675 |
|
|
*** <!-- struct reusable_data_vector *ptr; --> |
676 |
|
|
*** <!-- char *fn --> |
677 |
aw0a |
1 |
*** |
678 |
jds |
54 |
*** Implementation function for macros releasing reusable vectors. |
679 |
|
|
*** This includes: |
680 |
|
|
*** - mtx_null_vector_release() |
681 |
|
|
*** - mtx_null_col_vector_release() |
682 |
|
|
*** - mtx_null_row_vector_release() |
683 |
|
|
*** - mtx_null_mark_release() |
684 |
|
|
*** - mtx_null_sum_release() |
685 |
|
|
*** - mtx_null_index_release() |
686 |
|
|
*** |
687 |
|
|
*** Do not call this function directly - use the appropriate macro |
688 |
|
|
*** instead. |
689 |
|
|
*** Marks a vector as not in use, or whines if it wasn't. |
690 |
|
|
*** Does no other checking. Uses line, file and fn in error reporting. |
691 |
aw0a |
1 |
**/ |
692 |
|
|
|
693 |
|
|
extern void mtx_reset_null_vectors(void); |
694 |
jds |
54 |
/**< |
695 |
aw0a |
1 |
*** This resets the reusable arrays of zeroes to zero in the event |
696 |
|
|
*** that they may have been corrupted. |
697 |
|
|
**/ |
698 |
|
|
|
699 |
jds |
54 |
/* |
700 |
|
|
* INTERNAL element vector operations of some utility. |
701 |
aw0a |
1 |
*/ |
702 |
|
|
|
703 |
jds |
54 |
extern struct element_t **mtx_expand_row(mtx_matrix_t mtx, int32 orgrow); |
704 |
|
|
/**< |
705 |
|
|
*** <!-- buf = mtx_expand_row(mtx,orgrow); --> |
706 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
707 |
|
|
*** <!-- int32 orgrow; --> |
708 |
|
|
*** <!-- struct element_t **buf; --> |
709 |
|
|
*** |
710 |
aw0a |
1 |
*** Expands the given row into an array of pointers, indexed on original |
711 |
|
|
*** col number. The array is obtained from mtx_null_vector(). |
712 |
|
|
*** Be sure to call mtx_null_vector_release() when done with the vector and |
713 |
|
|
*** you have rezeroed it. |
714 |
jds |
54 |
*** You cannot call this twice without releasing first or call |
715 |
|
|
*** mtx_expand_col(). |
716 |
aw0a |
1 |
**/ |
717 |
|
|
|
718 |
jds |
54 |
extern struct element_t **mtx_expand_col(mtx_matrix_t mtx, int32 orgcol); |
719 |
|
|
/**< |
720 |
|
|
*** <!-- buf = mtx_expand_col(mtx,orgcol); --> |
721 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
722 |
|
|
*** <!-- int32 orgcol; --> |
723 |
|
|
*** <!-- struct element_t **buf; --> |
724 |
|
|
*** |
725 |
aw0a |
1 |
*** Expands the given col into an array of pointers, indexed on original |
726 |
|
|
*** row number. The array is obtained from mtx_null_vector(). |
727 |
|
|
*** Be sure to call mtx_null_vector_release() when done with the vector and |
728 |
|
|
*** you have rezeroed it. |
729 |
jds |
54 |
*** You cannot call this twice without releasing first or call |
730 |
|
|
*** mtx_expand_row(). |
731 |
aw0a |
1 |
**/ |
732 |
|
|
|
733 |
jds |
54 |
extern void mtx_renull_using_row(mtx_matrix_t mtx, |
734 |
|
|
int32 orgrow, |
735 |
|
|
struct element_t **arr); |
736 |
|
|
/**< |
737 |
|
|
*** <!-- mtx_renull_using_row(mtx,orgrow,arr) --> |
738 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
739 |
|
|
*** <!-- int32 orgrow; --> |
740 |
|
|
*** <!-- struct element_t **arr; --> |
741 |
aw0a |
1 |
*** |
742 |
|
|
*** Makes arr NULLed again, assuming that the only non-NULL elements |
743 |
|
|
*** must correspond to original col numbers that exist in the given |
744 |
|
|
*** orgrow. |
745 |
|
|
**/ |
746 |
|
|
|
747 |
jds |
54 |
extern void mtx_renull_using_col(mtx_matrix_t mtx, |
748 |
|
|
int32 orgcol, |
749 |
|
|
struct element_t **arr); |
750 |
|
|
/**< |
751 |
|
|
*** <!-- mtx_renull_using_row(mtx,orgcol,arr); --> |
752 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
753 |
|
|
*** <!-- int32 orgcol; --> |
754 |
|
|
*** <!-- struct element_t **arr; --> |
755 |
aw0a |
1 |
*** |
756 |
|
|
*** Makes arr NULLed again, assuming that the only non-NULL elements |
757 |
|
|
*** must correspond to original row numbers that exist in the given |
758 |
|
|
*** orgcol. |
759 |
|
|
**/ |
760 |
|
|
|
761 |
jds |
54 |
extern void mtx_renull_all(mtx_matrix_t mtx, struct element_t **arr); |
762 |
|
|
/**< |
763 |
|
|
*** <!-- mtx_renull_all(mtx,arr); --> |
764 |
|
|
*** <!-- mtx_matrix_t mtx; --> |
765 |
|
|
*** <!-- struct element_t **arr; --> |
766 |
aw0a |
1 |
*** |
767 |
|
|
*** Makes arr NULLed again, assuming it is size mtx->order. |
768 |
|
|
**/ |
769 |
|
|
|
770 |
jds |
54 |
#endif /* __MTX_INTERNAL_USE_ONLY_H__ */ |
771 |
|
|
#endif /* none of your business if you aren't mtx_*.c */ |
772 |
|
|
|