ascend4/solver/mtx_use_only.c

/*
 *  mtx: Ascend Sparse Matrix Package
 *  by Karl Michael Westerberg
 *  Created: 2/6/90
 *  Version: $Revision: 1.9 $
 *  Version control file: $RCSfile: mtx_use_only.c,v $
 *  Date last modified: $Date: 2000/01/25 02:27:12 $
 *  Last modified by: $Author: ballan $
 *
 *  This file is part of the SLV solver.
 *
 *  Copyright (C) 1990 Karl Michael Westerberg
 *  Copyright (C) 1993 Joseph Zaher
 *  Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
 *  Copyright (C) 1995 Kirk Andre Abbott, Benjamin Andrew Allan
 *  Copyright (C) 1996 Benjamin Andrew Allan
 *
 *  The SLV solver is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public License as
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version.
 *
 *  The SLV solver is distributed in hope that it will be
 *  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along with
 *  the program; if not, write to the Free Software Foundation, Inc., 675
 *  Mass Ave, Cambridge, MA 02139 USA.  Check the file named COPYING.
 *  COPYING is found in ../compiler.
 */
#include "utilities/ascConfig.h"
#include "utilities/ascMalloc.h"
#include "utilities/mem.h"
#include "solver/mtx.h"
/* grab our private parts */
#define __MTX_C_SEEN__
#include "solver/mtx_use_only.h"

#define free_unless_null(ptr) if( NOTNULL(ptr) ) ascfree(ptr)

FILE *g_mtxerr;

struct element_t *mtx_find_element(mtx_matrix_t mtx, int32 org_row,
                                   int32 org_col)
{
   register struct element_t *elt;

/* old way.
   for( elt = mtx->hdr.row[org_row] ; NOTNULL(elt) ; elt = elt->next.col )
      if( elt->col == org_col )
         break;
*/
   for( elt = mtx->hdr.row[org_row] ;
        NOTNULL(elt) && elt->col!=org_col ;
        elt = elt->next.col );
   return(elt);
}

struct element_t *mtx_create_element(mtx_matrix_t mtx,
                                     int32 org_row, int32 org_col)
/**
 ***  Creates the given element and returns a pointer to it.  The value is
 ***  initially zero.
 **/
{
  register struct element_t *elt;

#if MTX_DEBUG
  if( NOTNULL(mtx_find_element(mtx,org_row,org_col)) ) {
    FPRINTF(g_mtxerr,"ERROR:  (mtx) create_element\n");
    FPRINTF(g_mtxerr,"        Element (%d,%d) already exists.\n",
            org_row,org_col);
  }
#endif
  elt = (struct element_t *)mem_get_element(mtx->ms);
  /* guess who didn't check for the alloc return. */
  elt->value = 0.0;
  elt->row = org_row;
  elt->col = org_col;
  elt->next.col = mtx->hdr.row[org_row];
  elt->next.row = mtx->hdr.col[org_col];
  mtx->hdr.row[org_row] = mtx->hdr.col[org_col] = elt;
  return(elt);
}

struct element_t *mtx_create_element_value(mtx_matrix_t mtx,
                                           int32 org_row,
                                           int32 org_col,
                                           real64 val)
{
  register struct element_t *elt;

#if MTX_DEBUG
  if( NOTNULL(mtx_find_element(mtx,org_row,org_col)) ) {
    FPRINTF(g_mtxerr,"ERROR:  (mtx) create_element_value\n");
    FPRINTF(g_mtxerr,"        Element (%d,%d) already exists.\n",
            org_row,org_col);
  }
#endif
  elt = (struct element_t *)mem_get_element(mtx->ms);
  /* guess who didn't check for the alloc return. not needed here. */
  elt->value = val;
  elt->row = org_row;
  elt->col = org_col;
  elt->next.col = mtx->hdr.row[org_row];
  elt->next.row = mtx->hdr.col[org_col];
  mtx->hdr.row[org_row] = mtx->hdr.col[org_col] = elt;
  return(elt);
}

/* new version */
struct element_t *mtx_next_col(register struct element_t *elt,
                               mtx_range_t *rng, int32 *tocur)
{
   if( NOTNULL(elt) ) {
     elt = elt->next.col;
   } else {
     return elt;
   }
   /* use of not_in_range in the following is a loser on alphas */
   /* due to the price of assigning the range limits */
   if( rng != mtx_ALL_COLS ) {
      if( rng->high < rng->low ) return NULL;
      while( NOTNULL(elt) && !in_range(rng,tocur[elt->col]) )
         elt = elt->next.col;
   }
   return( elt );
}

/* original version */
struct element_t *mtx_next_row(register struct element_t *elt,
                               mtx_range_t *rng, int32 *tocur)
{
   if( NOTNULL(elt) ) elt = elt->next.row;
   if( rng != mtx_ALL_ROWS ) {
      if( rng->high < rng->low ) return NULL;
      while( NOTNULL(elt) && !in_range(rng,tocur[elt->row]) )
         elt = elt->next.row;
   }
   return( elt );
}


int32 *mtx_alloc_perm(int32 cap)
{
  int32 *perm;
  perm = (int32 *)ascmalloc( (cap+1)*sizeof(int32) );
  *perm = -1;
  return( perm+1 );
}

void mtx_copy_perm( int32 *tarperm, int32 *srcperm, int32 cap)
/**
 ***  Copies srcperm to tarperm given the capacity of srcperm.
 ***  If tarperm was obtained from
 ***  alloc_perm(), the -1 has already been copied
 **/
{
  mem_copy_cast(srcperm,tarperm,cap*sizeof(int32));
}

#define free_perm(perm) ascfree( (POINTER)((perm)-1) )
/**
 ***  Frees a permutation vector.
 **/

void mtx_free_perm(int32 *perm)
{
  free_perm(perm);
}

/* some local scope globals to keep memory so we aren't constantly
   reallocating */
struct reusable_data_vector
  g_mtx_null_index_data = {NULL,0,sizeof(int32),0},
  g_mtx_null_sum_data = {NULL,0,sizeof(real64),0},
  g_mtx_null_mark_data = {NULL,0,sizeof(char),0},
  g_mtx_null_vector_data = {NULL,0,sizeof(struct element_t *),0},
  g_mtx_null_col_vector_data = {NULL,0,sizeof(struct element_t *),0},
  g_mtx_null_row_vector_data = {NULL,0,sizeof(struct element_t *),0};

/* the psychologist function. resets data to zero after some checking */
static void mtx_sanify_reuseable(struct reusable_data_vector *r)
{
  if (r->last_line) {
    r->last_line = ZERO;
    if (r->capacity > 0 && r->arr != NULL) {
      /* rezero the buffer */
      switch (r->entry_size) {
      case sizeof(int32):
        mtx_zero_int32(r->arr,r->capacity);
        break;
      case sizeof(real64):
        mtx_zero_real64(r->arr,r->capacity);
        break;
      default:
        mtx_zero_char(r->arr, r->capacity*r->entry_size);
        break;
      }
    }
  }
  if (r->capacity && r->arr == NULL) {
    /* forget it exists (capacity && !arr --> insanity) */
    r->capacity = 0;
    r->arr = NULL;
  }
}

/* a function to insure the sanity of all the reuseable vectors */
void mtx_reset_null_vectors(void) {
  mtx_sanify_reuseable( &g_mtx_null_index_data);
  mtx_sanify_reuseable( &g_mtx_null_sum_data);
  mtx_sanify_reuseable( &g_mtx_null_mark_data);
  mtx_sanify_reuseable( &g_mtx_null_row_vector_data);
  mtx_sanify_reuseable( &g_mtx_null_col_vector_data);
  mtx_sanify_reuseable( &g_mtx_null_vector_data);
}

/* get a chunk of memory (probably already allocated) and return
   to user who will reset it to 0s before releasing it. */
void *mtx_null_vector_f(int32 cap, int line, CONST char *file,
                        struct reusable_data_vector *ptr, char *fn)
/**
 ***  Call only via the macros.
 ***  Bugs: does not check for null return.
 ***  but then if this returns null, a crash from some other cause
 ***  is imminent anyway.
 **/
{
  if (ptr->last_line) {
    FPRINTF(g_mtxerr,"Warning: (%s) mtx_%s called while data in use\n",file,fn);
    FPRINTF(g_mtxerr,"         Last called line: %d, this call line %d\n",
            ptr->last_line,line);
    if (cap > 0 && cap <= ptr->capacity) {
      switch (ptr->entry_size) {
      case sizeof(int32):
    mtx_zero_int32(ptr->arr,ptr->capacity);
    break;
      case sizeof(real64):
    mtx_zero_real64(ptr->arr,ptr->capacity);
    break;
      default:
        mtx_zero_char(ptr->arr, ptr->capacity*ptr->entry_size);
        break;
      }
    }
  }
  if (!cap) {
    free_unless_null( (POINTER)ptr->arr );
    ptr->capacity = ZERO;
    ptr->arr = NULL;
    ptr->last_line = ZERO;
    return( ptr->arr );
  }

  if( cap > ptr->capacity ) {
    free_unless_null( (POINTER)ptr->arr );
    ptr->arr = asccalloc( cap, ptr->entry_size );
    if (ISNULL(ptr->arr)) {
      ptr->capacity = ZERO;
      ptr->last_line = ZERO;
      FPRINTF(g_mtxerr,"ERROR: memory allocation failed in mtx_%s.\n",fn);
      return ptr->arr;
    }
    ptr->capacity = cap;
  }
  ptr->last_line = line;
  return( ptr->arr );
}

/**
 ***  Call only via the macros.
 **/
void mtx_null_vector_release_f(int line, CONST char *file,
                               struct reusable_data_vector *ptr, char *fn)
{
  if (ptr->last_line == 0) {
    FPRINTF(g_mtxerr,
      "Warning:  (%s) mtx_%s_release called (line %d) on vector not in use\n",
       file, fn, line);
  }
  ptr->last_line = 0;
}

/********************************************************************/

struct element_t **mtx_expand_row(mtx_matrix_t mtx, int32 org)
{
   struct element_t **arr;
   struct element_t *elt;

   arr = mtx_null_vector(mtx->order);
   for( elt=mtx->hdr.row[org]; NOTNULL(elt); elt=elt->next.col )
      arr[elt->col] = elt;
   return(arr);
}

struct element_t **mtx_expand_col(mtx_matrix_t mtx, int32 org)
{
   struct element_t **arr;
   struct element_t *elt;

   arr = mtx_null_vector(mtx->order);
   for( elt = mtx->hdr.col[org] ; NOTNULL(elt) ; elt = elt->next.row )
      arr[elt->row] = elt;
   return(arr);
}

void mtx_renull_using_row(mtx_matrix_t mtx, int32 org,
                          struct element_t **arr)
{
   struct element_t *elt;

   for( elt=mtx->hdr.row[org]; NOTNULL(elt); elt=elt->next.col )
      arr[elt->col] = NULL;
}

void mtx_renull_using_col(mtx_matrix_t mtx, int32 org,
                          struct element_t **arr)
{
   struct element_t *elt;

   for( elt=mtx->hdr.col[org]; NOTNULL(elt); elt=elt->next.row )
      arr[elt->row] = NULL;
}

void mtx_renull_all(mtx_matrix_t mtx, struct element_t **arr)
/**
 ***  Makes arr NULLed again, assuming it is size mtx->order.
 **/
{
  if (mtx && arr)
    zero(arr,mtx->order,struct element_t *);
}

#undef __MTX_C_SEEN__
1	aw0a	1	/*
2			* mtx: Ascend Sparse Matrix Package
3			* by Karl Michael Westerberg
4			* Created: 2/6/90
5			* Version: $Revision: 1.9 $
6			* Version control file: $RCSfile: mtx_use_only.c,v $
7			* Date last modified: $Date: 2000/01/25 02:27:12 $
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			#include "utilities/ascConfig.h"
34			#include "utilities/ascMalloc.h"
35			#include "utilities/mem.h"
36			#include "solver/mtx.h"
37			/* grab our private parts */
38			#define __MTX_C_SEEN__
39			#include "solver/mtx_use_only.h"
40
41			#define free_unless_null(ptr) if( NOTNULL(ptr) ) ascfree(ptr)
42
43			FILE *g_mtxerr;
44
45			struct element_t *mtx_find_element(mtx_matrix_t mtx, int32 org_row,
46			int32 org_col)
47			{
48			register struct element_t *elt;
49
50			/* old way.
51			for( elt = mtx->hdr.row[org_row] ; NOTNULL(elt) ; elt = elt->next.col )
52			if( elt->col == org_col )
53			break;
54			*/
55			for( elt = mtx->hdr.row[org_row] ;
56			NOTNULL(elt) && elt->col!=org_col ;
57			elt = elt->next.col );
58			return(elt);
59			}
60
61			struct element_t *mtx_create_element(mtx_matrix_t mtx,
62			int32 org_row, int32 org_col)
63			/**
64			*** Creates the given element and returns a pointer to it. The value is
65			*** initially zero.
66			**/
67			{
68			register struct element_t *elt;
69
70			#if MTX_DEBUG
71			if( NOTNULL(mtx_find_element(mtx,org_row,org_col)) ) {
72			FPRINTF(g_mtxerr,"ERROR: (mtx) create_element\n");
73			FPRINTF(g_mtxerr," Element (%d,%d) already exists.\n",
74			org_row,org_col);
75			}
76			#endif
77			elt = (struct element_t *)mem_get_element(mtx->ms);
78			/* guess who didn't check for the alloc return. */
79			elt->value = 0.0;
80			elt->row = org_row;
81			elt->col = org_col;
82			elt->next.col = mtx->hdr.row[org_row];
83			elt->next.row = mtx->hdr.col[org_col];
84			mtx->hdr.row[org_row] = mtx->hdr.col[org_col] = elt;
85			return(elt);
86			}
87
88			struct element_t *mtx_create_element_value(mtx_matrix_t mtx,
89			int32 org_row,
90			int32 org_col,
91			real64 val)
92			{
93			register struct element_t *elt;
94
95			#if MTX_DEBUG
96			if( NOTNULL(mtx_find_element(mtx,org_row,org_col)) ) {
97			FPRINTF(g_mtxerr,"ERROR: (mtx) create_element_value\n");
98			FPRINTF(g_mtxerr," Element (%d,%d) already exists.\n",
99			org_row,org_col);
100			}
101			#endif
102			elt = (struct element_t *)mem_get_element(mtx->ms);
103			/* guess who didn't check for the alloc return. not needed here. */
104			elt->value = val;
105			elt->row = org_row;
106			elt->col = org_col;
107			elt->next.col = mtx->hdr.row[org_row];
108			elt->next.row = mtx->hdr.col[org_col];
109			mtx->hdr.row[org_row] = mtx->hdr.col[org_col] = elt;
110			return(elt);
111			}
112
113			/* new version */
114			struct element_t mtx_next_col(register struct element_t elt,
115			mtx_range_t rng, int32 tocur)
116			{
117			if( NOTNULL(elt) ) {
118			elt = elt->next.col;
119			} else {
120			return elt;
121			}
122			/* use of not_in_range in the following is a loser on alphas */
123			/* due to the price of assigning the range limits */
124			if( rng != mtx_ALL_COLS ) {
125			if( rng->high < rng->low ) return NULL;
126			while( NOTNULL(elt) && !in_range(rng,tocur[elt->col]) )
127			elt = elt->next.col;
128			}
129			return( elt );
130			}
131
132			/* original version */
133			struct element_t mtx_next_row(register struct element_t elt,
134			mtx_range_t rng, int32 tocur)
135			{
136			if( NOTNULL(elt) ) elt = elt->next.row;
137			if( rng != mtx_ALL_ROWS ) {
138			if( rng->high < rng->low ) return NULL;
139			while( NOTNULL(elt) && !in_range(rng,tocur[elt->row]) )
140			elt = elt->next.row;
141			}
142			return( elt );
143			}
144
145
146			int32 *mtx_alloc_perm(int32 cap)
147			{
148			int32 *perm;
149			perm = (int32 )ascmalloc( (cap+1)sizeof(int32) );
150			*perm = -1;
151			return( perm+1 );
152			}
153
154			void mtx_copy_perm( int32 tarperm, int32 srcperm, int32 cap)
155			/**
156			*** Copies srcperm to tarperm given the capacity of srcperm.
157			*** If tarperm was obtained from
158			*** alloc_perm(), the -1 has already been copied
159			**/
160			{
161			mem_copy_cast(srcperm,tarperm,cap*sizeof(int32));
162			}
163
164			#define free_perm(perm) ascfree( (POINTER)((perm)-1) )
165			/**
166			*** Frees a permutation vector.
167			**/
168
169			void mtx_free_perm(int32 *perm)
170			{
171			free_perm(perm);
172			}
173
174			/* some local scope globals to keep memory so we aren't constantly
175			reallocating */
176			struct reusable_data_vector
177			g_mtx_null_index_data = {NULL,0,sizeof(int32),0},
178			g_mtx_null_sum_data = {NULL,0,sizeof(real64),0},
179			g_mtx_null_mark_data = {NULL,0,sizeof(char),0},
180			g_mtx_null_vector_data = {NULL,0,sizeof(struct element_t *),0},
181			g_mtx_null_col_vector_data = {NULL,0,sizeof(struct element_t *),0},
182			g_mtx_null_row_vector_data = {NULL,0,sizeof(struct element_t *),0};
183
184			/* the psychologist function. resets data to zero after some checking */
185			static void mtx_sanify_reuseable(struct reusable_data_vector *r)
186			{
187			if (r->last_line) {
188			r->last_line = ZERO;
189			if (r->capacity > 0 && r->arr != NULL) {
190			/* rezero the buffer */
191			switch (r->entry_size) {
192			case sizeof(int32):
193			mtx_zero_int32(r->arr,r->capacity);
194			break;
195			case sizeof(real64):
196			mtx_zero_real64(r->arr,r->capacity);
197			break;
198			default:
199			mtx_zero_char(r->arr, r->capacity*r->entry_size);
200			break;
201			}
202			}
203			}
204			if (r->capacity && r->arr == NULL) {
205			/* forget it exists (capacity && !arr --> insanity) */
206			r->capacity = 0;
207			r->arr = NULL;
208			}
209			}
210
211			/* a function to insure the sanity of all the reuseable vectors */
212			void mtx_reset_null_vectors(void) {
213			mtx_sanify_reuseable( &g_mtx_null_index_data);
214			mtx_sanify_reuseable( &g_mtx_null_sum_data);
215			mtx_sanify_reuseable( &g_mtx_null_mark_data);
216			mtx_sanify_reuseable( &g_mtx_null_row_vector_data);
217			mtx_sanify_reuseable( &g_mtx_null_col_vector_data);
218			mtx_sanify_reuseable( &g_mtx_null_vector_data);
219			}
220
221			/* get a chunk of memory (probably already allocated) and return
222			to user who will reset it to 0s before releasing it. */
223			void mtx_null_vector_f(int32 cap, int line, CONST char file,
224			struct reusable_data_vector ptr, char fn)
225			/**
226			*** Call only via the macros.
227			*** Bugs: does not check for null return.
228			*** but then if this returns null, a crash from some other cause
229			*** is imminent anyway.
230			**/
231			{
232			if (ptr->last_line) {
233			FPRINTF(g_mtxerr,"Warning: (%s) mtx_%s called while data in use\n",file,fn);
234			FPRINTF(g_mtxerr," Last called line: %d, this call line %d\n",
235			ptr->last_line,line);
236			if (cap > 0 && cap <= ptr->capacity) {
237			switch (ptr->entry_size) {
238			case sizeof(int32):
239			mtx_zero_int32(ptr->arr,ptr->capacity);
240			break;
241			case sizeof(real64):
242			mtx_zero_real64(ptr->arr,ptr->capacity);
243			break;
244			default:
245			mtx_zero_char(ptr->arr, ptr->capacity*ptr->entry_size);
246			break;
247			}
248			}
249			}
250			if (!cap) {
251			free_unless_null( (POINTER)ptr->arr );
252			ptr->capacity = ZERO;
253			ptr->arr = NULL;
254			ptr->last_line = ZERO;
255			return( ptr->arr );
256			}
257
258			if( cap > ptr->capacity ) {
259			free_unless_null( (POINTER)ptr->arr );
260			ptr->arr = asccalloc( cap, ptr->entry_size );
261			if (ISNULL(ptr->arr)) {
262			ptr->capacity = ZERO;
263			ptr->last_line = ZERO;
264			FPRINTF(g_mtxerr,"ERROR: memory allocation failed in mtx_%s.\n",fn);
265			return ptr->arr;
266			}
267			ptr->capacity = cap;
268			}
269			ptr->last_line = line;
270			return( ptr->arr );
271			}
272
273			/**
274			*** Call only via the macros.
275			**/
276			void mtx_null_vector_release_f(int line, CONST char *file,
277			struct reusable_data_vector ptr, char fn)
278			{
279			if (ptr->last_line == 0) {
280			FPRINTF(g_mtxerr,
281			"Warning: (%s) mtx_%s_release called (line %d) on vector not in use\n",
282			file, fn, line);
283			}
284			ptr->last_line = 0;
285			}
286
287			/********************************************************************/
288
289			struct element_t **mtx_expand_row(mtx_matrix_t mtx, int32 org)
290			{
291			struct element_t **arr;
292			struct element_t *elt;
293
294			arr = mtx_null_vector(mtx->order);
295			for( elt=mtx->hdr.row[org]; NOTNULL(elt); elt=elt->next.col )
296			arr[elt->col] = elt;
297			return(arr);
298			}
299
300			struct element_t **mtx_expand_col(mtx_matrix_t mtx, int32 org)
301			{
302			struct element_t **arr;
303			struct element_t *elt;
304
305			arr = mtx_null_vector(mtx->order);
306			for( elt = mtx->hdr.col[org] ; NOTNULL(elt) ; elt = elt->next.row )
307			arr[elt->row] = elt;
308			return(arr);
309			}
310
311			void mtx_renull_using_row(mtx_matrix_t mtx, int32 org,
312			struct element_t **arr)
313			{
314			struct element_t *elt;
315
316			for( elt=mtx->hdr.row[org]; NOTNULL(elt); elt=elt->next.col )
317			arr[elt->col] = NULL;
318			}
319
320			void mtx_renull_using_col(mtx_matrix_t mtx, int32 org,
321			struct element_t **arr)
322			{
323			struct element_t *elt;
324
325			for( elt=mtx->hdr.col[org]; NOTNULL(elt); elt=elt->next.row )
326			arr[elt->row] = NULL;
327			}
328
329			void mtx_renull_all(mtx_matrix_t mtx, struct element_t **arr)
330			/**
331			*** Makes arr NULLed again, assuming it is size mtx->order.
332			**/
333			{
334			if (mtx && arr)
335			zero(arr,mtx->order,struct element_t *);
336			}
337
338			#undef __MTX_C_SEEN__