Parent Directory
|
Revision Log
Setting up web subdirectory in repository
1 | /* |
2 | * mtx: Ascend Sparse Matrix Package |
3 | * by Karl Michael Westerberg |
4 | * Created: 2/6/90 |
5 | * Version: $Revision: 1.20 $ |
6 | * Version control file: $RCSfile: mtx_basic.c,v $ |
7 | * Date last modified: $Date: 2000/01/25 02:27:07 $ |
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 Kirk Andre Abbott, Benjamin Andrew Allan |
16 | * Copyright (C) 1996 Benjamin Andrew Allan |
17 | * |
18 | * The SLV solver is free software; you can redistribute |
19 | * it and/or modify it under the terms of the GNU General Public License as |
20 | * published by the Free Software Foundation; either version 2 of the |
21 | * License, or (at your option) any later version. |
22 | * |
23 | * The SLV solver is distributed in hope that it will be |
24 | * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of |
25 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
26 | * General Public License for more details. |
27 | * |
28 | * You should have received a copy of the GNU General Public License along with |
29 | * the program; if not, write to the Free Software Foundation, Inc., 675 |
30 | * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING. |
31 | * COPYING is found in ../compiler. |
32 | */ |
33 | |
34 | #include <math.h> |
35 | #include "utilities/ascConfig.h" |
36 | #include "utilities/ascMalloc.h" |
37 | #include "utilities/mem.h" |
38 | #include "solver/mtx.h" |
39 | /* grab our private parts */ |
40 | #define __MTX_C_SEEN__ |
41 | #include "solver/mtx_use_only.h" |
42 | |
43 | /* this returns an error count */ |
44 | int super_check_matrix( mtx_matrix_t mtx) |
45 | /** |
46 | *** Really check the matrix. |
47 | **/ |
48 | { |
49 | int32 ndx,errcnt=0; |
50 | int32 rowc,colc; |
51 | |
52 | /* Test consistency of permutation arrays */ |
53 | for( ndx=ZERO ; ndx < mtx->capacity ; ++ndx ) { |
54 | if( mtx->perm.row.org_to_cur[mtx->perm.row.cur_to_org[ndx]] != ndx ) { |
55 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
56 | FPRINTF(g_mtxerr," Permutation violation in row %d.\n", ndx); |
57 | errcnt++; |
58 | } |
59 | if( mtx->perm.col.org_to_cur[mtx->perm.col.cur_to_org[ndx]] != ndx ) { |
60 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
61 | FPRINTF(g_mtxerr," Permutation violation in col %d.\n",ndx); |
62 | errcnt++; |
63 | } |
64 | } |
65 | |
66 | if( mtx->order > mtx->capacity ) { |
67 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
68 | FPRINTF(g_mtxerr," Capacity %d is less than order %d\n", |
69 | mtx->capacity,mtx->order); |
70 | errcnt++; |
71 | } |
72 | |
73 | /* Make sure rows and columns which don't exist are blank */ |
74 | for( ndx = mtx->order ; ndx < mtx->capacity ; ++ndx ) { |
75 | int32 org_row,org_col; |
76 | org_row = mtx->perm.row.cur_to_org[ndx]; |
77 | org_col = mtx->perm.col.cur_to_org[ndx]; |
78 | if( NOTNULL(mtx->hdr.row[org_row]) ) { |
79 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
80 | FPRINTF(g_mtxerr," Non-zeros found in non-existent row %d.\n",ndx); |
81 | errcnt++; |
82 | } |
83 | if( NOTNULL(mtx->hdr.col[org_col]) ) { |
84 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
85 | FPRINTF(g_mtxerr," Non-zeros found in non-existent col %d.\n",ndx); |
86 | errcnt++; |
87 | } |
88 | } |
89 | |
90 | /* Verify fishnet structure */ |
91 | for( ndx=rowc=colc=ZERO ; ndx < mtx->capacity ; ++ndx ) { |
92 | struct element_t *elt; |
93 | int32 org_row,org_col; |
94 | org_row = mtx->perm.row.cur_to_org[ndx]; |
95 | org_col = mtx->perm.col.cur_to_org[ndx]; |
96 | for( elt = mtx->hdr.row[org_row] ; NOTNULL(elt) ; elt = elt->next.col ) { |
97 | if( elt->row != mtx_NONE ) { |
98 | ++rowc; |
99 | if( elt->row != org_row ) { |
100 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
101 | FPRINTF(g_mtxerr," Element mismatch in row %d.\n", ndx); |
102 | errcnt++; |
103 | } |
104 | } else { |
105 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
106 | FPRINTF(g_mtxerr," Disclaimed element in row %d.\n", ndx); |
107 | errcnt++; |
108 | } |
109 | } |
110 | for( elt = mtx->hdr.col[org_col] ; NOTNULL(elt) ; elt = elt->next.row ) { |
111 | if( elt->col != mtx_NONE ) { |
112 | ++colc; |
113 | if( elt->col != org_col ) { |
114 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
115 | FPRINTF(g_mtxerr," Element mismatch in col %d.\n", ndx); |
116 | errcnt++; |
117 | } |
118 | } else { |
119 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
120 | FPRINTF(g_mtxerr," Disclaimed element in col %d.\n", ndx); |
121 | errcnt++; |
122 | } |
123 | } |
124 | } |
125 | if( rowc != colc ) { |
126 | FPRINTF(g_mtxerr,"ERROR: (mtx) super_check_matrix\n"); |
127 | FPRINTF(g_mtxerr," Non-zero discrepancy, %d by row, %d by col.\n", |
128 | rowc, colc); |
129 | errcnt++; |
130 | } |
131 | return errcnt; |
132 | } |
133 | |
134 | boolean check_matrix( mtx_matrix_t mtx, char *file, int line) |
135 | /** |
136 | *** Checks the integrity flag of the matrix. |
137 | **/ |
138 | { |
139 | |
140 | if( ISNULL(mtx) ) { |
141 | FPRINTF(g_mtxerr,"ERROR: (mtx) check_matrix\n"); |
142 | FPRINTF(g_mtxerr," NULL matrix in %s line %d.\n",file,line); |
143 | return 0; |
144 | } |
145 | |
146 | switch( mtx->integrity ) { |
147 | case OK: |
148 | break; |
149 | case DESTROYED: |
150 | FPRINTF(g_mtxerr,"ERROR: (mtx) check_matrix\n"); |
151 | FPRINTF(g_mtxerr, |
152 | " Matrix deceased found in %s line %d.\n",file,line); |
153 | return 0; |
154 | default: |
155 | FPRINTF(g_mtxerr,"ERROR: (mtx) check_matrix\n"); |
156 | FPRINTF(g_mtxerr," Matrix garbage found in %s line %d.\n", |
157 | file,line); |
158 | return 0; |
159 | } |
160 | |
161 | #if MTX_DEBUG |
162 | super_check_matrix(mtx); |
163 | #endif |
164 | |
165 | if (ISSLAVE(mtx)) { |
166 | return (check_matrix(mtx->master,file,line)); |
167 | } |
168 | return 1; |
169 | } |
170 | |
171 | /* this function checks the consistency of a sparse as best it can. |
172 | returns FALSE if something wierd. */ |
173 | boolean check_sparse(const mtx_sparse_t * const sp, char *file, int line) |
174 | { |
175 | if ( ISNULL(sp) || |
176 | sp->len > sp->cap || |
177 | ( sp->cap > 0 && ( ISNULL(sp->idata) || ISNULL(sp->data) ) ) |
178 | ) { |
179 | FPRINTF(g_mtxerr,"ERROR: (mtx) check_sparse - Inconsistent or\n"); |
180 | FPRINTF(g_mtxerr," NULL sparse in %s line %d.\n",file,line); |
181 | return 0; |
182 | } |
183 | return 1; |
184 | } |
185 | |
186 | #define free_unless_null(ptr) if( NOTNULL(ptr) ) ascfree(ptr) |
187 | /* |
188 | static void free_unless_null(POINTER ptr) |
189 | { |
190 | if( NOTNULL(ptr) ) |
191 | ascfree(ptr); |
192 | } |
193 | */ |
194 | |
195 | void mtx_zero_int32(int32 *data, int len) |
196 | { |
197 | int i; |
198 | if (data==NULL) return; |
199 | for (i=0; i<len; i++) data[i]=0; |
200 | } |
201 | |
202 | void mtx_zero_real64(real64 *data, int len) |
203 | { |
204 | int i; |
205 | if (data==NULL) return; |
206 | for (i=0; i<len; i++) data[i]=0.0; |
207 | } |
208 | |
209 | void mtx_zero_ptr(void **data, int len) |
210 | { |
211 | int i; |
212 | if (data==NULL) return; |
213 | for (i=0; i<len; i++) data[i]=NULL; |
214 | } |
215 | |
216 | static mtx_matrix_t alloc_header(void) |
217 | /** |
218 | *** Allocates a matrix header and returns a pointer to it. |
219 | **/ |
220 | { |
221 | mtx_matrix_t mtx; |
222 | mtx = (mtx_matrix_t)ascmalloc( sizeof(struct mtx_header) ); |
223 | mtx->integrity = OK; |
224 | return(mtx); |
225 | } |
226 | |
227 | static void free_header(mtx_matrix_t mtx) |
228 | /** |
229 | *** Frees a matrix header. |
230 | **/ |
231 | { |
232 | mtx->integrity = DESTROYED; |
233 | ascfree(mtx); |
234 | } |
235 | |
236 | #ifdef DELETE_THIS_UNUSED_FUNCTION |
237 | /* not IN use? */ |
238 | static struct element_t **alloc_nz_header(int32 cap) |
239 | /** |
240 | *** Allocates a row or column header randomly initialized. |
241 | **/ |
242 | { |
243 | return(cap > 0 ? (struct element_t **) |
244 | ascmalloc(cap*sizeof(struct element_t *)) : NULL); |
245 | } |
246 | #endif /* DELETE_THIS_UNUSED_FUNCTION */ |
247 | |
248 | |
249 | static struct element_t **calloc_nz_header(int32 cap) |
250 | /** |
251 | *** Allocates a zeroed row or column header. |
252 | **/ |
253 | { |
254 | return(cap > 0 ? (struct element_t **) |
255 | asccalloc(cap,sizeof(struct element_t *)) : NULL); |
256 | } |
257 | |
258 | static void copy_nz_header(struct element_t **tarhdr, |
259 | struct element_t **srchdr, int32 cap) |
260 | /** |
261 | *** Copies srchdr to tarhdr given the capacity of srchdr. |
262 | **/ |
263 | { |
264 | mem_copy_cast(srchdr,tarhdr,cap*sizeof(struct element_t *)); |
265 | } |
266 | |
267 | #define free_nz_header(hdr) free_unless_null( (POINTER)(hdr) ) |
268 | /** |
269 | *** Frees a row or column header |
270 | **/ |
271 | |
272 | |
273 | /********************************************************************/ |
274 | /* |
275 | Any method of deallocating elements that doesn't use |
276 | delete_from_row/col or nuke_fishnet is DOOMED to die. |
277 | On ANY entry to delete_from_row or delete_from_col, |
278 | last_value_matrix should be already set to the mtx being deleted from. |
279 | What this really means is that before any of the following list of functions |
280 | is called, the last_value_matrix should be set: |
281 | |
282 | disclaim element |
283 | delete_from_row/col |
284 | blast_row/col |
285 | clean_row/col |
286 | |
287 | Note that the above functions call each other, with blast_/clean_ being |
288 | the top housekeeping/garbage collection pair. Users of delete, blast or |
289 | clean should set last_value_matrix at the scope where the matrix is known. |
290 | |
291 | Excuses for this kluge: |
292 | 1) mtx_value/set_value need it. |
293 | 2) the mem_store scheme REQUIRES it. |
294 | 3) it is faster than passing the mtx |
295 | up and down through the functions listed above, and we do a LOT of |
296 | these operations. |
297 | */ |
298 | static mtx_matrix_t last_value_matrix = NULL; |
299 | |
300 | static void disclaim_element(struct element_t *element) |
301 | /** |
302 | *** Indicates that one of the row or column (it doesn't matter) in which |
303 | *** this element appears will no longer point to this element (it has |
304 | *** already been removed from the list). Elements disclaimed once have |
305 | *** their row,col indices set to mtx_NONE. Elements disclaimed twice are |
306 | *** discarded. |
307 | *** NOTE WELL: last_value_mtx must be set before this function is called. |
308 | **/ |
309 | { |
310 | if( element->row == mtx_NONE ) { |
311 | mem_free_element(last_value_matrix->ms,(void *)element); |
312 | } else { |
313 | element->row = element->col = mtx_NONE; |
314 | } |
315 | } |
316 | |
317 | static void delete_from_row(struct element_t **link) |
318 | /** |
319 | *** Deletes the element from the row it is in, given a pointer to |
320 | *** the row link which leads to this element. |
321 | *** On return the pointer pointed to by link points to the new next element. |
322 | *** NOTE WELL: last_value_mtx must be set before this function is called. |
323 | **/ |
324 | { |
325 | struct element_t *p; |
326 | p = *link; |
327 | *link = p->next.col; |
328 | disclaim_element(p); |
329 | /* conservatively cause mtx_set_value to forget */ |
330 | last_value_matrix->last_value = NULL; |
331 | } |
332 | |
333 | static void delete_from_col(struct element_t **link) |
334 | /** |
335 | *** Deletes the element from the col it is in, given a pointer to |
336 | *** the col link which leads to this element. |
337 | *** NOTE WELL: last_value_mtx must be set before this function is called. |
338 | **/ |
339 | { |
340 | struct element_t *p; |
341 | p = *link; |
342 | *link = p->next.row; |
343 | disclaim_element(p); |
344 | /* conservatively cause mtx_set_value to forget */ |
345 | last_value_matrix->last_value = NULL; |
346 | } |
347 | |
348 | static void nuke_fishnet(mtx_matrix_t mtx) |
349 | /** |
350 | *** Deletes every element even remotely associated with a matrix |
351 | *** and nulls the headers. Also nulls the headers of slaves |
352 | *** and cleans out their elements. |
353 | **/ |
354 | { |
355 | int32 i; |
356 | |
357 | if ( ISNULL(mtx) || ISNULL(mtx->hdr.row) || |
358 | ISNULL(mtx->hdr.col) || ISSLAVE(mtx)) { |
359 | FPRINTF(g_mtxerr,"WARNING: (mtx) nuke_fishnet called\n"); |
360 | FPRINTF(g_mtxerr," with bad matrix.\n"); |
361 | return; |
362 | } |
363 | if (mtx->capacity<1) return; |
364 | mtx->last_value = NULL; |
365 | mem_clear_store(mtx->ms); |
366 | zero(mtx->hdr.row,mtx->capacity,struct element_t *); |
367 | zero(mtx->hdr.col,mtx->capacity,struct element_t *); |
368 | for (i = 0; i < mtx->nslaves; i++) { |
369 | zero(mtx->slaves[i]->hdr.row,mtx->capacity,struct element_t *); |
370 | zero(mtx->slaves[i]->hdr.col,mtx->capacity,struct element_t *); |
371 | } |
372 | } |
373 | |
374 | static void blast_row( mtx_matrix_t mtx, int32 org, mtx_range_t *rng) |
375 | /** |
376 | *** Destroys all elements in the given row with (current) |
377 | *** col index in the given col range: if rng == mtx_ALL_COLS, |
378 | *** entire row is destroyed. |
379 | *** NOTE WELL: last_value_mtx must be set before this function is called. |
380 | **/ |
381 | { |
382 | struct element_t **link; |
383 | |
384 | link = &(mtx->hdr.row[org]); |
385 | if( rng == mtx_ALL_COLS ) { |
386 | while( NOTNULL(*link) ) |
387 | delete_from_row(link); |
388 | } else { |
389 | int32 *tocur = mtx->perm.col.org_to_cur; |
390 | int32 col; |
391 | while( NOTNULL(*link) ) { |
392 | col = (*link)->col; |
393 | if( col == mtx_NONE || in_range(rng,tocur[col]) ) |
394 | delete_from_row(link); |
395 | else |
396 | link = &((*link)->next.col); |
397 | } |
398 | } |
399 | } |
400 | |
401 | static void blast_col( mtx_matrix_t mtx, int32 org, mtx_range_t *rng) |
402 | /** |
403 | *** Destroys all elements in the given col with (current) |
404 | *** row index in the given row range: if rng == mtx_ALL_ROWS, |
405 | *** entire col is destroyed. |
406 | *** NOTE WELL: last_value_mtx must be set before this function is called. |
407 | **/ |
408 | { |
409 | struct element_t **link; |
410 | |
411 | link = &(mtx->hdr.col[org]); |
412 | if( rng == mtx_ALL_ROWS ) { |
413 | while( NOTNULL(*link) ) |
414 | delete_from_col(link); |
415 | } else { |
416 | int32 *tocur = mtx->perm.row.org_to_cur, row; |
417 | |
418 | while( NOTNULL(*link) ) { |
419 | row = (*link)->row; |
420 | if( row == mtx_NONE || in_range(rng,tocur[row]) ) |
421 | delete_from_col(link); |
422 | else |
423 | link = &((*link)->next.row); |
424 | } |
425 | } |
426 | } |
427 | |
428 | static void clean_row(mtx_matrix_t mtx, int32 org) |
429 | /** |
430 | *** Disclaims all elements in the given row which have already been |
431 | *** disclaimed once. |
432 | *** NOTE WELL: last_value_mtx must be set before this function is called. |
433 | **/ |
434 | { |
435 | struct element_t **link; |
436 | |
437 | link = &(mtx->hdr.row[org]); |
438 | while( NOTNULL(*link) ) |
439 | if( (*link)->row == mtx_NONE ) |
440 | delete_from_row(link); |
441 | else |
442 | link = &((*link)->next.col); |
443 | } |
444 | |
445 | static void clean_col(mtx_matrix_t mtx, int32 org) |
446 | /** |
447 | *** Disclaims all elements in the given col which have already been |
448 | *** disclaimed once. |
449 | *** NOTE WELL: last_value_mtx must be set before this function is called. |
450 | **/ |
451 | { |
452 | struct element_t **link; |
453 | |
454 | link = &(mtx->hdr.col[org]); |
455 | while( NOTNULL(*link) ) |
456 | if( (*link)->col == mtx_NONE ) |
457 | delete_from_col(link); |
458 | else |
459 | link = &((*link)->next.row); |
460 | } |
461 | /********************************************************************/ |
462 | |
463 | mtx_coord_t *mtx_coord(mtx_coord_t *coordp, int32 row, int32 col) |
464 | { |
465 | coordp->row = row; |
466 | coordp->col = col; |
467 | return(coordp); |
468 | } |
469 | |
470 | mtx_range_t *mtx_range(mtx_range_t *rangep, int32 low, int32 high) |
471 | { |
472 | rangep->low = low; |
473 | rangep->high = high; |
474 | return(rangep); |
475 | } |
476 | |
477 | mtx_region_t *mtx_region( mtx_region_t *regionp, int32 rowlow, int32 rowhigh, |
478 | int32 collow, int32 colhigh) |
479 | { |
480 | mtx_range(&(regionp->row),rowlow,rowhigh); |
481 | mtx_range(&(regionp->col),collow,colhigh); |
482 | return(regionp); |
483 | } |
484 | |
485 | static int g_mtx_debug_redirect = 0; |
486 | |
487 | void mtx_debug_redirect_freeze() { |
488 | g_mtx_debug_redirect = 1; |
489 | } |
490 | |
491 | static void mtx_redirectErrors(FILE *f) { |
492 | if (!g_mtx_debug_redirect) { |
493 | assert(f != NULL); |
494 | g_mtxerr = f; |
495 | } |
496 | } |
497 | |
498 | mtx_matrix_t mtx_create(void) |
499 | { |
500 | mtx_matrix_t mtx; |
501 | |
502 | mtx_redirectErrors(stderr); /* call mtx_debug_redirect_freeze() to bypass */ |
503 | |
504 | mtx = alloc_header(); |
505 | mtx->order = ZERO; |
506 | mtx->capacity = ZERO; |
507 | mtx->hdr.row = mtx->hdr.col = NULL; |
508 | mtx->perm.row.org_to_cur = mtx_alloc_perm(ZERO); |
509 | mtx->perm.row.cur_to_org = mtx_alloc_perm(ZERO); |
510 | mtx->perm.row.parity = EVEN; |
511 | mtx->perm.col.org_to_cur = mtx_alloc_perm(ZERO); |
512 | mtx->perm.col.cur_to_org = mtx_alloc_perm(ZERO); |
513 | mtx->perm.col.parity = EVEN; |
514 | mtx->perm.transpose = mtx_NONE; |
515 | mtx->data = |
516 | (struct structural_data_t *)ascmalloc(sizeof(struct structural_data_t)); |
517 | mtx->data->symbolic_rank = -1; |
518 | mtx->data->nblocks = -1; |
519 | mtx->data->block = NULL; |
520 | mtx->master = NULL; |
521 | mtx->nslaves = ZERO; |
522 | mtx->slaves = NULL; |
523 | mtx->last_value = NULL; |
524 | mtx->ms = mem_create_store(LENMAGIC, WIDTHMAGIC/sizeof(struct element_t) - 1, |
525 | sizeof(struct element_t),10,2048); |
526 | return(mtx); |
527 | } |
528 | |
529 | /** |
530 | ** Slave headers share the master's: |
531 | ** perm.anything, data->anything, ms->anything, order and capacity. |
532 | ** Unique to the slave are: |
533 | ** hdr.anything |
534 | ** last_value. |
535 | ** All slaves of a master appear as pointers in the masters slaves list |
536 | ** which is nslaves long. slaves do not have slaves. |
537 | **/ |
538 | mtx_matrix_t mtx_create_slave(mtx_matrix_t master) |
539 | { |
540 | mtx_matrix_t mtx, *mlist; |
541 | int32 newcnt; |
542 | mtx_redirectErrors(stderr); /* call mtx_debug_redirect_freeze() to bypass */ |
543 | if(!mtx_check_matrix(master) || ISSLAVE(master)) return NULL; |
544 | |
545 | mtx = alloc_header(); |
546 | if (ISNULL(mtx)) return mtx; |
547 | mtx->order = master->order; |
548 | mtx->capacity = master->capacity; |
549 | mtx->nslaves = ZERO; |
550 | mtx->hdr.row = calloc_nz_header(master->capacity); |
551 | if (ISNULL(mtx->hdr.row)) { |
552 | ascfree(mtx); |
553 | FPRINTF(g_mtxerr,"mtx_create_slave: Insufficient memory.\n"); |
554 | return NULL; |
555 | } |
556 | mtx->hdr.col = calloc_nz_header(master->capacity); |
557 | if (ISNULL(mtx->hdr.col)) { |
558 | ascfree(mtx->hdr.row); |
559 | ascfree(mtx); |
560 | FPRINTF(g_mtxerr,"mtx_create_slave: Insufficient memory.\n"); |
561 | return NULL; |
562 | } |
563 | mtx->perm.row.org_to_cur = master->perm.row.org_to_cur; |
564 | mtx->perm.row.cur_to_org = master->perm.row.cur_to_org; |
565 | mtx->perm.row.parity = EVEN; /* dont look at this again*/ |
566 | mtx->perm.col.org_to_cur = master->perm.col.org_to_cur; |
567 | mtx->perm.col.cur_to_org = master->perm.col.cur_to_org; |
568 | mtx->perm.col.parity = EVEN; /* dont look at this again*/ |
569 | mtx->perm.transpose = mtx_NONE; /* dont look at this again*/ |
570 | mtx->slaves = NULL; |
571 | mtx->data = master->data; |
572 | mtx->master = master; |
573 | mtx->last_value = NULL; |
574 | mtx->ms = master->ms; |
575 | newcnt = master->nslaves + 1; |
576 | if (newcnt == 1) { |
577 | mlist = (mtx_matrix_t *)ascmalloc(sizeof(mtx_matrix_t)); |
578 | } else { |
579 | mlist = |
580 | (mtx_matrix_t *)ascrealloc(master->slaves,newcnt*sizeof(mtx_matrix_t)); |
581 | } |
582 | if (ISNULL(mlist)) { |
583 | ascfree(mtx->hdr.row); |
584 | ascfree(mtx->hdr.col); |
585 | ascfree(mtx); |
586 | FPRINTF(g_mtxerr,"mtx_create_slave: Insufficient memory.\n"); |
587 | return NULL; |
588 | } else { |
589 | master->slaves = mlist; |
590 | master->slaves[master->nslaves] = mtx; |
591 | master->nslaves = newcnt; |
592 | } |
593 | return(mtx); |
594 | } |
595 | |
596 | static void destroy_slave(mtx_matrix_t master, mtx_matrix_t mtx) |
597 | { |
598 | int32 sid,i; |
599 | |
600 | if(!mtx_check_matrix(mtx) || !mtx_check_matrix(master)) { |
601 | FPRINTF(g_mtxerr, |
602 | "destroy_slave(mtx.c) called with corrupt slave or master\n"); |
603 | return; |
604 | } |
605 | if (ISSLAVE(master) || !ISSLAVE(mtx)) { |
606 | FPRINTF(g_mtxerr, |
607 | "destroy_slave(mtx.c) called with mismatched slave/master pair\n"); |
608 | return; |
609 | } |
610 | /* get slave index */ |
611 | sid = -1; |
612 | for (i=0; i<master->nslaves; i++) { |
613 | if (master->slaves[i] == mtx) { |
614 | sid = i; |
615 | } |
616 | } |
617 | if (sid < 0) { |
618 | FPRINTF(g_mtxerr, |
619 | "destroy_slave(mtx.c) called with mismatched slave/master pair\n"); |
620 | return; |
621 | } |
622 | /* move the rest of the slaves up the list */ |
623 | for (i = sid; i < master->nslaves -1;) { |
624 | master->slaves[i] = master->slaves[++i]; |
625 | } |
626 | (master->nslaves)--; |
627 | /* we will not realloc smaller. */ |
628 | if (master->nslaves == 0 && NOTNULL(master->slaves)) { |
629 | ascfree(master->slaves); |
630 | master->slaves = NULL; |
631 | } |
632 | mtx_clear_region(mtx,mtx_ENTIRE_MATRIX); |
633 | free_nz_header(mtx->hdr.row); |
634 | free_nz_header(mtx->hdr.col); |
635 | free_header(mtx); |
636 | } |
637 | |
638 | void mtx_destroy(mtx_matrix_t mtx) |
639 | { |
640 | int32 i; |
641 | if(!mtx_check_matrix(mtx)) return; |
642 | if (ISSLAVE(mtx)) { |
643 | destroy_slave(mtx->master,mtx); |
644 | return; |
645 | } |
646 | |
647 | for (i=0; i<mtx->nslaves; i++) { |
648 | if (mtx_check_matrix(mtx->slaves[i])) { |
649 | free_nz_header(mtx->slaves[i]->hdr.row); |
650 | free_nz_header(mtx->slaves[i]->hdr.col); |
651 | free_header(mtx->slaves[i]); |
652 | mtx->slaves[i] = NULL; |
653 | } else { |
654 | FPRINTF(g_mtxerr, |
655 | "mtx_destroy: Corrupt slave found while destroying master.\n"); |
656 | FPRINTF(g_mtxerr, |
657 | " Slave %d being abandoned.\n",i); |
658 | } |
659 | } |
660 | mtx->nslaves = 0; |
661 | |
662 | mtx_clear_region(mtx,mtx_ENTIRE_MATRIX); |
663 | /* mtx_check_matrix is called again if MTX_DEBUG*/ |
664 | /** |
665 | *** The fishnet structure is |
666 | *** destroyed, just leaving |
667 | *** the headers and maybe ms. |
668 | **/ |
669 | mem_destroy_store(mtx->ms); |
670 | |
671 | free_nz_header(mtx->hdr.row); |
672 | free_nz_header(mtx->hdr.col); |
673 | |
674 | mtx_free_perm(mtx->perm.row.org_to_cur); |
675 | mtx_free_perm(mtx->perm.row.cur_to_org); |
676 | mtx_free_perm(mtx->perm.col.org_to_cur); |
677 | mtx_free_perm(mtx->perm.col.cur_to_org); |
678 | |
679 | if( NOTNULL(mtx->data->block) ) { |
680 | ascfree(mtx->data->block); |
681 | } |
682 | ascfree(mtx->data); |
683 | |
684 | free_header(mtx); |
685 | } |
686 | |
687 | mtx_sparse_t *mtx_create_sparse(int32 cap) { |
688 | mtx_sparse_t *ret; |
689 | ret = (mtx_sparse_t *)ascmalloc(sizeof(mtx_sparse_t)); |
690 | if (ISNULL(ret)) { |
691 | FPRINTF(g_mtxerr,"ERROR: (mtx_create_sparse) Insufficient memory.\n"); |
692 | return ret; |
693 | } |
694 | ret->data = (real64 *)ascmalloc(cap*sizeof(real64)); |
695 | if (ISNULL(ret->data)) { |
696 | FPRINTF(g_mtxerr,"ERROR: (mtx_create_sparse) Insufficient memory.\n"); |
697 | ascfree(ret); |
698 | return NULL; |
699 | } |
700 | ret->idata = (int32 *)ascmalloc(cap*sizeof(int32)); |
701 | if (ISNULL(ret->idata)) { |
702 | FPRINTF(g_mtxerr,"ERROR: (mtx_create_sparse) Insufficient memory.\n"); |
703 | ascfree(ret->data); |
704 | ascfree(ret); |
705 | return NULL; |
706 | } |
707 | ret->cap = cap; |
708 | ret->len = 0; |
709 | return ret; |
710 | } |
711 | |
712 | void mtx_destroy_sparse(mtx_sparse_t *ret) |
713 | { |
714 | if (ISNULL(ret)) return; |
715 | if (NOTNULL(ret->idata)) ascfree(ret->idata); |
716 | if (NOTNULL(ret->data)) ascfree(ret->data); |
717 | ascfree(ret); |
718 | } |
719 | |
720 | void mtx_destroy_blocklist(mtx_block_t *bl) |
721 | { |
722 | if (ISNULL(bl)) return; |
723 | if (NOTNULL(bl->block) && bl->nblocks>0) ascfree(bl->block); |
724 | ascfree(bl); |
725 | } |
726 | |
727 | mtx_matrix_t mtx_duplicate_region(mtx_matrix_t mtx, mtx_region_t *reg, |
728 | real64 drop) |
729 | { |
730 | mtx_matrix_t new; |
731 | int32 org_row; |
732 | struct element_t *elt; |
733 | |
734 | if (!mtx_check_matrix(mtx) || |
735 | (mtx->master != NULL && !mtx_check_matrix(mtx->master))) { |
736 | return NULL; |
737 | } |
738 | drop = fabs(drop); |
739 | if (ISSLAVE(mtx)) { |
740 | new = mtx_create_slave(mtx->master); |
741 | } else { |
742 | new = mtx_create_slave(mtx); |
743 | } |
744 | if (new == NULL) { |
745 | FPRINTF(g_mtxerr,"ERROR: (mtx_duplicate_region) Insufficient memory.\n"); |
746 | return NULL; |
747 | } |
748 | if (reg==mtx_ENTIRE_MATRIX){ |
749 | if (drop >0.0) { /* copy with drop tolerance */ |
750 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
751 | for( elt = mtx->hdr.row[org_row] ; |
752 | NOTNULL(elt) ; elt = elt->next.col ) { |
753 | if( elt->row != mtx_NONE && fabs(elt->value) > drop) { |
754 | mtx_create_element_value(new,elt->row,elt->col,elt->value); |
755 | } |
756 | } |
757 | } |
758 | } else { |
759 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
760 | for( elt = mtx->hdr.row[org_row] ; |
761 | NOTNULL(elt) ; elt = elt->next.col ) { |
762 | if( elt->row != mtx_NONE ) { |
763 | mtx_create_element_value(new,elt->row,elt->col,elt->value); |
764 | } |
765 | } |
766 | } |
767 | } |
768 | } else { |
769 | int32 *rowcur, *colcur, row,col; |
770 | mtx_range_t *colrng,*rowrng; |
771 | rowcur=mtx->perm.row.org_to_cur; |
772 | colcur=mtx->perm.col.org_to_cur; |
773 | rowrng=&(reg->row); |
774 | colrng=&(reg->col); |
775 | if (drop >0.0) { /* copy with drop tolerance */ |
776 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
777 | row=rowcur[org_row]; |
778 | if (in_range(rowrng,row)) { |
779 | for( elt = mtx->hdr.row[org_row] ; |
780 | NOTNULL(elt) ; elt = elt->next.col ) { |
781 | col=colcur[elt->col]; |
782 | if( in_range(colrng,col) && elt->row != mtx_NONE && |
783 | fabs(elt->value) > drop ) { |
784 | /* don't copy^bad^elements */ |
785 | mtx_create_element_value(new,elt->row,elt->col,elt->value); |
786 | } |
787 | } |
788 | } |
789 | } |
790 | } else { |
791 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
792 | row=rowcur[org_row]; |
793 | if (in_range(rowrng,row)) { |
794 | for( elt = mtx->hdr.row[org_row] ; |
795 | NOTNULL(elt) ; elt = elt->next.col ) { |
796 | col=colcur[elt->col]; |
797 | if( in_range(colrng,col) && elt->row != mtx_NONE ) { |
798 | /* don't copy^bad^elements */ |
799 | mtx_create_element_value(new,elt->row,elt->col,elt->value); |
800 | } |
801 | } |
802 | } |
803 | } |
804 | } |
805 | } |
806 | |
807 | return new; |
808 | } |
809 | |
810 | /* WARNING: this function assumes sufficient memory is available. |
811 | * it needs to be check for null returns in allocs and doesn't. |
812 | */ |
813 | mtx_matrix_t mtx_copy_options(mtx_matrix_t mtx, boolean blocks, |
814 | boolean incidence, mtx_region_t *reg, |
815 | real64 drop) |
816 | { |
817 | mtx_matrix_t copy; |
818 | int32 org_row; |
819 | struct element_t *elt; |
820 | int32 bnum; |
821 | |
822 | if (!mtx_check_matrix(mtx)) return NULL; |
823 | drop = fabs(drop); |
824 | |
825 | /* create same size matrix */ |
826 | copy = alloc_header(); |
827 | copy->order = mtx->order; |
828 | copy->capacity = mtx->capacity; |
829 | copy->master = NULL; |
830 | copy->nslaves = ZERO; |
831 | copy->slaves = NULL; |
832 | copy->last_value = NULL; |
833 | |
834 | /* copy headers and fishnet */ |
835 | copy->hdr.row = calloc_nz_header(copy->capacity); |
836 | copy->hdr.col = calloc_nz_header(copy->capacity); |
837 | if (incidence) { |
838 | struct mem_statistics stat; |
839 | int s_expool; |
840 | mem_get_stats(&stat,mtx->ms); |
841 | s_expool = MAX(2048,stat.elt_total/stat.str_wid); |
842 | copy->ms = mem_create_store(LENMAGIC, |
843 | WIDTHMAGIC/sizeof(struct element_t) - 1, |
844 | sizeof(struct element_t), |
845 | 10,s_expool-LENMAGIC); |
846 | /* copy of a slave or master matrix will end up with an |
847 | * initial mem_store that is the size cumulative of the |
848 | * master and all its slaves since they share the ms. |
849 | */ |
850 | if (reg==mtx_ENTIRE_MATRIX){ |
851 | if (drop >0.0) { /* copy with drop tolerance */ |
852 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
853 | for( elt = mtx->hdr.row[org_row] ; |
854 | NOTNULL(elt) ; elt = elt->next.col ) { |
855 | if( elt->row != mtx_NONE && fabs(elt->value) > drop) { |
856 | mtx_create_element_value(copy,elt->row,elt->col,elt->value); |
857 | } |
858 | } |
859 | } |
860 | } else { |
861 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
862 | for( elt = mtx->hdr.row[org_row] ; |
863 | NOTNULL(elt) ; elt = elt->next.col ) { |
864 | if( elt->row != mtx_NONE ) { |
865 | mtx_create_element_value(copy,elt->row,elt->col,elt->value); |
866 | } |
867 | } |
868 | } |
869 | } |
870 | } else { |
871 | int32 *rowcur, *colcur, row,col; |
872 | mtx_range_t *colrng,*rowrng; |
873 | rowcur=mtx->perm.row.org_to_cur; |
874 | colcur=mtx->perm.col.org_to_cur; |
875 | rowrng=&(reg->row); |
876 | colrng=&(reg->col); |
877 | if (drop >0.0) { /* copy with drop tolerance */ |
878 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
879 | row=rowcur[org_row]; |
880 | if (in_range(rowrng,row)) { |
881 | for( elt = mtx->hdr.row[org_row] ; |
882 | NOTNULL(elt) ; elt = elt->next.col ) { |
883 | col=colcur[elt->col]; |
884 | if( in_range(colrng,col) && elt->row != mtx_NONE && |
885 | fabs(elt->value) > drop ) { |
886 | /* don't copy^bad^elements */ |
887 | mtx_create_element_value(copy,elt->row,elt->col,elt->value); |
888 | } |
889 | } |
890 | } |
891 | } |
892 | } else { |
893 | for( org_row=0 ; org_row < mtx->order ; ++org_row ) { |
894 | row=rowcur[org_row]; |
895 | if (in_range(rowrng,row)) { |
896 | for( elt = mtx->hdr.row[org_row] ; |
897 | NOTNULL(elt) ; elt = elt->next.col ) { |
898 | col=colcur[elt->col]; |
899 | if( in_range(colrng,col) && elt->row != mtx_NONE ) { |
900 | /* don't copy^bad^elements */ |
901 | mtx_create_element_value(copy,elt->row,elt->col,elt->value); |
902 | } |
903 | } |
904 | } |
905 | } |
906 | } |
907 | } |
908 | } else { |
909 | copy->ms = mem_create_store(LENMAGIC, |
910 | WIDTHMAGIC/sizeof(struct element_t) - 1, |
911 | sizeof(struct element_t),10,2048); |
912 | /* copy of a slave or master matrix will end up with an |
913 | initial mem_store that is the default size */ |
914 | } |
915 | |
916 | /* copy permutation information */ |
917 | copy->perm.row.org_to_cur = mtx_alloc_perm(copy->capacity); |
918 | mtx_copy_perm(copy->perm.row.org_to_cur, |
919 | mtx->perm.row.org_to_cur,copy->capacity); |
920 | copy->perm.row.cur_to_org = mtx_alloc_perm(copy->capacity); |
921 | mtx_copy_perm(copy->perm.row.cur_to_org, |
922 | mtx->perm.row.cur_to_org,copy->capacity); |
923 | copy->perm.row.parity = mtx_row_parity(mtx); |
924 | |
925 | copy->perm.col.org_to_cur = mtx_alloc_perm(copy->capacity); |
926 | mtx_copy_perm(copy->perm.col.org_to_cur, |
927 | mtx->perm.col.org_to_cur,copy->capacity); |
928 | copy->perm.col.cur_to_org = mtx_alloc_perm(copy->capacity); |
929 | mtx_copy_perm(copy->perm.col.cur_to_org, |
930 | mtx->perm.col.cur_to_org,copy->capacity); |
931 | copy->perm.col.parity = mtx_col_parity(mtx); |
932 | |
933 | /* copy (or not) block data */ |
934 | copy->data = |
935 | (struct structural_data_t *)ascmalloc(sizeof(struct structural_data_t)); |
936 | if (blocks) { |
937 | copy->data->symbolic_rank = mtx->data->symbolic_rank; |
938 | copy->data->nblocks = mtx->data->nblocks; |
939 | copy->data->block = mtx->data->nblocks > 0 ? (mtx_region_t *) |
940 | ascmalloc( mtx->data->nblocks*sizeof(mtx_region_t) ) : NULL; |
941 | for( bnum=0; bnum < mtx->data->nblocks; bnum++ ) { |
942 | copy->data->block[bnum].row.low = mtx->data->block[bnum].row.low; |
943 | copy->data->block[bnum].row.high = mtx->data->block[bnum].row.high; |
944 | copy->data->block[bnum].col.low = mtx->data->block[bnum].col.low; |
945 | copy->data->block[bnum].col.high = mtx->data->block[bnum].col.high; |
946 | } |
947 | } else { |
948 | copy->data->symbolic_rank=0; |
949 | copy->data->nblocks=0; |
950 | copy->data->block =NULL; |
951 | } |
952 | |
953 | return(copy); |
954 | } |
955 | |
956 | int32 mtx_order( mtx_matrix_t mtx) |
957 | { |
958 | if (!mtx_check_matrix(mtx)) return -1; |
959 | return(mtx->order); |
960 | } |
961 | |
962 | int32 mtx_capacity( mtx_matrix_t mtx) |
963 | { |
964 | if (!mtx_check_matrix(mtx)) return -1; |
965 | return(mtx->capacity); |
966 | } |
967 | |
968 | /* this is only to be called from in mtx_set_order */ |
969 | static void trim_incidence(mtx_matrix_t mtx, int32 order) |
970 | { |
971 | int32 ndx; |
972 | int32 *toorg; |
973 | |
974 | last_value_matrix = mtx; |
975 | toorg = mtx->perm.col.cur_to_org; |
976 | for( ndx = order ; ndx < mtx->order ; ++ndx ) |
977 | blast_col(mtx,toorg[ndx],mtx_ALL_ROWS); |
978 | for( ndx = ZERO ; ndx < order ; ++ndx ) |
979 | clean_row(mtx,toorg[ndx]); |
980 | |
981 | toorg = mtx->perm.row.cur_to_org; |
982 | for( ndx = order ; ndx < mtx->order ; ++ndx ) |
983 | blast_row(mtx,toorg[ndx],mtx_ALL_COLS); |
984 | for( ndx = ZERO ; ndx < order ; ++ndx ) |
985 | clean_col(mtx,toorg[ndx]); |
986 | } |
987 | |
988 | /* this is only to be called from in mtx_set_order */ |
989 | static void enlarge_nzheaders(mtx_matrix_t mtx, int32 order) |
990 | { |
991 | struct element_t **newhdr; |
992 | |
993 | /* realloc does not initialize, so calloc is the best we can do here */ |
994 | |
995 | newhdr = calloc_nz_header(order); |
996 | copy_nz_header(newhdr,mtx->hdr.row,mtx->capacity); |
997 | free_nz_header(mtx->hdr.row); |
998 | mtx->hdr.row = newhdr; |
999 | |
1000 | newhdr = calloc_nz_header(order); |
1001 | copy_nz_header(newhdr,mtx->hdr.col,mtx->capacity); |
1002 | free_nz_header(mtx->hdr.col); |
1003 | mtx->hdr.col = newhdr; |
1004 | } |
1005 | |
1006 | void mtx_set_order( mtx_matrix_t mtx, int32 order) |
1007 | /** |
1008 | *** This function will preserve the fact that all of the arrays are |
1009 | *** "correct" out to the capacity of the arrays, not just out to the order |
1010 | *** of the matrix. In other words, extensions to orders which still do |
1011 | *** not exceed the capacity of the arrays are trivial. |
1012 | **/ |
1013 | { |
1014 | int32 i; |
1015 | if(!mtx_check_matrix(mtx)) return; |
1016 | if (ISSLAVE(mtx)) { |
1017 | mtx_set_order(mtx->master,order); |
1018 | return; |
1019 | } |
1020 | /* we now have a master matrix */ |
1021 | if( order < mtx->order ) { /* Truncate */ |
1022 | trim_incidence(mtx,order); /* clean master */ |
1023 | for (i = 0; i < mtx->nslaves; i++) { |
1024 | trim_incidence(mtx->slaves[i],order); |
1025 | } |
1026 | } |
1027 | if (mtx->perm.transpose == mtx_NONE) { |
1028 | mtx->perm.transpose = 0; |
1029 | } |
1030 | |
1031 | if( order > mtx->capacity ) { |
1032 | int32 *newperm; |
1033 | int32 ndx; |
1034 | |
1035 | enlarge_nzheaders(mtx,order); |
1036 | for (i = 0; i < mtx->nslaves; i++) { |
1037 | enlarge_nzheaders(mtx->slaves[i],order); |
1038 | } |
1039 | |
1040 | /* realloc not in order here. Happens only on the master. */ |
1041 | newperm = mtx_alloc_perm(order); |
1042 | mtx_copy_perm(newperm,mtx->perm.row.org_to_cur,mtx->capacity); |
1043 | for( ndx=mtx->capacity ; ndx < order ; ++ndx ) |
1044 | newperm[ndx] = ndx; |
1045 | mtx_free_perm(mtx->perm.row.org_to_cur); |
1046 | mtx->perm.row.org_to_cur = newperm; |
1047 | |
1048 | newperm = mtx_alloc_perm(order); |
1049 | mtx_copy_perm(newperm,mtx->perm.row.cur_to_org,mtx->capacity); |
1050 | for( ndx=mtx->capacity ; ndx < order ; ++ndx ) |
1051 | newperm[ndx] = ndx; |
1052 | mtx_free_perm(mtx->perm.row.cur_to_org); |
1053 | mtx->perm.row.cur_to_org = newperm; |
1054 | |
1055 | newperm = mtx_alloc_perm(order); |
1056 | mtx_copy_perm(newperm,mtx->perm.col.org_to_cur,mtx->capacity); |
1057 | for( ndx=mtx->capacity ; ndx < order ; ++ndx ) |
1058 | newperm[ndx] = ndx; |
1059 | mtx_free_perm(mtx->perm.col.org_to_cur); |
1060 | mtx->perm.col.org_to_cur = newperm; |
1061 | |
1062 | newperm = mtx_alloc_perm(order); |
1063 | mtx_copy_perm(newperm,mtx->perm.col.cur_to_org,mtx->capacity); |
1064 | for( ndx=mtx->capacity ; ndx < order ; ++ndx ) |
1065 | newperm[ndx] = ndx; |
1066 | mtx_free_perm(mtx->perm.col.cur_to_org); |
1067 | mtx->perm.col.cur_to_org = newperm; |
1068 | |
1069 | mtx->capacity = order; |
1070 | for (i = 0; i < mtx->nslaves; i++) { |
1071 | mtx->slaves[i]->capacity = order; |
1072 | /* point slaves at master perm again in case anything has changed */ |
1073 | mtx->slaves[i]->perm.row.org_to_cur = mtx->perm.row.org_to_cur; |
1074 | mtx->slaves[i]->perm.row.cur_to_org = mtx->perm.row.cur_to_org; |
1075 | mtx->slaves[i]->perm.col.org_to_cur = mtx->perm.col.org_to_cur; |
1076 | mtx->slaves[i]->perm.col.cur_to_org = mtx->perm.col.cur_to_org; |
1077 | } |
1078 | } |
1079 | mtx->order = order; |
1080 | for (i = 0; i < mtx->nslaves; i++) { |
1081 | mtx->slaves[i]->order = order; |
1082 | } |
1083 | } |
1084 | |
1085 | void mtx_clear_coord(mtx_matrix_t mtx, int32 row, int32 col) |
1086 | { |
1087 | struct element_t **rlink, **clink; |
1088 | int32 org_row, org_col; |
1089 | |
1090 | #if MTX_DEBUG |
1091 | if( !mtx_check_matrix(mtx) ) return; /*ben*/ |
1092 | #endif |
1093 | |
1094 | last_value_matrix = mtx; |
1095 | org_row = mtx->perm.row.cur_to_org[row]; |
1096 | org_col = mtx->perm.col.cur_to_org[col]; |
1097 | rlink = &(mtx->hdr.row[org_row]); |
1098 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1099 | rlink = &((*rlink)->next.col); |
1100 | if( ISNULL(*rlink) ) return; |
1101 | clink = &(mtx->hdr.col[org_col]); |
1102 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1103 | clink = &((*clink)->next.row); |
1104 | delete_from_row(rlink); |
1105 | delete_from_col(clink); |
1106 | } |
1107 | |
1108 | void mtx_clear_row( mtx_matrix_t mtx, int32 row, mtx_range_t *rng) |
1109 | { |
1110 | struct element_t **rlink, **clink; |
1111 | int32 org_row; |
1112 | |
1113 | #if MTX_DEBUG |
1114 | if( !mtx_check_matrix(mtx) ) return; /*ben*/ |
1115 | #endif |
1116 | |
1117 | last_value_matrix = mtx; |
1118 | org_row = mtx->perm.row.cur_to_org[row]; |
1119 | rlink = &(mtx->hdr.row[org_row]); |
1120 | if( rng == mtx_ALL_COLS ) { |
1121 | while( NOTNULL(*rlink) ) { |
1122 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1123 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1124 | clink = &((*clink)->next.row); |
1125 | delete_from_row(rlink); |
1126 | delete_from_col(clink); |
1127 | } |
1128 | } else if( rng->high >= rng->low ) { |
1129 | while( NOTNULL(*rlink) ) { |
1130 | if( in_range(rng,mtx->perm.col.org_to_cur[(*rlink)->col]) ) { |
1131 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1132 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1133 | clink = &((*clink)->next.row); |
1134 | delete_from_row(rlink); |
1135 | delete_from_col(clink); |
1136 | } else |
1137 | rlink = &((*rlink)->next.col); |
1138 | } |
1139 | } |
1140 | } |
1141 | |
1142 | void mtx_clear_col( mtx_matrix_t mtx, int32 col, mtx_range_t *rng) |
1143 | { |
1144 | struct element_t **clink, **rlink; |
1145 | int32 org_col; |
1146 | |
1147 | #if MTX_DEBUG |
1148 | if( !mtx_check_matrix(mtx) ) return; |
1149 | #endif |
1150 | |
1151 | last_value_matrix = mtx; |
1152 | org_col = mtx->perm.col.cur_to_org[col]; |
1153 | clink = &(mtx->hdr.col[org_col]); |
1154 | if( rng == mtx_ALL_ROWS ) { |
1155 | while( NOTNULL(*clink) ) { |
1156 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1157 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1158 | rlink = &((*rlink)->next.col); |
1159 | delete_from_col(clink); |
1160 | delete_from_row(rlink); |
1161 | } |
1162 | } else if( rng->high >= rng->low ) { |
1163 | while( NOTNULL(*clink) ) { |
1164 | if( in_range(rng,mtx->perm.row.org_to_cur[(*clink)->row]) ) { |
1165 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1166 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1167 | rlink = &((*rlink)->next.col); |
1168 | delete_from_col(clink); |
1169 | delete_from_row(rlink); |
1170 | } else |
1171 | clink = &((*clink)->next.row); |
1172 | } |
1173 | } |
1174 | } |
1175 | |
1176 | void mtx_clear_region(mtx_matrix_t mtx, mtx_region_t *region) |
1177 | { |
1178 | #if MTX_DEBUG |
1179 | if(!mtx_check_matrix(mtx)) return; |
1180 | #endif |
1181 | |
1182 | if( region == mtx_ENTIRE_MATRIX && !ISSLAVE(mtx)) { |
1183 | /* damn the torpedos, wipe that sucker and slaves fast */ |
1184 | nuke_fishnet(mtx); |
1185 | } else { |
1186 | int32 ndx; |
1187 | int32 *toorg; |
1188 | mtx_region_t reg; |
1189 | |
1190 | if (region == mtx_ENTIRE_MATRIX) { |
1191 | reg.row.low = reg.col.low = 0; |
1192 | reg.row.high = reg.col.high = mtx->order-1; |
1193 | } else { |
1194 | reg = *region; |
1195 | } |
1196 | last_value_matrix = mtx; |
1197 | toorg = mtx->perm.row.cur_to_org; |
1198 | for( ndx = reg.row.low ; ndx <= reg.row.high ; ++ndx ) |
1199 | blast_row(mtx,toorg[ndx],&(reg.col)); |
1200 | toorg = mtx->perm.col.cur_to_org; |
1201 | for( ndx = reg.col.low ; ndx <= reg.col.high ; ++ndx ) |
1202 | clean_col(mtx,toorg[ndx]); |
1203 | } |
1204 | } |
1205 | |
1206 | |
1207 | void mtx_reset_perm(mtx_matrix_t mtx) |
1208 | { |
1209 | int32 ndx; |
1210 | |
1211 | #if MTX_DEBUG |
1212 | if(!mtx_check_matrix(mtx)) return; |
1213 | #endif |
1214 | if (ISSLAVE(mtx)) { |
1215 | mtx_reset_perm(mtx->master); |
1216 | return; |
1217 | } |
1218 | |
1219 | for( ndx=ZERO ; ndx < mtx->capacity ; ++ndx ) { |
1220 | mtx->perm.row.org_to_cur[ndx] = ndx; |
1221 | mtx->perm.row.cur_to_org[ndx] = ndx; |
1222 | mtx->perm.col.org_to_cur[ndx] = ndx; |
1223 | mtx->perm.col.cur_to_org[ndx] = ndx; |
1224 | } |
1225 | mtx->perm.row.parity = mtx->perm.col.parity = EVEN; |
1226 | |
1227 | mtx->data->symbolic_rank = -1; |
1228 | mtx->data->nblocks = -1; |
1229 | if( NOTNULL(mtx->data->block) ) { |
1230 | ascfree(mtx->data->block); |
1231 | mtx->data->block = NULL; |
1232 | } |
1233 | } |
1234 | |
1235 | void mtx_clear(mtx_matrix_t mtx) |
1236 | { |
1237 | if (ISSLAVE(mtx)) { |
1238 | mtx_clear_region(mtx->master,mtx_ENTIRE_MATRIX); |
1239 | } |
1240 | mtx_clear_region(mtx,mtx_ENTIRE_MATRIX); |
1241 | mtx_reset_perm(mtx); |
1242 | } |
1243 | |
1244 | real64 mtx_value(mtx_matrix_t mtx, mtx_coord_t *coord) |
1245 | { |
1246 | struct element_t *elt; |
1247 | register int32 org_row,org_col; |
1248 | |
1249 | #if MTX_DEBUG |
1250 | if(!mtx_check_matrix(mtx)) return D_ZERO; |
1251 | #endif |
1252 | org_row = mtx->perm.row.cur_to_org[coord->row]; |
1253 | org_col = mtx->perm.col.cur_to_org[coord->col]; |
1254 | elt = mtx_find_element(mtx,org_row,org_col); |
1255 | mtx->last_value = elt; |
1256 | return( ISNULL(elt) ? D_ZERO : elt->value ); |
1257 | } |
1258 | |
1259 | void mtx_set_value( mtx_matrix_t mtx, mtx_coord_t *coord, real64 value) |
1260 | { |
1261 | register int32 org_row,org_col; |
1262 | |
1263 | #if MTX_DEBUG |
1264 | if(!mtx_check_matrix(mtx)) return; |
1265 | #endif |
1266 | org_row = mtx->perm.row.cur_to_org[coord->row]; |
1267 | org_col = mtx->perm.col.cur_to_org[coord->col]; |
1268 | if ( NOTNULL(mtx->last_value) && |
1269 | mtx->last_value->row==org_row && |
1270 | mtx->last_value->col==org_col ) { |
1271 | mtx->last_value->value = value; |
1272 | } else { |
1273 | struct element_t *elt; |
1274 | if( ISNULL(elt = mtx_find_element(mtx,org_row,org_col)) ) { |
1275 | if (value != D_ZERO ) { |
1276 | elt = mtx_create_element_value(mtx,org_row,org_col,value); |
1277 | } |
1278 | } else { |
1279 | elt->value = value; |
1280 | } |
1281 | } |
1282 | } |
1283 | |
1284 | void mtx_fill_value(mtx_matrix_t mtx, mtx_coord_t *coord, real64 value) |
1285 | { |
1286 | register int32 org_row,org_col; |
1287 | |
1288 | #if MTX_DEBUG |
1289 | if(!mtx_check_matrix(mtx)) return; |
1290 | #endif |
1291 | org_row = mtx->perm.row.cur_to_org[coord->row]; |
1292 | org_col = mtx->perm.col.cur_to_org[coord->col]; |
1293 | mtx_create_element_value(mtx,org_row,org_col,value); |
1294 | } |
1295 | |
1296 | void mtx_fill_org_value(mtx_matrix_t mtx, const mtx_coord_t *coord, |
1297 | real64 value) |
1298 | { |
1299 | #if MTX_DEBUG |
1300 | if(!mtx_check_matrix(mtx)) return; |
1301 | #endif |
1302 | mtx_create_element_value(mtx,coord->row,coord->col,value); |
1303 | } |
1304 | |
1305 | /* sparse matrix assembly of potentially duplicate fills */ |
1306 | /* takes a matrix, assumed to have redundant and otherwise insane incidences |
1307 | created by 'misusing mtx_fill_value' and sums all like entries, eliminating |
1308 | the duplicates and the zeroes. Returns -# of elements removed, or 1 if bad. |
1309 | returns 1 if fails for some reason, 0 otherwise. |
1310 | Could stand to have the error messages it emits improved. |
1311 | Could stand to take a rowrange or a rowlist, a colrange or a collist,droptol. |
1312 | algorithm: O(3tau) |
1313 | establish a vector mv of columns needing cleaning markers |
1314 | set lowmark=-1, highmark= -2 |
1315 | for rows |
1316 | if row empty continue |
1317 | establish a vector ev of elt pointers set null |
1318 | prevlink = &header 1st elt pointer |
1319 | for each elt in row |
1320 | if elt is not marked in ev |
1321 | mark ev with elt, update prevlink |
1322 | else |
1323 | add value to elt pointed at by ev |
1324 | mark col in mv as needing cleaning, updating lowmark,highmark |
1325 | delete elt |
1326 | endif |
1327 | endfor |
1328 | for each elt in row |
1329 | unmark ev with elt |
1330 | if elt is 0, delete and mark col in mv. |
1331 | endfor |
1332 | endfor |
1333 | for i = lowmark to highmark |
1334 | clean col |
1335 | endfor |
1336 | My god, the tricks we play on a linked list. |
1337 | */ |
1338 | int32 mtx_assemble(mtx_matrix_t mtx) |
1339 | { |
1340 | char *real_mv=NULL, *mv; |
1341 | struct element_t **real_ev=NULL, **ev, **prevlink, *elt; |
1342 | int32 orgrow, orgcol, lowmark,highmark,dup; |
1343 | /* note elt and prevlink could be combined, but the code is |
1344 | unreadable utterly if you do so. */ |
1345 | |
1346 | if (ISNULL(mtx)) return 1; |
1347 | if (mtx->order <1) return 0; |
1348 | real_mv = mtx_null_mark(mtx->order+1); |
1349 | mv = real_mv+1; |
1350 | if (ISNULL(real_mv)) return 1; |
1351 | real_ev = mtx_null_vector(mtx->order+1); |
1352 | ev = real_ev+1; |
1353 | if (ISNULL(real_ev)) { |
1354 | mtx_null_mark_release(); |
1355 | return 1; |
1356 | } |
1357 | /* we have allocated arrays which include a -1 element to buy |
1358 | ourselves an awful convenient lot of safety. */ |
1359 | lowmark=-1; |
1360 | highmark=-2; |
1361 | dup = 0; |
1362 | last_value_matrix = mtx; |
1363 | for (orgrow=0; orgrow < mtx->order; orgrow++) { |
1364 | elt = mtx->hdr.row[orgrow]; |
1365 | prevlink = &(mtx->hdr.row[orgrow]); |
1366 | while (NOTNULL(elt)) { |
1367 | if (ISNULL(ev[elt->col])) { |
1368 | ev[elt->col] = elt; /* mark first elt found for this col */ |
1369 | prevlink= &(elt->next.col); /* collect pointer to where we go next */ |
1370 | elt = elt->next.col; /* go there */ |
1371 | } else { |
1372 | /* elt is duplicate and must die */ |
1373 | dup++; |
1374 | ev[elt->col]->value += elt->value; /* accumulate value */ |
1375 | /* update lowmark, highmark . this is a debatable implementation. */ |
1376 | /* for small mods on large matrices, this is sane */ |
1377 | if (lowmark > -1) { /* not first mark */ |
1378 | if (elt->col < lowmark && elt->col > -1) { |
1379 | lowmark = elt->col; |
1380 | } |
1381 | if (elt->col > highmark && elt->col < mtx->order) { |
1382 | highmark = elt->col; |
1383 | } |
1384 | } else { /* very first mark */ |
1385 | if (elt->col > -1 && elt->col < mtx->order) { |
1386 | lowmark = highmark = elt->col; |
1387 | } |
1388 | } |
1389 | mv[elt->col] = 1; /* mark column as to be cleaned */ |
1390 | delete_from_row(prevlink); |
1391 | elt = *prevlink; |
1392 | } |
1393 | } |
1394 | elt = mtx->hdr.row[orgrow]; |
1395 | prevlink = &(mtx->hdr.row[orgrow]); |
1396 | while (NOTNULL(elt)) { /* renull the accumulator and trash 0s */ |
1397 | ev[elt->col] = NULL; /* regardless, reset accum */ |
1398 | if (elt->value != D_ZERO) { |
1399 | prevlink= &(elt->next.col); /* collect pointer to where we go next */ |
1400 | elt = elt->next.col; /* go there */ |
1401 | } else { |
1402 | /* this is still a debatable implementation. */ |
1403 | if (lowmark > -1) { /* not first mark */ |
1404 | if (elt->col < lowmark && elt->col > -1) { |
1405 | lowmark = elt->col; |
1406 | } |
1407 | if (elt->col > highmark && elt->col < mtx->order) { |
1408 | highmark = elt->col; |
1409 | } |
1410 | } else { /* very first mark */ |
1411 | if (elt->col > -1 && elt->col < mtx->order) { |
1412 | lowmark = highmark = elt->col; |
1413 | } |
1414 | } |
1415 | mv[elt->col] = 1; /* mark column as to be cleaned */ |
1416 | delete_from_row(prevlink); |
1417 | elt = *prevlink; |
1418 | } |
1419 | } |
1420 | } |
1421 | for (orgcol = lowmark; orgcol <= highmark; orgcol++) { |
1422 | if (mv[orgcol]) { |
1423 | clean_col(mtx,orgcol); /* scrap dups and 0s */ |
1424 | } |
1425 | } |
1426 | mtx_null_mark_release(); |
1427 | mtx_null_vector_release(); |
1428 | return -dup; |
1429 | } |
1430 | |
1431 | void mtx_del_zr_in_row(mtx_matrix_t mtx, int32 row) |
1432 | { |
1433 | register int32 org; |
1434 | struct element_t **rlink, **clink; |
1435 | |
1436 | #if MTX_DEBUG |
1437 | if(!mtx_check_matrix(mtx)) return; |
1438 | #endif |
1439 | org = mtx->perm.row.cur_to_org[row]; |
1440 | rlink = &(mtx->hdr.row[org]); |
1441 | |
1442 | last_value_matrix = mtx; |
1443 | while( NOTNULL(*rlink) ) |
1444 | if( (*rlink)->value == D_ZERO ) { |
1445 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1446 | for( ; NOTNULL(*clink) && (*clink)->row != org ; ) |
1447 | clink = &((*clink)->next.row); |
1448 | delete_from_row(rlink); |
1449 | delete_from_col(clink); |
1450 | } else |
1451 | rlink = &((*rlink)->next.col); |
1452 | } |
1453 | |
1454 | void mtx_del_zr_in_col(mtx_matrix_t mtx, int32 col) |
1455 | { |
1456 | register int32 org; |
1457 | struct element_t **clink, **rlink; |
1458 | |
1459 | #if MTX_DEBUG |
1460 | if(!mtx_check_matrix(mtx)) return; |
1461 | #endif |
1462 | org = mtx->perm.col.cur_to_org[col]; |
1463 | clink = &(mtx->hdr.col[org]); |
1464 | |
1465 | last_value_matrix = mtx; |
1466 | while( NOTNULL(*clink) ) |
1467 | if( (*clink)->value == D_ZERO ) { |
1468 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1469 | for( ; NOTNULL(*rlink) && (*rlink)->col != org ; ) |
1470 | rlink = &((*rlink)->next.col); |
1471 | delete_from_col(clink); |
1472 | delete_from_row(rlink); |
1473 | } else |
1474 | clink = &((*clink)->next.row); |
1475 | } |
1476 | |
1477 | void mtx_del_zr_in_rowrange(mtx_matrix_t mtx, mtx_range_t *rng) |
1478 | { |
1479 | register int32 org,row,rowhi, *toorg; |
1480 | struct element_t **rlink, **clink; |
1481 | |
1482 | #if MTX_DEBUG |
1483 | if(!mtx_check_matrix(mtx)) return; |
1484 | #endif |
1485 | rowhi=rng->high; |
1486 | toorg= mtx->perm.row.cur_to_org; |
1487 | last_value_matrix = mtx; |
1488 | for (row=rng->low; row <=rowhi; row++) { |
1489 | org = toorg[row]; |
1490 | rlink = &(mtx->hdr.row[org]); |
1491 | |
1492 | while( NOTNULL(*rlink) ) { |
1493 | if( (*rlink)->value == D_ZERO ) { |
1494 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1495 | for( ; NOTNULL(*clink) && (*clink)->row != org ; ) |
1496 | clink = &((*clink)->next.row); |
1497 | delete_from_row(rlink); |
1498 | delete_from_col(clink); |
1499 | } else { |
1500 | rlink = &((*rlink)->next.col); |
1501 | } |
1502 | } |
1503 | |
1504 | } |
1505 | } |
1506 | |
1507 | void mtx_del_zr_in_colrange(mtx_matrix_t mtx, mtx_range_t *rng) |
1508 | { |
1509 | register int32 org,col,colhi, *toorg; |
1510 | struct element_t **clink, **rlink; |
1511 | |
1512 | #if MTX_DEBUG |
1513 | if(!mtx_check_matrix(mtx)) return; |
1514 | #endif |
1515 | colhi=rng->high; |
1516 | toorg= mtx->perm.col.cur_to_org; |
1517 | last_value_matrix = mtx; |
1518 | for (col=rng->low; col <=colhi; col++) { |
1519 | org = toorg[col]; |
1520 | clink = &(mtx->hdr.col[org]); |
1521 | |
1522 | while( NOTNULL(*clink) ) { |
1523 | if( (*clink)->value == D_ZERO ) { |
1524 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1525 | for( ; NOTNULL(*rlink) && (*rlink)->col != org ; ) |
1526 | rlink = &((*rlink)->next.col); |
1527 | delete_from_col(clink); |
1528 | delete_from_row(rlink); |
1529 | } else { |
1530 | clink = &((*clink)->next.row); |
1531 | } |
1532 | } |
1533 | |
1534 | } |
1535 | } |
1536 | |
1537 | void mtx_steal_org_row_vec(mtx_matrix_t mtx, int32 row, |
1538 | real64 *vec, mtx_range_t *rng) |
1539 | { |
1540 | struct element_t **rlink, **clink; |
1541 | int32 org_row; |
1542 | |
1543 | #if MTX_DEBUG |
1544 | if( !mtx_check_matrix(mtx) ) return; |
1545 | #endif |
1546 | |
1547 | last_value_matrix = mtx; |
1548 | org_row = mtx->perm.row.cur_to_org[row]; |
1549 | rlink = &(mtx->hdr.row[org_row]); |
1550 | if( rng == mtx_ALL_COLS ) { |
1551 | while( NOTNULL(*rlink) ) { |
1552 | vec[(*rlink)->col] = (*rlink)->value; |
1553 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1554 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1555 | clink = &((*clink)->next.row); |
1556 | delete_from_row(rlink); |
1557 | delete_from_col(clink); |
1558 | } |
1559 | } else if( rng->high >= rng->low ) { |
1560 | int32 *tocur; |
1561 | tocur = mtx->perm.col.org_to_cur; |
1562 | while( NOTNULL(*rlink) ) { |
1563 | if( in_range(rng,tocur[(*rlink)->col]) ) { |
1564 | vec[(*rlink)->col] = (*rlink)->value; |
1565 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1566 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1567 | clink = &((*clink)->next.row); |
1568 | delete_from_row(rlink); |
1569 | delete_from_col(clink); |
1570 | } else { |
1571 | rlink = &((*rlink)->next.col); |
1572 | } |
1573 | } |
1574 | } |
1575 | } |
1576 | |
1577 | void mtx_steal_org_col_vec(mtx_matrix_t mtx, int32 col, |
1578 | real64 *vec, mtx_range_t *rng) |
1579 | { |
1580 | struct element_t **clink, **rlink; |
1581 | int32 org_col; |
1582 | |
1583 | #if MTX_DEBUG |
1584 | if( !mtx_check_matrix(mtx) ) return; |
1585 | #endif |
1586 | |
1587 | last_value_matrix = mtx; |
1588 | org_col = mtx->perm.col.cur_to_org[col]; |
1589 | clink = &(mtx->hdr.col[org_col]); |
1590 | if( rng == mtx_ALL_ROWS ) { |
1591 | while( NOTNULL(*clink) ) { |
1592 | vec[(*clink)->row] = (*clink)->value; |
1593 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1594 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1595 | rlink = &((*rlink)->next.col); |
1596 | delete_from_col(clink); |
1597 | delete_from_row(rlink); |
1598 | } |
1599 | } else if( rng->high >= rng->low ) { |
1600 | int32 *tocur; |
1601 | tocur = mtx->perm.row.org_to_cur; |
1602 | while( NOTNULL(*clink) ) { |
1603 | if( in_range(rng,tocur[(*clink)->row]) ) { |
1604 | vec[(*clink)->row] = (*clink)->value; |
1605 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1606 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1607 | rlink = &((*rlink)->next.col); |
1608 | delete_from_col(clink); |
1609 | delete_from_row(rlink); |
1610 | } else { |
1611 | clink = &((*clink)->next.row); |
1612 | } |
1613 | } |
1614 | } |
1615 | } |
1616 | |
1617 | void mtx_steal_cur_row_vec(mtx_matrix_t mtx, int32 row, |
1618 | real64 *vec, mtx_range_t *rng) |
1619 | { |
1620 | struct element_t **rlink, **clink; |
1621 | int32 org_row, *tocur; |
1622 | |
1623 | #if MTX_DEBUG |
1624 | if( !mtx_check_matrix(mtx) ) return; |
1625 | #endif |
1626 | |
1627 | tocur = mtx->perm.col.org_to_cur; |
1628 | last_value_matrix = mtx; |
1629 | org_row = mtx->perm.row.cur_to_org[row]; |
1630 | rlink = &(mtx->hdr.row[org_row]); |
1631 | if( rng == mtx_ALL_COLS ) { |
1632 | while( NOTNULL(*rlink) ) { |
1633 | vec[tocur[(*rlink)->col]] = (*rlink)->value; |
1634 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1635 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1636 | clink = &((*clink)->next.row); |
1637 | delete_from_row(rlink); |
1638 | delete_from_col(clink); |
1639 | } |
1640 | } else if( rng->high >= rng->low ) { |
1641 | while( NOTNULL(*rlink) ) { |
1642 | if( in_range(rng,tocur[(*rlink)->col]) ) { |
1643 | vec[tocur[(*rlink)->col]] = (*rlink)->value; |
1644 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1645 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1646 | clink = &((*clink)->next.row); |
1647 | delete_from_row(rlink); |
1648 | delete_from_col(clink); |
1649 | } else { |
1650 | rlink = &((*rlink)->next.col); |
1651 | } |
1652 | } |
1653 | } |
1654 | } |
1655 | |
1656 | void mtx_steal_cur_col_vec(mtx_matrix_t mtx, int32 col, |
1657 | real64 *vec, mtx_range_t *rng) |
1658 | { |
1659 | struct element_t **clink, **rlink; |
1660 | int32 org_col, *tocur; |
1661 | |
1662 | #if MTX_DEBUG |
1663 | if( !mtx_check_matrix(mtx) ) return; |
1664 | #endif |
1665 | |
1666 | tocur = mtx->perm.row.org_to_cur; |
1667 | last_value_matrix = mtx; |
1668 | org_col = mtx->perm.col.cur_to_org[col]; |
1669 | clink = &(mtx->hdr.col[org_col]); |
1670 | if( rng == mtx_ALL_ROWS ) { |
1671 | while( NOTNULL(*clink) ) { |
1672 | vec[tocur[(*clink)->row]] = (*clink)->value; |
1673 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1674 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1675 | rlink = &((*rlink)->next.col); |
1676 | delete_from_col(clink); |
1677 | delete_from_row(rlink); |
1678 | } |
1679 | } else if( rng->high >= rng->low ) { |
1680 | while( NOTNULL(*clink) ) { |
1681 | if( in_range(rng,tocur[(*clink)->row]) ) { |
1682 | vec[tocur[(*clink)->row]] = (*clink)->value; |
1683 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1684 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1685 | rlink = &((*rlink)->next.col); |
1686 | delete_from_col(clink); |
1687 | delete_from_row(rlink); |
1688 | } else { |
1689 | clink = &((*clink)->next.row); |
1690 | } |
1691 | } |
1692 | } |
1693 | } |
1694 | |
1695 | #ifdef DELETE_THIS_UNUSED_FUNCTION |
1696 | /* Little function to enlarge the capacity of a sparse to the len given. |
1697 | * New memory allocated, if any, is not initialized in any way. |
1698 | * The pointer given is not itself changed, only its data. |
1699 | * This function in no way shrinks the data in a sparse. |
1700 | * (exception: buggy realloc implementations may shrink it, ack!) |
1701 | * Data/idata values up to the old capacity (ret->cap) are preserved. |
1702 | */ |
1703 | static int enlarge_sparse(mtx_sparse_t *ret, int32 len) |
1704 | { |
1705 | int32 *inew; |
1706 | real64 *dnew; |
1707 | if (ISNULL(ret)) { |
1708 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) enlarge_sparse called with NULL.\n"); |
1709 | return 1; |
1710 | } |
1711 | if (len <= ret->cap || len <1) return 0; /* already big enough */ |
1712 | |
1713 | if (ret->idata == NULL) { |
1714 | inew = (int32 *)ascmalloc(sizeof(int32)*len); |
1715 | } else { |
1716 | inew = (int32 *)ascrealloc(ret->idata,sizeof(int32)*len); |
1717 | } |
1718 | if (ISNULL(inew)) { |
1719 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) enlarge_sparse\n"); |
1720 | FPRINTF(g_mtxerr," Insufficient memory.\n"); |
1721 | return 1; |
1722 | } |
1723 | ret->idata = inew; /* dnew can still fail without losing inew memory. */ |
1724 | |
1725 | if (ret->data == NULL) { |
1726 | dnew = (real64 *)ascmalloc(sizeof(real64)*len); |
1727 | } else { |
1728 | dnew = (real64 *)ascrealloc(ret->data,sizeof(real64)*len); |
1729 | } |
1730 | if (ISNULL(dnew)) { |
1731 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) enlarge_sparse\n"); |
1732 | FPRINTF(g_mtxerr," Insufficient memory.\n"); |
1733 | return 1; |
1734 | } |
1735 | ret->data = dnew; /* we succeeded */ |
1736 | ret->cap = len; |
1737 | return 0; |
1738 | } |
1739 | #endif /* DELETE_THIS_UNUSED_FUNCTION */ |
1740 | |
1741 | |
1742 | /* going to try to make steal also handle sparse creation ...*/ |
1743 | /* don't you dare, whoever you are! */ |
1744 | boolean mtx_steal_org_row_sparse(mtx_matrix_t mtx, int32 row, |
1745 | mtx_sparse_t *sp, mtx_range_t *rng) |
1746 | { |
1747 | struct element_t **rlink, **clink; |
1748 | int32 org_row; |
1749 | int32 len,k; |
1750 | |
1751 | #if MTX_DEBUG |
1752 | if( !mtx_check_matrix(mtx) ) return TRUE; |
1753 | #endif |
1754 | if (sp == mtx_CREATE_SPARSE) { |
1755 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_org_row_sparse called with\n"); |
1756 | FPRINTF(g_mtxerr," mtx_CREATE_SPARSE. Not supported.\n"); |
1757 | return TRUE; |
1758 | } |
1759 | if (rng == mtx_ALL_COLS) { |
1760 | len = mtx->order; |
1761 | } else { |
1762 | len = rng->high - rng->low +1; |
1763 | } |
1764 | if (sp->cap < len) { |
1765 | sp->len = 0; |
1766 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_org_row_sparse called with\n"); |
1767 | FPRINTF(g_mtxerr," sparse of insufficient capacity.\n"); |
1768 | return TRUE; |
1769 | } |
1770 | |
1771 | last_value_matrix = mtx; |
1772 | org_row = mtx->perm.row.cur_to_org[row]; |
1773 | rlink = &(mtx->hdr.row[org_row]); |
1774 | k = 0; |
1775 | |
1776 | if( rng == mtx_ALL_COLS ) { |
1777 | while( NOTNULL(*rlink) ) { |
1778 | sp->idata[k] = (*rlink)->col; |
1779 | sp->data[k] = (*rlink)->value; |
1780 | k++; |
1781 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1782 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1783 | clink = &((*clink)->next.row); |
1784 | delete_from_row(rlink); |
1785 | delete_from_col(clink); |
1786 | } |
1787 | } else if( rng->high >= rng->low ) { |
1788 | int32 *tocur; |
1789 | tocur = mtx->perm.col.org_to_cur; |
1790 | while( NOTNULL(*rlink) ) { |
1791 | if( in_range(rng,tocur[(*rlink)->col]) ) { |
1792 | sp->idata[k] = (*rlink)->col; |
1793 | sp->data[k] = (*rlink)->value; |
1794 | k++; |
1795 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1796 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1797 | clink = &((*clink)->next.row); |
1798 | delete_from_row(rlink); |
1799 | delete_from_col(clink); |
1800 | } else { |
1801 | rlink = &((*rlink)->next.col); |
1802 | } |
1803 | } |
1804 | } |
1805 | sp->len = k; |
1806 | return FALSE; |
1807 | } |
1808 | |
1809 | boolean mtx_steal_org_col_sparse(mtx_matrix_t mtx, int32 col, |
1810 | mtx_sparse_t *sp, mtx_range_t *rng) |
1811 | { |
1812 | struct element_t **clink, **rlink; |
1813 | int32 org_col; |
1814 | int32 len,k; |
1815 | |
1816 | #if MTX_DEBUG |
1817 | if( !mtx_check_matrix(mtx) ) return TRUE; |
1818 | #endif |
1819 | if (sp == mtx_CREATE_SPARSE) { |
1820 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_org_col_sparse called with\n"); |
1821 | FPRINTF(g_mtxerr," mtx_CREATE_SPARSE. Not supported.\n"); |
1822 | return TRUE; |
1823 | } |
1824 | if (rng == mtx_ALL_ROWS) { |
1825 | len = mtx->order; |
1826 | } else { |
1827 | len = rng->high - rng->low +1; |
1828 | } |
1829 | if (sp->cap < len) { |
1830 | sp->len = 0; |
1831 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_org_col_sparse called with\n"); |
1832 | FPRINTF(g_mtxerr," sparse of insufficient capacity.\n"); |
1833 | return TRUE; |
1834 | } |
1835 | |
1836 | last_value_matrix = mtx; |
1837 | org_col = mtx->perm.col.cur_to_org[col]; |
1838 | clink = &(mtx->hdr.col[org_col]); |
1839 | k = 0; |
1840 | |
1841 | if( rng == mtx_ALL_ROWS ) { |
1842 | while( NOTNULL(*clink) ) { |
1843 | sp->idata[k] = (*clink)->row; |
1844 | sp->data[k] = (*clink)->value; |
1845 | k++; |
1846 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1847 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1848 | rlink = &((*rlink)->next.col); |
1849 | delete_from_col(clink); |
1850 | delete_from_row(rlink); |
1851 | } |
1852 | } else if( rng->high >= rng->low ) { |
1853 | int32 *tocur; |
1854 | tocur = mtx->perm.row.org_to_cur; |
1855 | while( NOTNULL(*clink) ) { |
1856 | if( in_range(rng,tocur[(*clink)->row]) ) { |
1857 | sp->idata[k] = (*clink)->row; |
1858 | sp->data[k] = (*clink)->value; |
1859 | k++; |
1860 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1861 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1862 | rlink = &((*rlink)->next.col); |
1863 | delete_from_col(clink); |
1864 | delete_from_row(rlink); |
1865 | } else { |
1866 | clink = &((*clink)->next.row); |
1867 | } |
1868 | } |
1869 | } |
1870 | sp->len = k; |
1871 | return FALSE; |
1872 | } |
1873 | |
1874 | boolean mtx_steal_cur_row_sparse(mtx_matrix_t mtx, int32 row, |
1875 | mtx_sparse_t *sp, mtx_range_t *rng) |
1876 | { |
1877 | struct element_t **rlink, **clink; |
1878 | int32 org_row, *tocur; |
1879 | int32 len,k; |
1880 | |
1881 | #if MTX_DEBUG |
1882 | if( !mtx_check_matrix(mtx) ) return TRUE; |
1883 | #endif |
1884 | if (sp == mtx_CREATE_SPARSE) { |
1885 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_cur_row_sparse called with\n"); |
1886 | FPRINTF(g_mtxerr," mtx_CREATE_SPARSE. Not supported.\n"); |
1887 | return TRUE; |
1888 | } |
1889 | if (rng == mtx_ALL_COLS) { |
1890 | len = mtx->order; |
1891 | } else { |
1892 | len = rng->high - rng->low +1; |
1893 | } |
1894 | if (sp->cap < len) { |
1895 | sp->len = 0; |
1896 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_cur_row_sparse called with\n"); |
1897 | FPRINTF(g_mtxerr," sparse of insufficient capacity.\n"); |
1898 | return TRUE; |
1899 | } |
1900 | |
1901 | tocur = mtx->perm.col.org_to_cur; |
1902 | last_value_matrix = mtx; |
1903 | org_row = mtx->perm.row.cur_to_org[row]; |
1904 | rlink = &(mtx->hdr.row[org_row]); |
1905 | k = 0; |
1906 | |
1907 | if( rng == mtx_ALL_COLS ) { |
1908 | while( NOTNULL(*rlink) ) { |
1909 | sp->idata[k] = tocur[(*rlink)->col]; |
1910 | sp->data[k] = (*rlink)->value; |
1911 | k++; |
1912 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1913 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1914 | clink = &((*clink)->next.row); |
1915 | delete_from_row(rlink); |
1916 | delete_from_col(clink); |
1917 | } |
1918 | } else if( rng->high >= rng->low ) { |
1919 | while( NOTNULL(*rlink) ) { |
1920 | if( in_range(rng,tocur[(*rlink)->col]) ) { |
1921 | sp->idata[k] = tocur[(*rlink)->col]; |
1922 | sp->data[k] = (*rlink)->value; |
1923 | k++; |
1924 | clink = &(mtx->hdr.col[(*rlink)->col]); |
1925 | for( ; NOTNULL(*clink) && (*clink)->row != org_row ; ) |
1926 | clink = &((*clink)->next.row); |
1927 | delete_from_row(rlink); |
1928 | delete_from_col(clink); |
1929 | } else { |
1930 | rlink = &((*rlink)->next.col); |
1931 | } |
1932 | } |
1933 | } |
1934 | sp->len = k; |
1935 | return FALSE; |
1936 | } |
1937 | |
1938 | boolean mtx_steal_cur_col_sparse(mtx_matrix_t mtx, int32 col, |
1939 | mtx_sparse_t *sp, mtx_range_t *rng) |
1940 | { |
1941 | struct element_t **clink, **rlink; |
1942 | int32 org_col, *tocur; |
1943 | int32 len,k; |
1944 | |
1945 | #if MTX_DEBUG |
1946 | if( !mtx_check_matrix(mtx) ) return TRUE; |
1947 | #endif |
1948 | if (sp == mtx_CREATE_SPARSE) { |
1949 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_cur_col_sparse called with\n"); |
1950 | FPRINTF(g_mtxerr," mtx_CREATE_SPARSE. Not supported.\n"); |
1951 | return TRUE; |
1952 | } |
1953 | if (rng == mtx_ALL_ROWS) { |
1954 | len = mtx->order; |
1955 | } else { |
1956 | len = rng->high - rng->low +1; |
1957 | } |
1958 | if (sp->cap < len) { |
1959 | sp->len = 0; |
1960 | FPRINTF(g_mtxerr,"ERROR: (mtx.c) mtx_steal_cur_col_sparse called with\n"); |
1961 | FPRINTF(g_mtxerr," sparse of insufficient capacity.\n"); |
1962 | return TRUE; |
1963 | } |
1964 | |
1965 | tocur = mtx->perm.row.org_to_cur; |
1966 | last_value_matrix = mtx; |
1967 | org_col = mtx->perm.col.cur_to_org[col]; |
1968 | clink = &(mtx->hdr.col[org_col]); |
1969 | k = 0; |
1970 | |
1971 | if( rng == mtx_ALL_ROWS ) { |
1972 | while( NOTNULL(*clink) ) { |
1973 | sp->idata[k] = tocur[(*clink)->row]; |
1974 | sp->data[k] = (*clink)->value; |
1975 | k++; |
1976 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1977 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1978 | rlink = &((*rlink)->next.col); |
1979 | delete_from_col(clink); |
1980 | delete_from_row(rlink); |
1981 | } |
1982 | } else if( rng->high >= rng->low ) { |
1983 | while( NOTNULL(*clink) ) { |
1984 | if( in_range(rng,tocur[(*clink)->row]) ) { |
1985 | sp->idata[k] = tocur[(*clink)->row]; |
1986 | sp->data[k] = (*clink)->value; |
1987 | k++; |
1988 | rlink = &(mtx->hdr.row[(*clink)->row]); |
1989 | for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; ) |
1990 | rlink = &((*rlink)->next.col); |
1991 | delete_from_col(clink); |
1992 | delete_from_row(rlink); |
1993 | } else { |
1994 | clink = &((*clink)->next.row); |
1995 | } |
1996 | } |
1997 | } |
1998 | sp->len = k; |
1999 | return FALSE; |
2000 | } |
2001 | |
2002 | void mtx_fill_org_row_vec(mtx_matrix_t mtx, int32 row, |
2003 | real64 *vec, mtx_range_t *rng) |
2004 | { |
2005 | int32 org_row,highcol, *toorg; |
2006 | register int32 org_col; |
2007 | register real64 value; |
2008 | |
2009 | #if MTX_DEBUG |
2010 | if(!mtx_check_matrix(mtx)) return; |
2011 | #endif |
2012 | |
2013 | org_row = mtx->perm.row.cur_to_org[row]; |
2014 | if( rng == mtx_ALL_COLS ) { |
2015 | highcol=mtx->order; |
2016 | for( org_col=0 ; org_col <highcol ; org_col++ ) { |
2017 | if ((value=vec[org_col]) != D_ZERO) |
2018 | mtx_create_element_value(mtx,org_row,org_col,value); |
2019 | } |
2020 | } else { |
2021 | register int32 cur_col; |
2022 | |
2023 | toorg = mtx->perm.col.cur_to_org; |
2024 | highcol= rng->high; |
2025 | for ( cur_col=rng->low; cur_col<=highcol; cur_col++) { |
2026 | if ((value=vec[(org_col=toorg[cur_col])]) != D_ZERO) |
2027 | mtx_create_element_value(mtx,org_row,org_col,value); |
2028 | } |
2029 | } |
2030 | } |
2031 | |
2032 | void mtx_fill_org_col_vec(mtx_matrix_t mtx, int32 col, |
2033 | real64 *vec, mtx_range_t *rng) |
2034 | { |
2035 | int32 org_col,highrow, *toorg; |
2036 | register int32 org_row; |
2037 | register real64 value; |
2038 | |
2039 | #if MTX_DEBUG |
2040 | if(!mtx_check_matrix(mtx)) return; |
2041 | #endif |
2042 | |
2043 | org_col = mtx->perm.col.cur_to_org[col]; |
2044 | if( rng == mtx_ALL_ROWS ) { |
2045 | highrow=mtx->order; |
2046 | for( org_row=0 ; org_row <highrow ; org_row++ ) { |
2047 | if ((value=vec[org_row]) != D_ZERO) |
2048 | mtx_create_element_value(mtx,org_row,org_col,value); |
2049 | } |
2050 | } else { |
2051 | register int32 cur_row; |
2052 | |
2053 | toorg = mtx->perm.row.cur_to_org; |
2054 | highrow= rng->high; |
2055 | for ( cur_row=rng->low; cur_row<=highrow; cur_row++) { |
2056 | if ((value=vec[(org_row=toorg[cur_row])]) != D_ZERO) |
2057 | mtx_create_element_value(mtx,org_row,org_col,value); |
2058 | } |
2059 | } |
2060 | } |
2061 | |
2062 | void mtx_fill_cur_row_vec(mtx_matrix_t mtx, int32 row, |
2063 | real64 *vec, mtx_range_t *rng) |
2064 | { |
2065 | int32 org_row,highcol,lowcol, *toorg; |
2066 | register int32 cur_col; |
2067 | register real64 value; |
2068 | |
2069 | #if MTX_DEBUG |
2070 | if(!mtx_check_matrix(mtx)) return; |
2071 | #endif |
2072 | |
2073 | org_row = mtx->perm.row.cur_to_org[row]; |
2074 | toorg=mtx->perm.col.cur_to_org; |
2075 | if( rng == mtx_ALL_COLS ) { |
2076 | highcol = mtx->order-1; |
2077 | lowcol = 0; |
2078 | } else { |
2079 | highcol= rng->high; |
2080 | lowcol=rng->low; |
2081 | } |
2082 | for( cur_col=lowcol ; cur_col <= highcol ; cur_col++ ) { |
2083 | if ((value=vec[cur_col]) != D_ZERO) |
2084 | mtx_create_element_value(mtx,org_row,toorg[cur_col],value); |
2085 | } |
2086 | } |
2087 | |
2088 | void mtx_fill_cur_col_vec(mtx_matrix_t mtx, int32 col, |
2089 | real64 *vec, mtx_range_t *rng) |
2090 | { |
2091 | int32 org_col,highrow,lowrow, *toorg; |
2092 | register int32 cur_row; |
2093 | register real64 value; |
2094 | |
2095 | #if MTX_DEBUG |
2096 | if(!mtx_check_matrix(mtx)) return; |
2097 | #endif |
2098 | |
2099 | org_col = mtx->perm.col.cur_to_org[col]; |
2100 | toorg=mtx->perm.row.cur_to_org; |
2101 | if( rng == mtx_ALL_ROWS ) { |
2102 | highrow=mtx->order-1; |
2103 | lowrow=0; |
2104 | } else { |
2105 | highrow= rng->high; |
2106 | lowrow=rng->low; |
2107 | } |
2108 | for( cur_row=lowrow ; cur_row <= highrow ; cur_row++ ) { |
2109 | if ((value=vec[cur_row]) != D_ZERO) |
2110 | mtx_create_element_value(mtx,toorg[cur_row],org_col,value); |
2111 | } |
2112 | } |
2113 | |
2114 | void mtx_dropfill_cur_col_vec(mtx_matrix_t mtx, int32 col, |
2115 | real64 *vec, mtx_range_t *rng, |
2116 | real64 tol) |
2117 | { |
2118 | int32 org_col,highrow,lowrow, *toorg; |
2119 | register int32 cur_row; |
2120 | register real64 value; |
2121 | |
2122 | if (tol==0.0) { |
2123 | mtx_fill_cur_col_vec(mtx,col,vec,rng); |
2124 | return; |
2125 | } |
2126 | |
2127 | #if MTX_DEBUG |
2128 | if(!mtx_check_matrix(mtx)) return; |
2129 | #endif |
2130 | |
2131 | org_col = mtx->perm.col.cur_to_org[col]; |
2132 | toorg=mtx->perm.row.cur_to_org; |
2133 | if( rng == mtx_ALL_ROWS ) { |
2134 | highrow=mtx->order-1; |
2135 | lowrow=0; |
2136 | } else { |
2137 | highrow= rng->high; |
2138 | lowrow=rng->low; |
2139 | } |
2140 | for( cur_row=lowrow ; cur_row <= highrow ; cur_row++ ) { |
2141 | if (fabs(value=vec[cur_row]) > tol) |
2142 | mtx_create_element_value(mtx,toorg[cur_row],org_col,value); |
2143 | } |
2144 | } |
2145 | |
2146 | void mtx_dropfill_cur_row_vec(mtx_matrix_t mtx, int32 row, |
2147 | real64 *vec, mtx_range_t *rng, |
2148 | real64 tol) |
2149 | { |
2150 | int32 org_row,highcol,lowcol, *toorg; |
2151 | register int32 cur_col; |
2152 | register real64 value; |
2153 | |
2154 | if (tol==0.0) { |
2155 | mtx_fill_cur_row_vec(mtx,row,vec,rng); |
2156 | return; |
2157 | } |
2158 | |
2159 | #if MTX_DEBUG |
2160 | if(!mtx_check_matrix(mtx)) return; |
2161 | #endif |
2162 | |
2163 | org_row = mtx->perm.row.cur_to_org[row]; |
2164 | toorg=mtx->perm.col.cur_to_org; |
2165 | if( rng == mtx_ALL_COLS ) { |
2166 | highcol = mtx->order-1; |
2167 | lowcol = 0; |
2168 | } else { |
2169 | highcol= rng->high; |
2170 | lowcol=rng->low; |
2171 | } |
2172 | for( cur_col=lowcol ; cur_col <= highcol ; cur_col++ ) { |
2173 | if (fabs(value=vec[cur_col]) > tol) |
2174 | mtx_create_element_value(mtx,org_row,toorg[cur_col],value); |
2175 | } |
2176 | } |
2177 | |
2178 | void mtx_fill_org_row_sparse(mtx_matrix_t mtx, int32 row, |
2179 | const mtx_sparse_t *sp) |
2180 | { |
2181 | int32 orgrow,i; |
2182 | #if MTX_DEBUG |
2183 | if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return; |
2184 | #endif |
2185 | orgrow = mtx->perm.row.cur_to_org[row]; |
2186 | for (i=0; i < sp->len; i++) { |
2187 | if (sp->data[i] != D_ZERO |
2188 | #if MTX_DEBUG |
2189 | && sp->idata[i] >=0 && sp->idata[i] < mtx->order |
2190 | #endif |
2191 | ) { |
2192 | mtx_create_element_value(mtx,orgrow,sp->idata[i],sp->data[i]); |
2193 | } |
2194 | } |
2195 | } |
2196 | |
2197 | void mtx_fill_org_col_sparse(mtx_matrix_t mtx, int32 col, |
2198 | const mtx_sparse_t *sp) |
2199 | { |
2200 | int32 orgcol,i; |
2201 | #if MTX_DEBUG |
2202 | if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return; |
2203 | #endif |
2204 | orgcol = mtx->perm.col.cur_to_org[col]; |
2205 | for (i=0; i < sp->len; i++) { |
2206 | if (sp->data[i] != D_ZERO |
2207 | #if MTX_DEBUG |
2208 | && sp->idata[i] >=0 && sp->idata[i] < mtx->order |
2209 | #endif |
2210 | ) { |
2211 | mtx_create_element_value(mtx,sp->idata[i],orgcol,sp->data[i]); |
2212 | } |
2213 | } |
2214 | } |
2215 | |
2216 | void mtx_fill_cur_row_sparse(mtx_matrix_t mtx, int32 row, |
2217 | const mtx_sparse_t *sp) |
2218 | { |
2219 | int32 orgrow,i; |
2220 | #if MTX_DEBUG |
2221 | if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return; |
2222 | #endif |
2223 | orgrow = mtx->perm.row.cur_to_org[row]; |
2224 | for (i=0; i < sp->len; i++) { |
2225 | if ( |
2226 | #if MTX_DEBUG |
2227 | sp->idata[i] >=0 && sp->idata[i] < mtx->order && |
2228 | #endif |
2229 | sp->data[i] != D_ZERO) { |
2230 | mtx_create_element_value(mtx,orgrow, |
2231 | mtx->perm.col.cur_to_org[sp->idata[i]], |
2232 | sp->data[i]); |
2233 | } |
2234 | } |
2235 | } |
2236 | |
2237 | void mtx_fill_cur_col_sparse(mtx_matrix_t mtx, int32 col, |
2238 | const mtx_sparse_t *sp) |
2239 | { |
2240 | int32 orgcol,i; |
2241 | #if MTX_DEBUG |
2242 | if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return; |
2243 | #endif |
2244 | orgcol = mtx->perm.col.cur_to_org[col]; |
2245 | for (i=0; i < sp->len; i++) { |
2246 | if ( |
2247 | #if MTX_DEBUG |
2248 | sp->idata[i] >=0 && sp->idata[i] < mtx->order && |
2249 | #endif |
2250 | sp->data[i] != D_ZERO) { |
2251 | mtx_create_element_value(mtx, |
2252 | mtx->perm.col.cur_to_org[sp->idata[i]], |
2253 | orgcol, |
2254 | sp->data[i]); |
2255 | } |
2256 | } |
2257 | } |
2258 | |
2259 | void mtx_mult_row(mtx_matrix_t mtx, int32 row, real64 factor, mtx_range_t *rng) |
2260 | { |
2261 | struct element_t *elt; |
2262 | int32 *tocur; |
2263 | |
2264 | if( factor == D_ZERO ) { |
2265 | mtx_clear_row(mtx,row,rng); |
2266 | return; |
2267 | } |
2268 | |
2269 | if (factor == D_ONE) return; |
2270 | #if MTX_DEBUG |
2271 | if(!mtx_check_matrix(mtx)) return; |
2272 | #endif |
2273 | |
2274 | tocur = mtx->perm.col.org_to_cur; |
2275 | elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]]; |
2276 | if( rng == mtx_ALL_COLS ) |
2277 | for( ; NOTNULL(elt); elt = elt->next.col ) elt->value *= factor; |
2278 | else if( rng->high >= rng->low ) |
2279 | for( ; NOTNULL(elt); elt = elt->next.col ) |
2280 | if( in_range(rng,tocur[elt->col]) ) elt->value *= factor; |
2281 | } |
2282 | |
2283 | void mtx_mult_col(mtx_matrix_t mtx, int32 col, real64 factor, mtx_range_t *rng) |
2284 | { |
2285 | struct element_t *elt; |
2286 | int32 *tocur; |
2287 | |
2288 | if( factor == D_ZERO ) { |
2289 | mtx_clear_col(mtx,col,rng); |
2290 | return; |
2291 | } |
2292 | if (factor == D_ONE) return; |
2293 | |
2294 | #if MTX_DEBUG |
2295 | if(!mtx_check_matrix(mtx)) return; |
2296 | #endif |
2297 | |
2298 | tocur = mtx->perm.row.org_to_cur; |
2299 | elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]]; |
2300 | if( rng == mtx_ALL_ROWS ) |
2301 | for( ; NOTNULL(elt); elt = elt->next.row ) elt->value *= factor; |
2302 | else if( rng->high >= rng->low ) |
2303 | for( ; NOTNULL(elt); elt = elt->next.row ) |
2304 | if( in_range(rng,tocur[elt->row]) ) elt->value *= factor; |
2305 | } |
2306 | |
2307 | void mtx_mult_row_zero(mtx_matrix_t mtx, int32 row, mtx_range_t *rng) |
2308 | { |
2309 | struct element_t *elt; |
2310 | int32 *tocur; |
2311 | |
2312 | #if MTX_DEBUG |
2313 | if(!mtx_check_matrix(mtx)) return; |
2314 | #endif |
2315 | |
2316 | tocur = mtx->perm.col.org_to_cur; |
2317 | elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]]; |
2318 | if( rng == mtx_ALL_COLS ) |
2319 | for( ; NOTNULL(elt); elt = elt->next.col ) elt->value = D_ZERO; |
2320 | else if( rng->high >= rng->low ) |
2321 | for( ; NOTNULL(elt); elt = elt->next.col ) |
2322 | if( in_range(rng,tocur[elt->col]) ) elt->value = D_ZERO; |
2323 | } |
2324 | |
2325 | void mtx_mult_col_zero(mtx_matrix_t mtx, int32 col, mtx_range_t *rng) |
2326 | { |
2327 | struct element_t *elt; |
2328 | int32 *tocur; |
2329 | |
2330 | #if MTX_DEBUG |
2331 | if(!mtx_check_matrix(mtx)) return; |
2332 | #endif |
2333 | |
2334 | tocur = mtx->perm.row.org_to_cur; |
2335 | elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]]; |
2336 | if( rng == mtx_ALL_ROWS ) |
2337 | for( ; NOTNULL(elt); elt = elt->next.row ) elt->value = D_ZERO; |
2338 | else if( rng->high >= rng->low ) |
2339 | for( ; NOTNULL(elt); elt = elt->next.row ) |
2340 | if( in_range(rng,tocur[elt->row]) ) elt->value = D_ZERO; |
2341 | } |
2342 | |
2343 | void mtx_add_row(mtx_matrix_t mtx, int32 s_cur, int32 t_cur, real64 factor, |
2344 | mtx_range_t *rng) |
2345 | { |
2346 | register int32 org_col; |
2347 | int32 t_org,*tocur; |
2348 | struct element_t **arr,*elt; |
2349 | |
2350 | #if MTX_DEBUG |
2351 | if(!mtx_check_matrix(mtx)) return; |
2352 | #endif |
2353 | |
2354 | t_org = mtx->perm.row.cur_to_org[t_cur]; |
2355 | if( rng == mtx_ALL_COLS ) { |
2356 | arr = mtx_expand_row(mtx,t_org); /* Expand the target row */ |
2357 | elt = mtx->hdr.row[mtx->perm.row.cur_to_org[s_cur]]; |
2358 | for( ; NOTNULL(elt); elt = elt->next.col ) { |
2359 | if( ISNULL(arr[(org_col=elt->col)]) ) { |
2360 | arr[org_col] = |
2361 | mtx_create_element_value(mtx,t_org,org_col,(factor * elt->value)); |
2362 | } else { |
2363 | arr[org_col]->value += factor * elt->value; |
2364 | } |
2365 | } |
2366 | mtx_renull_using_row(mtx,t_org,arr); |
2367 | } else if( rng->high >= rng->low ) { |
2368 | register int32 cur_col; |
2369 | |
2370 | tocur = mtx->perm.col.org_to_cur; |
2371 | arr = mtx_expand_row(mtx,t_org); /* Expand the target row */ |
2372 | elt = mtx->hdr.row[mtx->perm.row.cur_to_org[s_cur]]; |
2373 | for( ; NOTNULL(elt); elt = elt->next.col ) { |
2374 | cur_col=tocur[(org_col=elt->col)]; |
2375 | if( in_range(rng,cur_col) ) |
2376 | if( NOTNULL(arr[org_col]) ) { |
2377 | arr[org_col]->value += factor * elt->value; |
2378 | } else { |
2379 | arr[org_col] = |
2380 | mtx_create_element_value(mtx,t_org,org_col, |
2381 | (factor * elt->value)); |
2382 | /* hit rate usually lower */ |
2383 | } |
2384 | } |
2385 | mtx_renull_using_row(mtx,t_org,arr); |
2386 | } |
2387 | mtx_null_vector_release(); |
2388 | } |
2389 | |
2390 | void mtx_add_col(mtx_matrix_t mtx, int32 s_cur, int32 t_cur, real64 factor, |
2391 | mtx_range_t *rng) |
2392 | { |
2393 | int32 t_org,*tocur; |
2394 | struct element_t **arr,*elt; |
2395 | |
2396 | #if MTX_DEBUG |
2397 | if(!mtx_check_matrix(mtx)) return; |
2398 | #endif |
2399 | |
2400 | if( rng == mtx_ALL_ROWS ) { |
2401 | t_org = mtx->perm.col.cur_to_org[t_cur]; |
2402 | arr = mtx_expand_col(mtx,t_org); /* Expand the target col */ |
2403 | elt = mtx->hdr.col[mtx->perm.col.cur_to_org[s_cur]]; |
2404 | for( ; NOTNULL(elt); elt = elt->next.row ) |
2405 | if( ISNULL(arr[elt->row]) ) { |
2406 | arr[elt->row] = |
2407 | mtx_create_element_value(mtx,elt->row,t_org, |
2408 | (factor * elt->value)); |
2409 | } else arr[elt->row]->value += factor * elt->value; |
2410 | mtx_renull_using_col(mtx,t_org,arr); |
2411 | } else if( rng->high >= rng->low ) { |
2412 | tocur = mtx->perm.row.org_to_cur; |
2413 | t_org = mtx->perm.col.cur_to_org[t_cur]; |
2414 | arr = mtx_expand_col(mtx,t_org); /* Expand the target col */ |
2415 | elt = mtx->hdr.col[mtx->perm.col.cur_to_org[s_cur]]; |
2416 | for( ; NOTNULL(elt); elt = elt->next.row ) |
2417 | if( in_range(rng,tocur[elt->row]) ) |
2418 | if( NOTNULL(arr[elt->row]) ) { |
2419 | arr[elt->row]->value += factor * elt->value; |
2420 | } else { |
2421 | arr[elt->row] = |
2422 | mtx_create_element_value(mtx,elt->row,t_org, |
2423 | (factor * elt->value)); |
2424 | } |
2425 | mtx_renull_using_col(mtx,t_org,arr); |
2426 | } |
2427 | mtx_null_vector_release(); |
2428 | } |
2429 | |
2430 | static struct element_t **mtx_expand_row_series( mtx_matrix_t mtx, int32 org) |
2431 | /** |
2432 | *** Expands the given row into an array of pointers, indexed on original |
2433 | *** col number. The array is obtained from mtx_null_row_vector(). |
2434 | *** Be sure to call mtx_null_row_vector_release() when done with the vector and |
2435 | *** you have rezeroed it. |
2436 | **/ |
2437 | { |
2438 | struct element_t **arr; |
2439 | struct element_t *elt; |
2440 | |
2441 | arr = mtx_null_row_vector(mtx->order); |
2442 | for( elt=mtx->hdr.row[org]; NOTNULL(elt); elt=elt->next.col ) |
2443 | arr[elt->col] = elt; |
2444 | return(arr); |
2445 | } |
2446 | |
2447 | static struct element_t **mtx_expand_col_series(mtx_matrix_t mtx, int32 org) |
2448 | /** |
2449 | *** Expands the given col into an array of pointers, indexed on original |
2450 | *** row number. The array is obtained from mtx_null_col_vector(). |
2451 | *** Be sure to call mtx_null_col_vector_release() when done with the vector and |
2452 | *** you have rezeroed it. |
2453 | **/ |
2454 | { |
2455 | struct element_t **arr; |
2456 | struct element_t *elt; |
2457 | |
2458 | arr = mtx_null_col_vector(mtx->order); |
2459 | for( elt = mtx->hdr.col[org] ; NOTNULL(elt) ; elt = elt->next.row ) |
2460 | arr[elt->row] = elt; |
2461 | return(arr); |
2462 | } |
2463 | |
2464 | struct add_series_data { |
2465 | mtx_matrix_t mtx; /* matrix we're operating on */ |
2466 | struct element_t **arr; /* target row/col expansion array */ |
2467 | int32 *tocur; /* col/row permutation vector */ |
2468 | int32 t_org; /* org row/col which is expanded */ |
2469 | }; |
2470 | |
2471 | static struct add_series_data |
2472 | rsdata={NULL,NULL,NULL,mtx_NONE}, |
2473 | csdata={NULL,NULL,NULL,mtx_NONE}; |
2474 | |
2475 | static void add_series_data_release(void) |
2476 | { |
2477 | /* if apparently sane, and have array release the arr */ |
2478 | if (NOTNULL(rsdata.mtx) && NOTNULL(rsdata.tocur) |
2479 | && rsdata.t_org >= 0 && NOTNULL(rsdata.arr)) { |
2480 | mtx_null_row_vector_release(); |
2481 | } |
2482 | if (NOTNULL(csdata.mtx) && NOTNULL(csdata.tocur) |
2483 | && csdata.t_org >= 0 && NOTNULL(csdata.arr)) { |
2484 | mtx_null_col_vector_release(); |
2485 | } |
2486 | /* reinit data */ |
2487 | |
2488 | rsdata.mtx = NULL; |
2489 | rsdata.arr = NULL; |
2490 | rsdata.tocur = NULL; |
2491 | rsdata.t_org = mtx_NONE; |
2492 | |
2493 | csdata.mtx = NULL; |
2494 | csdata.arr = NULL; |
2495 | csdata.tocur = NULL; |
2496 | csdata.t_org = mtx_NONE; |
2497 | } |
2498 | |
2499 | void mtx_add_row_series_init(mtx_matrix_t mtx,int32 t_cur,boolean use) |
2500 | { |
2501 | |
2502 | #if MTX_DEBUG |
2503 | if(!mtx_check_matrix(mtx)) return; |
2504 | #endif |
2505 | |
2506 | if ( !(rsdata.mtx) ) { |
2507 | /* this is a grabbing call due to memory reuse */ |
2508 | if ( mtx && (t_cur >= 0) && (t_cur < mtx->order) ) { |
2509 | rsdata.mtx = mtx; |
2510 | rsdata.tocur = mtx->perm.col.org_to_cur; |
2511 | rsdata.t_org = mtx->perm.row.cur_to_org[t_cur]; |
2512 | rsdata.arr = mtx_expand_row_series(mtx,rsdata.t_org); |
2513 | } else { |
2514 | FPRINTF(g_mtxerr,"ERROR: mtx_add_row_series_init called\n"); |
2515 | FPRINTF(g_mtxerr," for grab with invalid column or mtx.\n"); |
2516 | FPRINTF(g_mtxerr," col number %d.\n",t_cur); |
2517 | } |
2518 | return; |
2519 | } else { |
2520 | /* this is supposed to be a releasing call */ |
2521 | if (t_cur != mtx_NONE) { |
2522 | FPRINTF(g_mtxerr,"ERROR: mtx_add_row_series_init called for\n"); |
2523 | FPRINTF(g_mtxerr," release without mtx_NONE.\n"); |
2524 | return; |
2525 | } |
2526 | if (mtx != rsdata.mtx) { |
2527 | FPRINTF(g_mtxerr,"ERROR: mtx_add_row_series_init called for\n"); |
2528 | FPRINTF(g_mtxerr," release with ungrabbed matrix.\n"); |
2529 | return; |
2530 | } |
2531 | if (use) { |
2532 | mtx_renull_using_row(rsdata.mtx, rsdata.t_org, rsdata.arr); |
2533 | } else { |
2534 | mtx_renull_all(rsdata.mtx, rsdata.arr); |
2535 | } |
2536 | mtx_null_row_vector_release(); |
2537 | rsdata.mtx = NULL; |
2538 | rsdata.arr = NULL; |
2539 | rsdata.tocur = NULL; |
2540 | rsdata.t_org = mtx_NONE; |
2541 | } |
2542 | } |
2543 | |
2544 | void mtx_add_row_series(int32 s_cur,real64 factor,mtx_range_t *rng) |
2545 | { |
2546 | register int32 org_col; |
2547 | int32 t_org,*tocur; |
2548 | struct element_t **arr,*elt; |
2549 | mtx_matrix_t mtx; |
2550 | |
2551 | if ( !(rsdata.mtx) ) { |
2552 | FPRINTF(g_mtxerr,"ERROR: mtx_add_row_series called for\n"); |
2553 | FPRINTF(g_mtxerr," without grabbing target row first.\n"); |
2554 | return; |
2555 | } |
2556 | mtx = rsdata.mtx; |
2557 | arr = rsdata.arr; |
2558 | tocur = rsdata.tocur; |
2559 | t_org = rsdata.t_org; |
2560 | |
2561 | #if MTX_DEBUG |
2562 | if(!mtx_check_matrix(mtx)) return; |
2563 | #endif |
2564 | |
2565 | elt = mtx->hdr.row[mtx->perm.row.cur_to_org[s_cur]]; |
2566 | if( rng == mtx_ALL_COLS ) { |
2567 | for( ; NOTNULL(elt); elt = elt->next.col ) { |
2568 | if( ISNULL(arr[(org_col=elt->col)]) ) { |
2569 | arr[org_col] = |
2570 | mtx_create_element_value(mtx,t_org,org_col,(factor * elt->value)); |
2571 | } else { |
2572 | arr[org_col]->value += factor * elt->value; |
2573 | } |
2574 | } |
2575 | return; |
2576 | } |
2577 | /* fast_in_range is a 10% winner on the alpha, and should be even more |
2578 | on the other platforms. */ |
2579 | if( rng->high >= rng->low ) { |
2580 | register int32 cur_col, lo,hi; |
2581 | lo=rng->low; hi= rng->high; |
2582 | for( ; NOTNULL(elt); elt = elt->next.col ) { |
2583 | cur_col=tocur[(org_col=elt->col)]; |
2584 | if( fast_in_range(lo,hi,cur_col) ) { |
2585 | if( NOTNULL(arr[org_col]) ) { |
2586 | arr[org_col]->value += factor * elt->value; |
2587 | } else { |
2588 | arr[org_col] = |
2589 | mtx_create_element_value(mtx,t_org,org_col,(factor*elt->value)); |
2590 | } |
2591 | } |
2592 | } |
2593 | } |
2594 | } |
2595 | |
2596 | void mtx_add_col_series_init(mtx_matrix_t mtx,int32 t_cur,boolean use) |
2597 | { |
2598 | |
2599 | #if MTX_DEBUG |
2600 | if(!mtx_check_matrix(mtx)) return; |
2601 | #endif |
2602 | |
2603 | if ( !(csdata.mtx) ) { |
2604 | /* this is a grabbing call */ |
2605 | if ( mtx && (t_cur >= 0) && (t_cur < mtx->order) ) { |
2606 | csdata.mtx = mtx; |
2607 | csdata.tocur = mtx->perm.row.org_to_cur; |
2608 | csdata.t_org = mtx->perm.col.cur_to_org[t_cur]; |
2609 | csdata.arr = mtx_expand_col_series(mtx,csdata.t_org); |
2610 | } else { |
2611 | FPRINTF(g_mtxerr,"ERROR: mtx_add_col_series_init called\n"); |
2612 | FPRINTF(g_mtxerr," for grab with invalid row or mtx.\n"); |
2613 | FPRINTF(g_mtxerr," row number %d.\n",t_cur); |
2614 | } |
2615 | return; |
2616 | } else { |
2617 | /* this is supposed to be a releasing call */ |
2618 | if (t_cur != mtx_NONE) { |
2619 | FPRINTF(g_mtxerr,"ERROR: mtx_add_col_series_init called for\n"); |
2620 | FPRINTF(g_mtxerr," release without mtx_NONE.\n"); |
2621 | return; |
2622 | } |
2623 | if (mtx != csdata.mtx) { |
2624 | FPRINTF(g_mtxerr,"ERROR: mtx_add_col_series_init called for\n"); |
2625 | FPRINTF(g_mtxerr," release with ungrabbed matrix.\n"); |
2626 | return; |
2627 | } |
2628 | if (use) { |
2629 | mtx_renull_using_col(csdata.mtx, csdata.t_org, csdata.arr); |
2630 | } else { |
2631 | mtx_renull_all(csdata.mtx, csdata.arr); |
2632 | } |
2633 | mtx_null_col_vector_release(); |
2634 | csdata.mtx = NULL; |
2635 | csdata.arr = NULL; |
2636 | csdata.tocur = NULL; |
2637 | csdata.t_org = mtx_NONE; |
2638 | } |
2639 | } |
2640 | |
2641 | void mtx_add_col_series(int32 s_cur,real64 factor,mtx_range_t *rng) |
2642 | { |
2643 | register int32 org_row; |
2644 | int32 t_org,*tocur; |
2645 | struct element_t **arr,*elt; |
2646 | mtx_matrix_t mtx; |
2647 | |
2648 | if ( !(csdata.mtx) ) { |
2649 | FPRINTF(g_mtxerr,"ERROR: mtx_add_col_series called for\n"); |
2650 | FPRINTF(g_mtxerr," without grabbing target col first.\n"); |
2651 | return; |
2652 | } |
2653 | mtx = csdata.mtx; |
2654 | arr = csdata.arr; |
2655 | tocur = csdata.tocur; |
2656 | t_org = csdata.t_org; |
2657 | |
2658 | #if MTX_DEBUG |
2659 | if(!mtx_check_matrix(mtx)) return; |
2660 | #endif |
2661 | |
2662 | elt = mtx->hdr.col[mtx->perm.col.cur_to_org[s_cur]]; |
2663 | if( rng == mtx_ALL_ROWS ) { |
2664 | for( ; NOTNULL(elt); elt = elt->next.row ) { |
2665 | if( ISNULL(arr[(org_row=elt->row)]) ) { |
2666 | arr[org_row] = |
2667 | mtx_create_element_value(mtx,org_row,t_org,(factor * elt->value)); |
2668 | } else { |
2669 | arr[org_row]->value += factor * elt->value; |
2670 | } |
2671 | } |
2672 | return; |
2673 | } |
2674 | if( rng->high >= rng->low ) { |
2675 | register int32 cur_row; |
2676 | |
2677 | for( ; NOTNULL(elt); elt = elt->next.row ) { |
2678 | cur_row=tocur[(org_row=elt->row)]; |
2679 | if( in_range(rng,cur_row) ) { |
2680 | if( NOTNULL(arr[org_row]) ) { |
2681 | arr[org_row]->value += factor * elt->value; |
2682 | } else { |
2683 | arr[org_row] = |
2684 | mtx_create_element_value(mtx,org_row,t_org,(factor*elt->value)); |
2685 | } |
2686 | } |
2687 | } |
2688 | } |
2689 | } |
2690 | |
2691 | void mtx_old_add_row_sparse(mtx_matrix_t mtx, |
2692 | int32 t_cur, /* cur index of target row */ |
2693 | real64 *s_cur, /* dense source row, curindexed */ |
2694 | real64 factor, /* coefficient */ |
2695 | mtx_range_t *rng, |
2696 | int32 *ilist) /* list to add over */ |
2697 | { |
2698 | int32 t_org,*toorg,cindex,orgcindex,hilim; |
2699 | struct element_t **arr; |
2700 | real64 value; |
2701 | |
2702 | #if MTX_DEBUG |
2703 | if(!mtx_check_matrix(mtx)) return; |
2704 | #endif |
2705 | |
2706 | if (factor==D_ZERO) return; /* adding 0 rather silly */ |
2707 | |
2708 | if( rng == mtx_ALL_COLS ) { |
2709 | t_org = mtx->perm.row.cur_to_org[t_cur]; /* org of target row */ |
2710 | arr = mtx_expand_row(mtx,t_org); /* Expand the target row */ |
2711 | toorg = mtx->perm.col.cur_to_org; /* current col perm */ |
2712 | hilim=mtx->order; |
2713 | |
2714 | if (!ilist) { /* FULL ROW CASE no ilist */ |
2715 | for (cindex=0; cindex< hilim; cindex++) { |
2716 | if ((value=s_cur[cindex])!=D_ZERO) { |
2717 | if ( ISNULL(arr[(orgcindex=toorg[cindex])]) ) { /* add element */ |
2718 | arr[orgcindex] = |
2719 | mtx_create_element_value(mtx,t_org,orgcindex,(factor * value)); |
2720 | } else { /* increment element */ |
2721 | arr[orgcindex]->value += factor * value; |
2722 | } |
2723 | } |
2724 | } |
2725 | } else { /* SPARSE ROW CASE with ilist */ |
2726 | int32 i; |
2727 | i=0; |
2728 | while ((cindex=ilist[i])>=0) { |
2729 | value=s_cur[cindex]; |
2730 | if ( ISNULL(arr[(orgcindex=toorg[cindex])]) ) { /* add element */ |
2731 | arr[orgcindex] = |
2732 | mtx_create_element_value(mtx,t_org,orgcindex,(factor * value)); |
2733 | } else { /* increment element */ |
2734 | arr[orgcindex]->value += factor * value; |
2735 | } |
2736 | i++; |
2737 | } |
2738 | } |
2739 | mtx_renull_using_row(mtx,t_org,arr); |
2740 | } else if( rng->high >= rng->low ) { /* DENSE RANGE CASE */ |
2741 | t_org = mtx->perm.row.cur_to_org[t_cur]; /* org of target row */ |
2742 | arr = mtx_expand_row(mtx,t_org); /* Expand the target row */ |
2743 | toorg = mtx->perm.col.cur_to_org; /* current col perm */ |
2744 | hilim = rng->high; |
2745 | |
2746 | for (cindex=rng->low; cindex<= hilim; cindex++) { |
2747 | if ((value=s_cur[cindex])!=D_ZERO) { |
2748 | if ( ISNULL(arr[(orgcindex=toorg[cindex])]) ) { /* add element */ |
2749 | arr[orgcindex] = |
2750 | mtx_create_element_value(mtx,t_org,orgcindex,(factor * value)); |
2751 | } else { /* increment element */ |
2752 | arr[orgcindex]->value += factor * value; |
2753 | } |
2754 | } |
2755 | } |
2756 | mtx_renull_using_row(mtx,t_org,arr); |
2757 | } |
2758 | mtx_null_vector_release(); |
2759 | } |
2760 | |
2761 | void mtx_old_add_col_sparse(mtx_matrix_t mtx, |
2762 | int32 t_cur, /* cur index of target col */ |
2763 | real64 *s_cur, /* dense source col, curindexed */ |
2764 | real64 factor, /* coefficient */ |
2765 | mtx_range_t *rng, /* range to add over or */ |
2766 | int32 *ilist) /* list to add over */ |
2767 | { |
2768 | int32 t_org,*toorg,rindex,orgrindex,hilim; |
2769 | struct element_t **arr; |
2770 | real64 value; |
2771 | |
2772 | #if MTX_DEBUG |
2773 | if(!mtx_check_matrix(mtx)) return; |
2774 | #endif |
2775 | |
2776 | if (factor==D_ZERO) return; |
2777 | |
2778 | if( rng == mtx_ALL_ROWS ) { |
2779 | t_org = mtx->perm.col.cur_to_org[t_cur]; /* org of target col */ |
2780 | arr = mtx_expand_col(mtx,t_org); /* Expand the target col */ |
2781 | toorg = mtx->perm.row.cur_to_org; /* current row perm */ |
2782 | hilim=mtx->order; |
2783 | |
2784 | if (!ilist) { /* FULL COL CASE no ilist */ |
2785 | for (rindex=0; rindex< hilim; rindex++) { |
2786 | if ((value=s_cur[rindex])!=D_ZERO) { |
2787 | if ( ISNULL(arr[(orgrindex=toorg[rindex])]) ) { /* add element */ |
2788 | arr[orgrindex] = |
2789 | mtx_create_element_value(mtx,orgrindex,t_org,(factor * value)); |
2790 | } else { /* increment element */ |
2791 | arr[orgrindex]->value += factor * value; |
2792 | } |
2793 | } |
2794 | } |
2795 | } else { /* SPARSE COL CASE with ilist */ |
2796 | int32 i; |
2797 | i=0; |
2798 | while ((rindex=ilist[i])>=0) { |
2799 | value=s_cur[rindex]; |
2800 | if ( ISNULL(arr[(orgrindex=toorg[rindex])]) ) { /* add element */ |
2801 | arr[orgrindex] = |
2802 | mtx_create_element_value(mtx,orgrindex,t_org,(factor * value)); |
2803 | } else { /* increment element */ |
2804 | arr[orgrindex]->value += factor * value; |
2805 | } |
2806 | i++; |
2807 | } |
2808 | } |
2809 | mtx_renull_using_col(mtx,t_org,arr); |
2810 | |
2811 | } else if( rng->high >= rng->low ) { /* DENSE RANGE CASE */ |
2812 | t_org = mtx->perm.col.cur_to_org[t_cur]; /* org of target col */ |
2813 | arr = mtx_expand_col(mtx,t_org); /* Expand the target col */ |
2814 | toorg = mtx->perm.row.cur_to_org; /* current row perm */ |
2815 | hilim = rng->high; |
2816 | |
2817 | for (rindex=rng->low; rindex<= hilim; rindex++) { |
2818 | if ((value=s_cur[rindex])!=D_ZERO) { |
2819 | if ( ISNULL(arr[(orgrindex=toorg[rindex])]) ) { /* add element */ |
2820 | arr[orgrindex] = |
2821 | mtx_create_element_value(mtx,orgrindex,t_org,(factor * value)); |
2822 | } else { /* increment element */ |
2823 | arr[orgrindex]->value += factor * value; |
2824 | } |
2825 | } |
2826 | } |
2827 | mtx_renull_using_col(mtx,t_org,arr); |
2828 | } |
2829 | mtx_null_vector_release(); |
2830 | } |
2831 | |
2832 | size_t mtx_size(mtx_matrix_t mtx) { |
2833 | size_t size=0; |
2834 | if (ISNULL(mtx)) return size; |
2835 | size += (1+mtx->nslaves)*sizeof(struct mtx_header); |
2836 | if (ISSLAVE(mtx)) { |
2837 | return 2*sizeof(struct element_t *)*(size_t)mtx->capacity +size; |
2838 | } |
2839 | size += mem_sizeof_store(mtx->ms); |
2840 | /* headers */ |
2841 | size += (1 + mtx->nslaves) * (size_t)mtx->capacity * |
2842 | (size_t)2 * sizeof(struct element_t *); |
2843 | /* block data */ |
2844 | if (mtx->data->nblocks >0) |
2845 | size += sizeof(mtx_region_t)*(size_t)(mtx->data->nblocks-1); |
2846 | /* permutations */ |
2847 | size += (size_t)4*(sizeof(int32)*(size_t)mtx->capacity); |
2848 | return size; |
2849 | } |
2850 | |
2851 | size_t mtx_chattel_size(mtx_matrix_t mtx) { |
2852 | size_t size=0; |
2853 | int32 i; |
2854 | if (ISNULL(mtx) || ISSLAVE(mtx)) return size; |
2855 | /* headers */ |
2856 | size += (mtx->nslaves)*sizeof(struct mtx_header); |
2857 | /* incidence */ |
2858 | for (i=0; i <mtx->nslaves; i++) { |
2859 | size += (size_t)mtx_nonzeros_in_region(mtx->slaves[i],mtx_ENTIRE_MATRIX); |
2860 | } |
2861 | /*nz headers */ |
2862 | size += mtx->nslaves * (size_t)mtx->capacity * |
2863 | ( |