| 1 |
aw0a |
1 |
/* |
| 2 |
|
|
* General routines for Slv interface. |
| 3 |
|
|
* by Ben Allan |
| 4 |
|
|
* Created 1/95 |
| 5 |
|
|
* Version: $Revision: 1.25 $ |
| 6 |
|
|
* Version control file: $RCSfile: slv_common.c,v $ |
| 7 |
|
|
* Date last modified: $Date: 1998/02/02 23:16:22 $ |
| 8 |
|
|
* Last modified by: $Author: ballan $ |
| 9 |
|
|
* This file is part of the SLV solver. |
| 10 |
|
|
* |
| 11 |
|
|
* Copyright (C) 1990 Karl Michael Westerberg |
| 12 |
|
|
* Copyright (C) 1993 Joseph Zaher |
| 13 |
|
|
* Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan |
| 14 |
|
|
* |
| 15 |
|
|
* The SLV solver is free software; you can redistribute |
| 16 |
|
|
* it and/or modify it under the terms of the GNU General Public License as |
| 17 |
|
|
* published by the Free Software Foundation; either version 2 of the |
| 18 |
|
|
* License, or (at your option) any later version. |
| 19 |
|
|
* |
| 20 |
|
|
* The SLV solver is distributed in hope that it will be |
| 21 |
|
|
* useful, but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 22 |
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 23 |
|
|
* General Public License for more details. |
| 24 |
|
|
* |
| 25 |
|
|
* You should have received a copy of the GNU General Public License |
| 26 |
|
|
* along with the program; if not, write to the Free Software Foundation, |
| 27 |
|
|
* Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named |
| 28 |
|
|
* COPYING. COPYING is found in ../compiler. |
| 29 |
|
|
*/ |
| 30 |
|
|
|
| 31 |
|
|
/* |
| 32 |
|
|
* general C utility routines for slv class interfaces. Abstracted from |
| 33 |
|
|
* slvX.c January 1995. Ben Allan. |
| 34 |
|
|
*/ |
| 35 |
|
|
|
| 36 |
|
|
#include <math.h> |
| 37 |
|
|
#include "utilities/ascConfig.h" |
| 38 |
|
|
#include "utilities/ascSignal.h" |
| 39 |
|
|
#include "compiler/compiler.h" |
| 40 |
|
|
#include "utilities/ascMalloc.h" |
| 41 |
|
|
#include "solver/mtx.h" |
| 42 |
|
|
#include "solver/slv_types.h" |
| 43 |
|
|
#include "solver/rel.h" |
| 44 |
|
|
#include "solver/var.h" |
| 45 |
|
|
#include "solver/discrete.h" |
| 46 |
|
|
#include "solver/logrel.h" |
| 47 |
|
|
#include "solver/slv_common.h" |
| 48 |
|
|
#include "utilities/mem.h" |
| 49 |
|
|
/* if libasc.a running around, the following: */ |
| 50 |
|
|
#if SLV_INSTANCES |
| 51 |
|
|
#include "compiler/fractions.h" |
| 52 |
|
|
#include "compiler/dimen.h" |
| 53 |
|
|
#include "compiler/functype.h" |
| 54 |
|
|
#include "compiler/func.h" |
| 55 |
|
|
#include "solver/relman.h" |
| 56 |
|
|
#include "solver/logrelman.h" |
| 57 |
|
|
#else |
| 58 |
|
|
#ifdef NDEBUG |
| 59 |
|
|
#define ascnint(a) (((int) (a)>=0.0 ? floor((a) + 0.5) : -floor(0.5 - (a)))) |
| 60 |
|
|
#else |
| 61 |
|
|
#define ascnint(a) ascnintF(a) |
| 62 |
|
|
#endif /* NDEBUG */ |
| 63 |
|
|
#endif /* instances */ |
| 64 |
|
|
|
| 65 |
|
|
#define SAFE_FIX_ME 0 |
| 66 |
|
|
|
| 67 |
|
|
/** |
| 68 |
|
|
*** Array/vector operations |
| 69 |
|
|
*** ---------------------------- |
| 70 |
|
|
*** slv_zero_vector(vec) |
| 71 |
|
|
*** slv_copy_vector(vec1,vec2) |
| 72 |
|
|
*** prod = slv_inner_product(vec1,vec2) |
| 73 |
|
|
*** norm2 = slv_square_norm(vec) |
| 74 |
|
|
*** slv_matrix_product(mtx,vec,prod,scale,transpose) |
| 75 |
|
|
*** slv_write_vector(file,vec) |
| 76 |
|
|
**/ |
| 77 |
|
|
|
| 78 |
|
|
void slv_zero_vector( struct vector_data *vec) |
| 79 |
|
|
{ |
| 80 |
|
|
real64 *p; |
| 81 |
|
|
int32 len; |
| 82 |
|
|
|
| 83 |
|
|
p = vec->vec + vec->rng->low; |
| 84 |
|
|
len = vec->rng->high - vec->rng->low + 1; |
| 85 |
|
|
mtx_zero_real64(p,len); |
| 86 |
|
|
} |
| 87 |
|
|
|
| 88 |
|
|
void slv_copy_vector( struct vector_data *vec1,struct vector_data *vec2) |
| 89 |
|
|
{ |
| 90 |
|
|
real64 *p1,*p2; |
| 91 |
|
|
int32 len; |
| 92 |
|
|
|
| 93 |
|
|
p1 = vec1->vec + vec1->rng->low; |
| 94 |
|
|
p2 = vec2->vec + vec2->rng->low; |
| 95 |
|
|
len = vec1->rng->high - vec1->rng->low + 1; |
| 96 |
|
|
/*mem_copy_cast not in order here, probably */ |
| 97 |
|
|
mem_move_cast(p1,p2,len*sizeof(real64)); |
| 98 |
|
|
} |
| 99 |
|
|
|
| 100 |
|
|
#define USEDOT TRUE |
| 101 |
|
|
/* USEDOT = TRUE is a winner on alphas, hps, and sparc20 */ |
| 102 |
|
|
real64 slv_inner_product(struct vector_data *vec1 , |
| 103 |
|
|
struct vector_data *vec2) |
| 104 |
|
|
/** |
| 105 |
|
|
*** Computes inner product between vec1 and vec2, returning result. |
| 106 |
|
|
*** vec1 and vec2 may overlap or even be identical. |
| 107 |
|
|
**/ |
| 108 |
|
|
{ |
| 109 |
|
|
real64 *p1,*p2; |
| 110 |
|
|
#if !USEDOT |
| 111 |
|
|
real64 sum; |
| 112 |
|
|
#endif |
| 113 |
|
|
int32 len; |
| 114 |
|
|
|
| 115 |
|
|
p1 = vec1->vec + vec1->rng->low; |
| 116 |
|
|
p2 = vec2->vec + vec2->rng->low; |
| 117 |
|
|
len = vec1->rng->high - vec1->rng->low + 1; |
| 118 |
|
|
#if !USEDOT |
| 119 |
|
|
if (p1 != p2) { |
| 120 |
|
|
for( sum=0.0 ; --len >= 0 ; ++p1,++p2 ) { |
| 121 |
|
|
sum += (*p1) * (*p2); |
| 122 |
|
|
} |
| 123 |
|
|
return(sum); |
| 124 |
|
|
} else { |
| 125 |
|
|
for( sum=0.0 ; --len >= 0 ; ++p1 ) { |
| 126 |
|
|
sum += (*p1) * (*p1); |
| 127 |
|
|
} |
| 128 |
|
|
return(sum); |
| 129 |
|
|
} |
| 130 |
|
|
#else |
| 131 |
|
|
return slv_dot(len,p1,p2); |
| 132 |
|
|
#endif |
| 133 |
|
|
} |
| 134 |
|
|
|
| 135 |
|
|
real64 slv_square_norm(struct vector_data *vec) |
| 136 |
|
|
/** |
| 137 |
|
|
*** Computes norm^2 of vector, assigning the result to vec->norm2 |
| 138 |
|
|
*** and returning the result as well. |
| 139 |
|
|
**/ |
| 140 |
|
|
{ |
| 141 |
|
|
vec->norm2 = slv_inner_product(vec,vec); |
| 142 |
|
|
return vec->norm2; |
| 143 |
|
|
} |
| 144 |
|
|
|
| 145 |
|
|
void slv_matrix_product(mtx_matrix_t mtx, struct vector_data *vec, |
| 146 |
|
|
struct vector_data *prod, real64 scale, |
| 147 |
|
|
boolean transpose) |
| 148 |
|
|
/** |
| 149 |
|
|
*** Stores prod := (scale)*(mtx)(vec) or (scale)*(mtx-transpose)(vec). |
| 150 |
|
|
*** vec and prod must be completely different. |
| 151 |
|
|
**/ |
| 152 |
|
|
{ |
| 153 |
|
|
mtx_coord_t nz; |
| 154 |
|
|
real64 value, *vvec, *pvec; |
| 155 |
|
|
int32 lim; |
| 156 |
|
|
|
| 157 |
|
|
lim=prod->rng->high; |
| 158 |
|
|
pvec=prod->vec; |
| 159 |
|
|
vvec=vec->vec; |
| 160 |
|
|
if( transpose ) { |
| 161 |
|
|
for(nz.col = prod->rng->low ; nz.col <= lim ; ++(nz.col) ) { |
| 162 |
|
|
pvec[nz.col] = 0.0; |
| 163 |
|
|
nz.row = mtx_FIRST; |
| 164 |
|
|
while( value = mtx_next_in_col(mtx,&nz,vec->rng), |
| 165 |
|
|
nz.row != mtx_LAST ) |
| 166 |
|
|
pvec[nz.col] += value*vvec[nz.row]; |
| 167 |
|
|
pvec[nz.col] *= scale; |
| 168 |
|
|
} |
| 169 |
|
|
} else { |
| 170 |
|
|
for(nz.row = prod->rng->low ; nz.row <= lim ; ++(nz.row) ) { |
| 171 |
|
|
pvec[nz.row] = 0.0; |
| 172 |
|
|
nz.col = mtx_FIRST; |
| 173 |
|
|
while( value = mtx_next_in_row(mtx,&nz,vec->rng), |
| 174 |
|
|
nz.col != mtx_LAST ) |
| 175 |
|
|
pvec[nz.row] += value*vvec[nz.col]; |
| 176 |
|
|
pvec[nz.row] *= scale; |
| 177 |
|
|
} |
| 178 |
|
|
} |
| 179 |
|
|
} |
| 180 |
|
|
|
| 181 |
|
|
void slv_write_vector(FILE *fp,struct vector_data *vec) |
| 182 |
|
|
/** |
| 183 |
|
|
*** Outputs a vector. |
| 184 |
|
|
**/ |
| 185 |
|
|
{ |
| 186 |
|
|
int32 ndx,hi; |
| 187 |
|
|
real64 *vvec; |
| 188 |
|
|
vvec = vec->vec; |
| 189 |
|
|
hi = vec->rng->high; |
| 190 |
|
|
FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n", |
| 191 |
|
|
sqrt(fabs(vec->norm2)), vec->accurate?"TRUE":"FALSE", |
| 192 |
|
|
vec->rng->low,vec->rng->high); |
| 193 |
|
|
FPRINTF(fp,"Vector --> "); |
| 194 |
|
|
for( ndx=vec->rng->low ; ndx<=hi ; ++ndx ) |
| 195 |
|
|
FPRINTF(fp, "%g ", vvec[ndx]); |
| 196 |
|
|
PUTC('\n',fp); |
| 197 |
|
|
} |
| 198 |
|
|
|
| 199 |
|
|
/* Dot product for loop unrolled vector norms */ |
| 200 |
|
|
real64 slv_dot(int32 len, const real64 *p1, const real64 *p2) |
| 201 |
|
|
{ |
| 202 |
|
|
register double sum,lsum; |
| 203 |
|
|
int m,n; |
| 204 |
|
|
|
| 205 |
|
|
/* AVMAGIC in fact isn't magic. |
| 206 |
|
|
* only goes to 2-10. change the code below if you mess with it. |
| 207 |
|
|
* Default AVMAGIC is 4, which works well on alphas, tika. |
| 208 |
|
|
*/ |
| 209 |
|
|
#define AVMAGIC 4 |
| 210 |
|
|
#ifdef sun |
| 211 |
|
|
#undef AVMAGIC |
| 212 |
|
|
#define AVMAGIC 6 |
| 213 |
|
|
#endif |
| 214 |
|
|
#ifdef __hpux |
| 215 |
|
|
#undef AVMAGIC |
| 216 |
|
|
#define AVMAGIC 4 |
| 217 |
|
|
/* 2 was best value on ranier(9000/720) but tika (9000/715) likes 4 |
| 218 |
|
|
* there are no recognizable (defines) compiler differences, so tika |
| 219 |
|
|
* will set the hp default |
| 220 |
|
|
*/ |
| 221 |
|
|
#endif |
| 222 |
|
|
/* hands down best avmagic on atlas is 4, ranier 2, unxi21 6 |
| 223 |
|
|
* under native compilers. no bets on the billions of gcc variants. |
| 224 |
|
|
* needs to be tried on sgi. |
| 225 |
|
|
* Note, these are tuned for speed, not for accuracy. on very large |
| 226 |
|
|
* vectors something like dnrm2 may be more appropriate, though |
| 227 |
|
|
* it is very slow. upping avmagic to 10 will help accuracy some. |
| 228 |
|
|
*/ |
| 229 |
|
|
|
| 230 |
|
|
#if (AVMAGIC>10) |
| 231 |
|
|
#undef AVMAGIC |
| 232 |
|
|
#define AVMAGIC 10 |
| 233 |
|
|
#endif |
| 234 |
|
|
m = len / AVMAGIC; |
| 235 |
|
|
n = len % AVMAGIC; |
| 236 |
|
|
if (p1!=p2) { |
| 237 |
|
|
/* get leading leftovers */ |
| 238 |
|
|
for( sum=0.0 ; --n >= 0 ; ) { |
| 239 |
|
|
sum += (*p1) * (*p2); |
| 240 |
|
|
++p1; ++p2; |
| 241 |
|
|
} |
| 242 |
|
|
/* p1,p2 now point at first unadded location, or just after the end (m=0)*/ |
| 243 |
|
|
/* eat m chunks */ |
| 244 |
|
|
for( n=0; n <m; n++) { |
| 245 |
|
|
/* now, as i am too lazy to figure out a macro that expands itself */ |
| 246 |
|
|
lsum = (*p1) * (*p2); /* zeroth term is assigned */ |
| 247 |
|
|
p1++; p2++; |
| 248 |
|
|
lsum += (*p1) * (*p2); /* add 1th */ |
| 249 |
|
|
#if (AVMAGIC>2) |
| 250 |
|
|
p1++; p2++; |
| 251 |
|
|
lsum += (*p1) * (*p2); /* add 2th */ |
| 252 |
|
|
#if (AVMAGIC>3) |
| 253 |
|
|
p1++; p2++; |
| 254 |
|
|
lsum += (*p1) * (*p2); /* add 3th */ |
| 255 |
|
|
#if (AVMAGIC>4) |
| 256 |
|
|
p1++; p2++; |
| 257 |
|
|
lsum += (*p1) * (*p2); /* add 4th */ |
| 258 |
|
|
#if (AVMAGIC>5) |
| 259 |
|
|
p1++; p2++; |
| 260 |
|
|
lsum += (*p1) * (*p2); /* add 5th */ |
| 261 |
|
|
#if (AVMAGIC>6) |
| 262 |
|
|
p1++; p2++; |
| 263 |
|
|
lsum += (*p1) * (*p2); /* add 6th */ |
| 264 |
|
|
#if (AVMAGIC>7) |
| 265 |
|
|
p1++; p2++; |
| 266 |
|
|
lsum += (*p1) * (*p2); /* add 7th */ |
| 267 |
|
|
#if (AVMAGIC>8) |
| 268 |
|
|
p1++; p2++; |
| 269 |
|
|
lsum += (*p1) * (*p2); /* add 8th */ |
| 270 |
|
|
#if (AVMAGIC>9) |
| 271 |
|
|
p1++; p2++; |
| 272 |
|
|
lsum += (*p1) * (*p2); /* add 9th */ |
| 273 |
|
|
#endif |
| 274 |
|
|
#endif |
| 275 |
|
|
#endif |
| 276 |
|
|
#endif |
| 277 |
|
|
#endif |
| 278 |
|
|
#endif |
| 279 |
|
|
#endif |
| 280 |
|
|
#endif |
| 281 |
|
|
p1++; p2++; /* leave p1,p2 pointing at next zeroth */ |
| 282 |
|
|
sum += lsum; |
| 283 |
|
|
} |
| 284 |
|
|
} else { |
| 285 |
|
|
/* get leading leftovers */ |
| 286 |
|
|
for( sum=0.0 ; --n >= 0 ; ++p1 ) { |
| 287 |
|
|
sum += (*p1) * (*p1); |
| 288 |
|
|
} |
| 289 |
|
|
/* p1 now points at first unadded location, or just after the end (m=0)*/ |
| 290 |
|
|
/* eat m chunks */ |
| 291 |
|
|
for( n=0; n <m; n++) { |
| 292 |
|
|
/* now, as i am too lazy to figure out a macro that expands itself */ |
| 293 |
|
|
lsum = (*p1) * (*p1); /* zeroth term is assigned */ |
| 294 |
|
|
p1++; |
| 295 |
|
|
lsum += (*p1) * (*p1); /* add 1th */ |
| 296 |
|
|
#if (AVMAGIC>2) |
| 297 |
|
|
p1++; |
| 298 |
|
|
lsum += (*p1) * (*p1); /* add 2th */ |
| 299 |
|
|
#if (AVMAGIC>3) |
| 300 |
|
|
p1++; |
| 301 |
|
|
lsum += (*p1) * (*p1); /* add 3th */ |
| 302 |
|
|
#if (AVMAGIC>4) |
| 303 |
|
|
p1++; |
| 304 |
|
|
lsum += (*p1) * (*p1); /* add 4th */ |
| 305 |
|
|
#if (AVMAGIC>5) |
| 306 |
|
|
p1++; |
| 307 |
|
|
lsum += (*p1) * (*p1); /* add 5th */ |
| 308 |
|
|
#if (AVMAGIC>6) |
| 309 |
|
|
p1++; |
| 310 |
|
|
lsum += (*p1) * (*p1); /* add 6th */ |
| 311 |
|
|
#if (AVMAGIC>7) |
| 312 |
|
|
p1++; |
| 313 |
|
|
lsum += (*p1) * (*p1); /* add 7th */ |
| 314 |
|
|
#if (AVMAGIC>8) |
| 315 |
|
|
p1++; |
| 316 |
|
|
lsum += (*p1) * (*p1); /* add 8th */ |
| 317 |
|
|
#if (AVMAGIC>9) |
| 318 |
|
|
p1++; |
| 319 |
|
|
lsum += (*p1) * (*p1); /* add 9th */ |
| 320 |
|
|
#endif |
| 321 |
|
|
#endif |
| 322 |
|
|
#endif |
| 323 |
|
|
#endif |
| 324 |
|
|
#endif |
| 325 |
|
|
#endif |
| 326 |
|
|
#endif |
| 327 |
|
|
#endif |
| 328 |
|
|
p1++; /* leave p1 pointing at next zeroth */ |
| 329 |
|
|
sum += lsum; |
| 330 |
|
|
} |
| 331 |
|
|
} |
| 332 |
|
|
return(sum); |
| 333 |
|
|
#undef AVMAGIC |
| 334 |
|
|
} |
| 335 |
|
|
|
| 336 |
|
|
/** |
| 337 |
|
|
*** General input/output routines |
| 338 |
|
|
*** ----------------------------- |
| 339 |
|
|
*** fp = MIF(sys) |
| 340 |
|
|
*** fp = LIF(sys) |
| 341 |
|
|
*** slv_print_var_name(out,sys,var) |
| 342 |
|
|
*** slv_print_rel_name(out,sys,rel) |
| 343 |
|
|
*** slv_print_dis_name(out,sys,dvar) |
| 344 |
|
|
*** slv_print_logrel_name(out,sys,lrel) |
| 345 |
|
|
*** NOt yet implemented correctly: |
| 346 |
|
|
*** slv_print_obj_name(out,obj) |
| 347 |
|
|
*** slv_print_var_sindex(out,var) |
| 348 |
|
|
*** slv_print_rel_sindex(out,rel) |
| 349 |
|
|
*** slv_print_dis_sindex(out,dvar) |
| 350 |
|
|
*** slv_print_logrel_sindex(out,lrel) |
| 351 |
|
|
*** slv_print_obj_index(out,obj) |
| 352 |
|
|
**/ |
| 353 |
|
|
|
| 354 |
|
|
FILE *slv_get_output_file(fp) |
| 355 |
|
|
FILE *fp; |
| 356 |
|
|
/** |
| 357 |
|
|
*** Returns fp if fp!=NULL, or a file pointer |
| 358 |
|
|
*** open to nul device if fp == NULL. |
| 359 |
|
|
**/ |
| 360 |
|
|
{ |
| 361 |
|
|
static FILE *nuldev = NULL; |
| 362 |
|
|
static char fname[] = "/dev/null"; |
| 363 |
|
|
|
| 364 |
|
|
if( fp==NULL ) { |
| 365 |
|
|
if(nuldev==NULL) |
| 366 |
|
|
if( (nuldev=fopen(fname,"w")) == NULL ) { |
| 367 |
|
|
FPRINTF(stderr,"ERROR: (slv) get_output_file\n"); |
| 368 |
|
|
FPRINTF(stderr," Unable to open %s.\n",fname); |
| 369 |
|
|
} |
| 370 |
|
|
fp=nuldev; |
| 371 |
|
|
} |
| 372 |
|
|
return(fp); |
| 373 |
|
|
} |
| 374 |
|
|
#define MIF(sys) slv_get_output_file( (sys)->p.output.more_important ) |
| 375 |
|
|
#define LIF(sys) slv_get_output_file( (sys)->p.output.less_important ) |
| 376 |
|
|
|
| 377 |
|
|
#if SLV_INSTANCES |
| 378 |
|
|
|
| 379 |
|
|
void slv_print_var_name( FILE *out,slv_system_t sys, struct var_variable *var) |
| 380 |
|
|
{ |
| 381 |
|
|
char *name=NULL; |
| 382 |
|
|
if (out == NULL || sys == NULL || var == NULL) return; |
| 383 |
|
|
name = var_make_name(sys,var); |
| 384 |
|
|
if( *name == '?' ) FPRINTF(out,"%d",var_sindex(var)); |
| 385 |
|
|
else FPRINTF(out,"%s",name); |
| 386 |
|
|
if (name) ascfree(name); |
| 387 |
|
|
} |
| 388 |
|
|
|
| 389 |
|
|
void slv_print_rel_name( FILE *out, slv_system_t sys, struct rel_relation *rel) |
| 390 |
|
|
{ |
| 391 |
|
|
char *name; |
| 392 |
|
|
name = rel_make_name(sys,rel); |
| 393 |
|
|
if( *name == '?' ) FPRINTF(out,"%d",rel_sindex(rel)); |
| 394 |
|
|
else FPRINTF(out,"%s",name); |
| 395 |
|
|
ascfree(name); |
| 396 |
|
|
} |
| 397 |
|
|
|
| 398 |
|
|
void slv_print_dis_name( FILE *out,slv_system_t sys, struct dis_discrete *dvar) |
| 399 |
|
|
{ |
| 400 |
|
|
char *name=NULL; |
| 401 |
|
|
if (out == NULL || sys == NULL || dvar == NULL) return; |
| 402 |
|
|
name = dis_make_name(sys,dvar); |
| 403 |
|
|
if( *name == '?' ) FPRINTF(out,"%d",dis_sindex(dvar)); |
| 404 |
|
|
else FPRINTF(out,"%s",name); |
| 405 |
|
|
if (name) ascfree(name); |
| 406 |
|
|
} |
| 407 |
|
|
|
| 408 |
|
|
void slv_print_logrel_name( FILE *out, slv_system_t sys, |
| 409 |
|
|
struct logrel_relation *lrel) |
| 410 |
|
|
{ |
| 411 |
|
|
char *name; |
| 412 |
|
|
name = logrel_make_name(sys,lrel); |
| 413 |
|
|
if( *name == '?' ) FPRINTF(out,"%d",logrel_sindex(lrel)); |
| 414 |
|
|
else FPRINTF(out,"%s",name); |
| 415 |
|
|
ascfree(name); |
| 416 |
|
|
} |
| 417 |
|
|
|
| 418 |
|
|
#ifdef NEWSTUFF |
| 419 |
|
|
void slv_print_obj_name(out,obj) |
| 420 |
|
|
FILE *out; |
| 421 |
|
|
obj_objective_t obj; |
| 422 |
|
|
{ |
| 423 |
|
|
char *name; |
| 424 |
|
|
name = obj_make_name(obj); |
| 425 |
|
|
if( *name == '?' ) FPRINTF(out,"%d",obj_index(obj)); |
| 426 |
|
|
else FPRINTF(out,"%s",name); |
| 427 |
|
|
ascfree(name); |
| 428 |
|
|
} |
| 429 |
|
|
#endif |
| 430 |
|
|
|
| 431 |
|
|
void slv_print_var_sindex(out,var) |
| 432 |
|
|
FILE *out; |
| 433 |
|
|
struct var_variable *var; |
| 434 |
|
|
{ |
| 435 |
|
|
FPRINTF(out,"%d",var_sindex(var)); |
| 436 |
|
|
} |
| 437 |
|
|
|
| 438 |
|
|
void slv_print_rel_sindex(out,rel) |
| 439 |
|
|
FILE *out; |
| 440 |
|
|
struct rel_relation *rel; |
| 441 |
|
|
{ |
| 442 |
|
|
FPRINTF(out,"%d",rel_sindex(rel)); |
| 443 |
|
|
} |
| 444 |
|
|
|
| 445 |
|
|
void slv_print_dis_sindex(out,dvar) |
| 446 |
|
|
FILE *out; |
| 447 |
|
|
struct dis_discrete *dvar; |
| 448 |
|
|
{ |
| 449 |
|
|
FPRINTF(out,"%d",dis_sindex(dvar)); |
| 450 |
|
|
} |
| 451 |
|
|
|
| 452 |
|
|
void slv_print_logrel_sindex(out,lrel) |
| 453 |
|
|
FILE *out; |
| 454 |
|
|
struct logrel_relation *lrel; |
| 455 |
|
|
{ |
| 456 |
|
|
FPRINTF(out,"%d",logrel_sindex(lrel)); |
| 457 |
|
|
} |
| 458 |
|
|
|
| 459 |
|
|
#ifdef NEWSTUFF |
| 460 |
|
|
void slv_print_obj_index(out,obj) |
| 461 |
|
|
FILE *out; |
| 462 |
|
|
struct rel_relation *rel; |
| 463 |
|
|
{ |
| 464 |
|
|
FPRINTF(out,"%d",obj_index(obj)); |
| 465 |
|
|
} |
| 466 |
|
|
#endif |
| 467 |
|
|
|
| 468 |
|
|
#define destroy_array(p) \ |
| 469 |
|
|
if( (p) != NULL ) ascfree((p)) |
| 470 |
|
|
|
| 471 |
|
|
int slv_direct_solve(slv_system_t server, struct rel_relation *rel, |
| 472 |
|
|
struct var_variable *var, FILE *fp, |
| 473 |
|
|
real64 epsilon, int ignore_bounds, int scaled) |
| 474 |
|
|
/** |
| 475 |
|
|
*** Attempt to directly solve the given relation (equality constraint) for |
| 476 |
|
|
*** the given variable, leaving the others fixed. Returns an integer |
| 477 |
|
|
*** signifying the status as one of the following three: |
| 478 |
|
|
*** |
| 479 |
|
|
*** 0 ==> Unable to determine anything, possibly due to fp exception. |
| 480 |
|
|
*** 1 ==> Solution found, variable value set to this solution. |
| 481 |
|
|
*** -1 ==> No solution found. |
| 482 |
|
|
*** |
| 483 |
|
|
*** The variable bounds will be upheld, unless they are to be ignored. |
| 484 |
|
|
**/ |
| 485 |
|
|
{ |
| 486 |
|
|
int32 able,status; |
| 487 |
|
|
int nsolns,allsolns; |
| 488 |
|
|
real64 *slist,save; |
| 489 |
|
|
|
| 490 |
|
|
slist = relman_directly_solve_new(rel,var,&able,&nsolns,epsilon); |
| 491 |
|
|
if( !able ) { |
| 492 |
|
|
return(0); |
| 493 |
|
|
} |
| 494 |
|
|
allsolns = nsolns; |
| 495 |
|
|
while( --nsolns >= 0 ) { |
| 496 |
|
|
if( ignore_bounds ) { |
| 497 |
|
|
var_set_value(var,slist[nsolns]); |
| 498 |
|
|
break; |
| 499 |
|
|
} |
| 500 |
|
|
if( var_lower_bound(var) > slist[nsolns] ) { |
| 501 |
|
|
save = var_value(var); |
| 502 |
|
|
var_set_value(var,var_lower_bound(var)); |
| 503 |
|
|
Asc_SignalHandlerPush(SIGFPE,SIG_IGN); |
| 504 |
|
|
relman_eval(rel,&status,SAFE_FIX_ME); |
| 505 |
|
|
Asc_SignalHandlerPop(SIGFPE,SIG_IGN); |
| 506 |
|
|
if (scaled) { |
| 507 |
|
|
if( relman_calc_satisfied_scaled(rel,epsilon) ) break; |
| 508 |
|
|
} else { |
| 509 |
|
|
if( relman_calc_satisfied(rel,epsilon) ) break; |
| 510 |
|
|
} |
| 511 |
|
|
var_set_value(var,save); |
| 512 |
|
|
} else if( var_upper_bound(var) < slist[nsolns] ) { |
| 513 |
|
|
save = var_value(var); |
| 514 |
|
|
var_set_value(var,var_upper_bound(var)); |
| 515 |
|
|
Asc_SignalHandlerPush(SIGFPE,SIG_IGN); |
| 516 |
|
|
relman_eval(rel,&status,SAFE_FIX_ME); |
| 517 |
|
|
Asc_SignalHandlerPop(SIGFPE,SIG_IGN); |
| 518 |
|
|
if (scaled) { |
| 519 |
|
|
if( relman_calc_satisfied_scaled(rel,epsilon) ) break; |
| 520 |
|
|
} else { |
| 521 |
|
|
if( relman_calc_satisfied(rel,epsilon) ) break; |
| 522 |
|
|
} |
| 523 |
|
|
var_set_value(var,save); |
| 524 |
|
|
} else { |
| 525 |
|
|
var_set_value(var,slist[nsolns]); |
| 526 |
|
|
Asc_SignalHandlerPush(SIGFPE,SIG_IGN); |
| 527 |
|
|
relman_eval(rel,&status,SAFE_FIX_ME); |
| 528 |
|
|
Asc_SignalHandlerPop(SIGFPE,SIG_IGN); |
| 529 |
|
|
break; |
| 530 |
|
|
} |
| 531 |
|
|
} |
| 532 |
|
|
if (nsolns<0 && allsolns>0 && fp !=NULL) { |
| 533 |
|
|
/* dump the rejected solutions to give the user a clue */ |
| 534 |
|
|
FPRINTF(fp,"Ignoring potential solutions for variable\n"); |
| 535 |
|
|
var_write_name(server,var,fp); |
| 536 |
|
|
FPRINTF(fp,"\nin equation\n"); |
| 537 |
|
|
rel_write_name(server,rel,fp); |
| 538 |
|
|
FPRINTF(fp,"\n"); |
| 539 |
|
|
for (--allsolns; allsolns >= 0; allsolns--) { |
| 540 |
|
|
FPRINTF(fp,"Rejected solution: %.18g\n",slist[allsolns]); |
| 541 |
|
|
} |
| 542 |
|
|
} |
| 543 |
|
|
/* destroy_array(slist); do not do this */ |
| 544 |
|
|
return( nsolns >= 0 ? 1 : -1 ); |
| 545 |
|
|
} |
| 546 |
|
|
|
| 547 |
|
|
|
| 548 |
|
|
int slv_direct_log_solve(slv_system_t server, struct logrel_relation *lrel, |
| 549 |
|
|
struct dis_discrete *dvar, FILE *fp, int perturb, |
| 550 |
|
|
struct gl_list_t *insts) |
| 551 |
|
|
/** |
| 552 |
|
|
*** Attempt to directly solve the given relation (equality constraint) for |
| 553 |
|
|
*** the given variable, leaving the others fixed. Returns an integer |
| 554 |
|
|
*** signifying the status as one of the following three: |
| 555 |
|
|
*** |
| 556 |
|
|
*** 0 ==> Unable to determine anything. Bad logrelation or dvar |
| 557 |
|
|
*** 1 ==> Solution found, discrete variable value set to this solution. |
| 558 |
|
|
*** 2 ==> More than one solution, Do not modify the current value |
| 559 |
|
|
*** of the variable and display an error message. |
| 560 |
|
|
*** -1 ==> No solution found. |
| 561 |
|
|
**/ |
| 562 |
|
|
{ |
| 563 |
|
|
int32 able; |
| 564 |
|
|
int32 nsolns, c; |
| 565 |
|
|
int32 *slist; |
| 566 |
|
|
|
| 567 |
|
|
slist = logrelman_directly_solve(lrel,dvar,&able,&nsolns,perturb,insts); |
| 568 |
|
|
if( !able ) return(0); |
| 569 |
|
|
if(nsolns == -1) return (-1); |
| 570 |
|
|
|
| 571 |
|
|
if (nsolns == 1) { |
| 572 |
|
|
dis_set_boolean_value(dvar,slist[1]); |
| 573 |
|
|
logrel_set_residual(lrel,TRUE); |
| 574 |
|
|
destroy_array(slist); |
| 575 |
|
|
return 1; |
| 576 |
|
|
} else { |
| 577 |
|
|
FPRINTF(fp,"Ignoring potential solutions for discrete variable\n"); |
| 578 |
|
|
dis_write_name(server,dvar,fp); |
| 579 |
|
|
FPRINTF(fp,"\nin equation\n"); |
| 580 |
|
|
logrel_write_name(server,lrel,fp); |
| 581 |
|
|
FPRINTF(fp,"\n"); |
| 582 |
|
|
for (c = nsolns; c >= 1; c--) { |
| 583 |
|
|
FPRINTF(fp,"Rejected solution: %d \n",slist[c]); |
| 584 |
|
|
} |
| 585 |
|
|
destroy_array(slist); /* should we have to do this? */ |
| 586 |
|
|
return 2; |
| 587 |
|
|
} |
| 588 |
|
|
} |
| 589 |
|
|
|
| 590 |
|
|
|
| 591 |
|
|
#endif |
| 592 |
|
|
|
| 593 |
|
|
int32 **slv_lnkmap_from_mtx(mtx_matrix_t mtx, int32 len, int32 m) |
| 594 |
|
|
{ |
| 595 |
|
|
int32 **map,*data,rl; |
| 596 |
|
|
real64 val; |
| 597 |
|
|
mtx_coord_t coord; |
| 598 |
|
|
|
| 599 |
|
|
if (mtx==NULL) { |
| 600 |
|
|
FPRINTF(stderr,"Warning: slv_lnkmap_from_mtx called with null mtx\n"); |
| 601 |
|
|
return NULL; |
| 602 |
|
|
} |
| 603 |
|
|
data=(int32 *)ascmalloc((m+2*len)*sizeof(int32)); |
| 604 |
|
|
map=(int **)ascmalloc(m*sizeof(void *)); |
| 605 |
|
|
for(coord.row=0; coord.row < m; coord.row ++) { |
| 606 |
|
|
rl=mtx_nonzeros_in_row(mtx,coord.row,mtx_ALL_COLS); |
| 607 |
|
|
map[coord.row]=data; |
| 608 |
|
|
data[0]=rl; |
| 609 |
|
|
data++; |
| 610 |
|
|
coord.col=mtx_FIRST; |
| 611 |
|
|
while( val=mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), |
| 612 |
|
|
coord.col != mtx_LAST) { |
| 613 |
|
|
data[0]=coord.col; |
| 614 |
|
|
data[1]= (int32)ascnint(val); |
| 615 |
|
|
data+=2; |
| 616 |
|
|
} |
| 617 |
|
|
} |
| 618 |
|
|
return map; |
| 619 |
|
|
} |
| 620 |
|
|
int32 **slv_create_lnkmap(int m, int n, int len, int *hi, int *hj) { |
| 621 |
|
|
mtx_matrix_t mtx; |
| 622 |
|
|
mtx_coord_t coord; |
| 623 |
|
|
int32 i, **map; |
| 624 |
|
|
|
| 625 |
|
|
mtx=mtx_create(); |
| 626 |
|
|
mtx_set_order(mtx,MAX(m,n)); |
| 627 |
|
|
for (i=0;i <len; i++) { |
| 628 |
|
|
coord.row=hi[i]; |
| 629 |
|
|
coord.col=hj[i]; |
| 630 |
|
|
mtx_fill_value(mtx,&coord,(real64)i); |
| 631 |
|
|
} |
| 632 |
|
|
map=slv_lnkmap_from_mtx(mtx,len,m); |
| 633 |
|
|
mtx_destroy(mtx); |
| 634 |
|
|
return map; |
| 635 |
|
|
} |
| 636 |
|
|
|
| 637 |
|
|
void slv_destroy_lnkmap(int32 **map) { |
| 638 |
|
|
if (NOTNULL(map)) { |
| 639 |
|
|
ascfree(map[0]); |
| 640 |
|
|
ascfree(map); |
| 641 |
|
|
} |
| 642 |
|
|
} |
| 643 |
|
|
|
| 644 |
|
|
void slv_write_lnkmap(FILE *fp, int32 m, int32 **map) { |
| 645 |
|
|
int32 i,j,nv, *v; |
| 646 |
|
|
if (ISNULL(map)) { |
| 647 |
|
|
FPRINTF(stderr,"slv_write_lnkmap: Cannot write NULL lnkmap\n"); |
| 648 |
|
|
return; |
| 649 |
|
|
} |
| 650 |
|
|
for(i=0; i<m; i++) { |
| 651 |
|
|
if (NOTNULL(map[i])) { |
| 652 |
|
|
v=map[i]; |
| 653 |
|
|
nv=v[0]; |
| 654 |
|
|
v++; |
| 655 |
|
|
FPRINTF(fp,"Row %d len = %d\n",i,nv); |
| 656 |
|
|
for (j=0;j<nv;j++) { |
| 657 |
|
|
FPRINTF(fp," Col %d, lnk location %d\n",v[0],v[1]); |
| 658 |
|
|
v+=2; |
| 659 |
|
|
} |
| 660 |
|
|
} |
| 661 |
|
|
} |
| 662 |
|
|
} |
| 663 |
|
|
|
Parent Directory
| 
Revision Log