/[ascend]/trunk/ascend4/solver/mtx_basic.c
ViewVC logotype

Annotation of /trunk/ascend4/solver/mtx_basic.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 106688 byte(s)
Setting up web subdirectory in repository
1 aw0a 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