/[ascend]/trunk/base/generic/solver/slv_common.c
ViewVC logotype

Contents of /trunk/base/generic/solver/slv_common.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 108 - (show annotations) (download) (as text)
Tue Dec 13 05:53:20 2005 UTC (14 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 22821 byte(s)
Moved some solver error messages to FPRINTF(ASCERR) convention.
Trying to get IAPWS-95 to work two-phase.
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 "utilities/ascPanic.h"
42 #include "solver/mtx.h"
43 #include "solver/slv_types.h"
44 #include "solver/rel.h"
45 #include "solver/var.h"
46 #include "solver/discrete.h"
47 #include "solver/logrel.h"
48 #include "solver/slv_common.h"
49 #include "utilities/mem.h"
50 /* if libasc.a running around, the following: */
51 #if SLV_INSTANCES
52 #include "compiler/fractions.h"
53 #include "compiler/dimen.h"
54 #include "compiler/functype.h"
55 #include "compiler/func.h"
56 #include "solver/relman.h"
57 #include "solver/logrelman.h"
58 #else
59 #ifdef NDEBUG
60 #define ascnint(a) (((int) (a)>=0.0 ? floor((a) + 0.5) : -floor(0.5 - (a))))
61 #else
62 #define ascnint(a) ascnintF(a)
63 #endif /* NDEBUG */
64 #endif /* instances */
65
66 #define SAFE_FIX_ME 0
67
68 /**
69 *** Array/vector operations
70 *** ----------------------------
71 *** slv_create_vector(low,high)
72 *** slv_init_vector(vec,low,high)
73 *** slv_destroy_vector(vec)
74 *** slv_zero_vector(vec)
75 *** slv_copy_vector(vec1,vec2)
76 *** prod = slv_inner_product(vec1,vec2)
77 *** norm2 = slv_square_norm(vec)
78 *** slv_matrix_product(mtx,vec,prod,scale,transpose)
79 *** slv_write_vector(file,vec)
80 **/
81
82 struct vector_data *slv_create_vector(int32 low, int32 high)
83 {
84 struct vector_data *result;
85
86 result = (struct vector_data *)ascmalloc(sizeof(struct vector_data));
87 if (NULL == result)
88 return NULL;
89
90 result->rng = NULL;
91 result->vec = NULL;
92 if (0 != slv_init_vector(result, low, high)) {
93 ascfree(result);
94 result = NULL;
95 }
96 return result;
97 }
98
99 int slv_init_vector(struct vector_data *vec, int32 low, int32 high)
100 {
101 int32 new_size;
102
103 if ((low < 0) || (high < low))
104 return 1;
105
106 if (NULL == vec)
107 return 2;
108
109 if (NULL == vec->rng) {
110 vec->rng = (mtx_range_t *)ascmalloc(sizeof(mtx_range_t));
111 if (NULL == vec->rng)
112 return 3;
113 }
114 vec->rng = mtx_range(vec->rng, low, high);
115
116 new_size = high + 1;
117 if (NULL == vec->vec) {
118 vec->vec = (real64 *)ascmalloc((new_size)*sizeof(real64));
119 if (NULL == vec->vec) {
120 ascfree(vec->rng);
121 vec->rng = NULL;
122 return 3;
123 }
124 }
125 else {
126 vec->vec = (real64 *)ascrealloc(vec->vec, (new_size)*sizeof(real64));
127 }
128
129 vec->accurate = FALSE;
130 return 0;
131 }
132
133 void slv_destroy_vector(struct vector_data *vec)
134 {
135 if (NULL != vec) {
136 if (NULL != vec->rng)
137 ascfree(vec->rng);
138 if (NULL != vec->vec)
139 ascfree(vec->vec);
140 ascfree(vec);
141 }
142 }
143
144 void slv_zero_vector( struct vector_data *vec)
145 {
146 real64 *p;
147 int32 len;
148
149 asc_assert((NULL != vec) &&
150 (NULL != vec->rng) &&
151 (NULL != vec->vec) &&
152 (vec->rng->low >= 0) &&
153 (vec->rng->low <= vec->rng->high));
154
155 p = vec->vec + vec->rng->low;
156 len = vec->rng->high - vec->rng->low + 1;
157 mtx_zero_real64(p,len);
158 }
159
160 void slv_copy_vector( struct vector_data *vec1,struct vector_data *vec2)
161 {
162 real64 *p1,*p2;
163 int32 len;
164
165 asc_assert((NULL != vec1) &&
166 (NULL != vec1->rng) &&
167 (NULL != vec1->vec) &&
168 (vec1->rng->low >= 0) &&
169 (vec1->rng->low <= vec1->rng->high) &&
170 (NULL != vec2) &&
171 (NULL != vec2->rng) &&
172 (NULL != vec2->vec) &&
173 (vec2->rng->low >= 0));
174
175 p1 = vec1->vec + vec1->rng->low;
176 p2 = vec2->vec + vec2->rng->low;
177 len = vec1->rng->high - vec1->rng->low + 1;
178 /*mem_copy_cast not in order here, probably */
179 mem_move_cast(p1,p2,len*sizeof(real64));
180 }
181
182 #define USEDOT TRUE
183 /* USEDOT = TRUE is a winner on alphas, hps, and sparc20 */
184 real64 slv_inner_product(struct vector_data *vec1 ,
185 struct vector_data *vec2)
186 /**
187 *** Computes inner product between vec1 and vec2, returning result.
188 *** vec1 and vec2 may overlap or even be identical.
189 **/
190 {
191 real64 *p1,*p2;
192 #if !USEDOT
193 real64 sum;
194 #endif
195 int32 len;
196
197 asc_assert((NULL != vec1) &&
198 (NULL != vec1->rng) &&
199 (NULL != vec1->vec) &&
200 (vec1->rng->low >= 0) &&
201 (vec1->rng->low <= vec1->rng->high) &&
202 (NULL != vec2) &&
203 (NULL != vec2->rng) &&
204 (NULL != vec2->vec) &&
205 (vec2->rng->low >= 0));
206
207 p1 = vec1->vec + vec1->rng->low;
208 p2 = vec2->vec + vec2->rng->low;
209 len = vec1->rng->high - vec1->rng->low + 1;
210 #if !USEDOT
211 if (p1 != p2) {
212 for( sum=0.0 ; --len >= 0 ; ++p1,++p2 ) {
213 sum += (*p1) * (*p2);
214 }
215 return(sum);
216 } else {
217 for( sum=0.0 ; --len >= 0 ; ++p1 ) {
218 sum += (*p1) * (*p1);
219 }
220 return(sum);
221 }
222 #else
223 return slv_dot(len,p1,p2);
224 #endif
225 }
226
227 real64 slv_square_norm(struct vector_data *vec)
228 /**
229 *** Computes norm^2 of vector, assigning the result to vec->norm2
230 *** and returning the result as well.
231 **/
232 {
233 vec->norm2 = slv_inner_product(vec,vec);
234 return vec->norm2;
235 }
236
237 void slv_matrix_product(mtx_matrix_t mtx, struct vector_data *vec,
238 struct vector_data *prod, real64 scale,
239 boolean transpose)
240 /**
241 *** Stores prod := (scale)*(mtx)(vec) or (scale)*(mtx-transpose)(vec).
242 *** vec and prod must be completely different.
243 **/
244 {
245 mtx_coord_t nz;
246 real64 value, *vvec, *pvec;
247 int32 lim;
248
249 asc_assert((NULL != vec) &&
250 (NULL != vec->rng) &&
251 (NULL != vec->vec) &&
252 (vec->rng->low >= 0) &&
253 (vec->rng->low <= vec->rng->high) &&
254 (NULL != prod) &&
255 (NULL != prod->rng) &&
256 (NULL != prod->vec) &&
257 (prod->rng->low >= 0) &&
258 (prod->rng->low <= prod->rng->high) &&
259 (NULL != mtx));
260
261 lim = prod->rng->high;
262 pvec = prod->vec;
263 vvec = vec->vec;
264 if( transpose ) {
265 for(nz.col = prod->rng->low ; nz.col <= lim ; ++(nz.col) ) {
266 pvec[nz.col] = 0.0;
267 nz.row = mtx_FIRST;
268 while( value = mtx_next_in_col(mtx,&nz,vec->rng),
269 nz.row != mtx_LAST )
270 pvec[nz.col] += value*vvec[nz.row];
271 pvec[nz.col] *= scale;
272 }
273 } else {
274 for(nz.row = prod->rng->low ; nz.row <= lim ; ++(nz.row) ) {
275 pvec[nz.row] = 0.0;
276 nz.col = mtx_FIRST;
277 while( value = mtx_next_in_row(mtx,&nz,vec->rng),
278 nz.col != mtx_LAST )
279 pvec[nz.row] += value*vvec[nz.col];
280 pvec[nz.row] *= scale;
281 }
282 }
283 }
284
285 void slv_write_vector(FILE *fp, struct vector_data *vec)
286 /**
287 *** Outputs a vector.
288 **/
289 {
290 int32 ndx,hi;
291 real64 *vvec;
292
293 if (NULL == fp) {
294 FPRINTF(ASCERR, "Error writing vector in slv_write_vector: NULL file pointer.\n");
295 return;
296 }
297 if ((NULL == vec) ||
298 (NULL == vec->rng) ||
299 (NULL == vec->vec) ||
300 (vec->rng->low < 0) ||
301 (vec->rng->low > vec->rng->high)) {
302 FPRINTF(ASCERR, "Error writing vector in slv_write_vector: uninitialized vector.\n");
303 return;
304 }
305
306 vvec = vec->vec;
307 hi = vec->rng->high;
308 FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n",
309 sqrt(fabs(vec->norm2)), vec->accurate?"TRUE":"FALSE",
310 vec->rng->low,vec->rng->high);
311 FPRINTF(fp,"Vector --> ");
312 for( ndx=vec->rng->low ; ndx<=hi ; ++ndx )
313 FPRINTF(fp, "%g ", vvec[ndx]);
314 PUTC('\n',fp);
315 }
316
317 /* Dot product for loop unrolled vector norms */
318 real64 slv_dot(int32 len, const real64 *p1, const real64 *p2)
319 {
320 register double sum,lsum;
321 int m,n;
322
323 /* AVMAGIC in fact isn't magic.
324 * only goes to 2-10. change the code below if you mess with it.
325 * Default AVMAGIC is 4, which works well on alphas, tika.
326 */
327 #define AVMAGIC 4
328 #ifdef sun
329 #undef AVMAGIC
330 #define AVMAGIC 6
331 #endif
332 #ifdef __hpux
333 #undef AVMAGIC
334 #define AVMAGIC 4
335 /* 2 was best value on ranier(9000/720) but tika (9000/715) likes 4
336 * there are no recognizable (defines) compiler differences, so tika
337 * will set the hp default
338 */
339 #endif
340 /* hands down best avmagic on atlas is 4, ranier 2, unxi21 6
341 * under native compilers. no bets on the billions of gcc variants.
342 * needs to be tried on sgi.
343 * Note, these are tuned for speed, not for accuracy. on very large
344 * vectors something like dnrm2 may be more appropriate, though
345 * it is very slow. upping avmagic to 10 will help accuracy some.
346 */
347
348 #if (AVMAGIC>10)
349 #undef AVMAGIC
350 #define AVMAGIC 10
351 #endif
352
353 asc_assert((NULL != p1) && (NULL != p2) && (len >= 0));
354
355 m = len / AVMAGIC;
356 n = len % AVMAGIC;
357 if (p1!=p2) {
358 /* get leading leftovers */
359 for( sum=0.0 ; --n >= 0 ; ) {
360 sum += (*p1) * (*p2);
361 ++p1; ++p2;
362 }
363 /* p1,p2 now point at first unadded location, or just after the end (m=0)*/
364 /* eat m chunks */
365 for( n=0; n <m; n++) {
366 /* now, as i am too lazy to figure out a macro that expands itself */
367 lsum = (*p1) * (*p2); /* zeroth term is assigned */
368 p1++; p2++;
369 lsum += (*p1) * (*p2); /* add 1th */
370 #if (AVMAGIC>2)
371 p1++; p2++;
372 lsum += (*p1) * (*p2); /* add 2th */
373 #if (AVMAGIC>3)
374 p1++; p2++;
375 lsum += (*p1) * (*p2); /* add 3th */
376 #if (AVMAGIC>4)
377 p1++; p2++;
378 lsum += (*p1) * (*p2); /* add 4th */
379 #if (AVMAGIC>5)
380 p1++; p2++;
381 lsum += (*p1) * (*p2); /* add 5th */
382 #if (AVMAGIC>6)
383 p1++; p2++;
384 lsum += (*p1) * (*p2); /* add 6th */
385 #if (AVMAGIC>7)
386 p1++; p2++;
387 lsum += (*p1) * (*p2); /* add 7th */
388 #if (AVMAGIC>8)
389 p1++; p2++;
390 lsum += (*p1) * (*p2); /* add 8th */
391 #if (AVMAGIC>9)
392 p1++; p2++;
393 lsum += (*p1) * (*p2); /* add 9th */
394 #endif
395 #endif
396 #endif
397 #endif
398 #endif
399 #endif
400 #endif
401 #endif
402 p1++; p2++; /* leave p1,p2 pointing at next zeroth */
403 sum += lsum;
404 }
405 } else {
406 /* get leading leftovers */
407 for( sum=0.0 ; --n >= 0 ; ++p1 ) {
408 sum += (*p1) * (*p1);
409 }
410 /* p1 now points at first unadded location, or just after the end (m=0)*/
411 /* eat m chunks */
412 for( n=0; n <m; n++) {
413 /* now, as i am too lazy to figure out a macro that expands itself */
414 lsum = (*p1) * (*p1); /* zeroth term is assigned */
415 p1++;
416 lsum += (*p1) * (*p1); /* add 1th */
417 #if (AVMAGIC>2)
418 p1++;
419 lsum += (*p1) * (*p1); /* add 2th */
420 #if (AVMAGIC>3)
421 p1++;
422 lsum += (*p1) * (*p1); /* add 3th */
423 #if (AVMAGIC>4)
424 p1++;
425 lsum += (*p1) * (*p1); /* add 4th */
426 #if (AVMAGIC>5)
427 p1++;
428 lsum += (*p1) * (*p1); /* add 5th */
429 #if (AVMAGIC>6)
430 p1++;
431 lsum += (*p1) * (*p1); /* add 6th */
432 #if (AVMAGIC>7)
433 p1++;
434 lsum += (*p1) * (*p1); /* add 7th */
435 #if (AVMAGIC>8)
436 p1++;
437 lsum += (*p1) * (*p1); /* add 8th */
438 #if (AVMAGIC>9)
439 p1++;
440 lsum += (*p1) * (*p1); /* add 9th */
441 #endif
442 #endif
443 #endif
444 #endif
445 #endif
446 #endif
447 #endif
448 #endif
449 p1++; /* leave p1 pointing at next zeroth */
450 sum += lsum;
451 }
452 }
453 return(sum);
454 #undef AVMAGIC
455 }
456
457 /**
458 *** General input/output routines
459 *** -----------------------------
460 *** fp = MIF(sys)
461 *** fp = LIF(sys)
462 *** slv_print_var_name(out,sys,var)
463 *** slv_print_rel_name(out,sys,rel)
464 *** slv_print_dis_name(out,sys,dvar)
465 *** slv_print_logrel_name(out,sys,lrel)
466 *** NOT yet implemented correctly:
467 *** slv_print_obj_name(out,obj)
468 *** slv_print_var_sindex(out,var)
469 *** slv_print_rel_sindex(out,rel)
470 *** slv_print_dis_sindex(out,dvar)
471 *** slv_print_logrel_sindex(out,lrel)
472 *** slv_print_obj_index(out,obj)
473 **/
474
475 /**
476 *** Returns fp if fp!=NULL, or a file pointer
477 *** open to nul device if fp == NULL.
478 **/
479 FILE *slv_get_output_file(FILE *fp)
480 {
481 static FILE *nuldev = NULL;
482 #ifndef __WIN32__
483 const char fname[] = "/dev/null";
484 #else
485 const char fname[] = "nul";
486 #endif
487
488 if( fp == NULL ) {
489 if(nuldev == NULL)
490 if( (nuldev = fopen(fname,"w")) == NULL ) {
491 FPRINTF(stderr,"ERROR: slv_get_output_file\n");
492 FPRINTF(stderr," Unable to open %s.\n",fname);
493 }
494 fp = nuldev;
495 }
496 return(fp);
497 }
498 #define MIF(sys) slv_get_output_file( (sys)->p.output.more_important )
499 #define LIF(sys) slv_get_output_file( (sys)->p.output.less_important )
500
501 #if SLV_INSTANCES
502
503 void slv_print_var_name( FILE *out,slv_system_t sys, struct var_variable *var)
504 {
505 char *name = NULL;
506 if (out == NULL || sys == NULL || var == NULL) return;
507 name = var_make_name(sys,var);
508 if( *name == '?' ) FPRINTF(out,"%d",var_sindex(var));
509 else FPRINTF(out,"%s",name);
510 if (name) ascfree(name);
511 }
512
513 void slv_print_rel_name( FILE *out, slv_system_t sys, struct rel_relation *rel)
514 {
515 char *name;
516 name = rel_make_name(sys,rel);
517 if( *name == '?' ) FPRINTF(out,"%d",rel_sindex(rel));
518 else FPRINTF(out,"%s",name);
519 ascfree(name);
520 }
521
522 void slv_print_dis_name( FILE *out, slv_system_t sys, struct dis_discrete *dvar)
523 {
524 char *name=NULL;
525 if (out == NULL || sys == NULL || dvar == NULL) return;
526 name = dis_make_name(sys,dvar);
527 if( *name == '?' ) FPRINTF(out,"%d",dis_sindex(dvar));
528 else FPRINTF(out,"%s",name);
529 if (name) ascfree(name);
530 }
531
532 void slv_print_logrel_name( FILE *out, slv_system_t sys,
533 struct logrel_relation *lrel)
534 {
535 char *name;
536 name = logrel_make_name(sys,lrel);
537 if( *name == '?' ) FPRINTF(out,"%d",logrel_sindex(lrel));
538 else FPRINTF(out,"%s",name);
539 ascfree(name);
540 }
541
542 #ifdef NEWSTUFF
543 void slv_print_obj_name(FILE *out, obj_objective_t obj)
544 {
545 char *name;
546 name = obj_make_name(obj);
547 if( *name == '?' ) FPRINTF(out,"%d",obj_index(obj));
548 else FPRINTF(out,"%s",name);
549 ascfree(name);
550 }
551 #endif
552
553 void slv_print_var_sindex(FILE *out, struct var_variable *var)
554 {
555 FPRINTF(out,"%d",var_sindex(var));
556 }
557
558 void slv_print_rel_sindex(FILE *out, struct rel_relation *rel)
559 {
560 FPRINTF(out,"%d",rel_sindex(rel));
561 }
562
563 void slv_print_dis_sindex(FILE *out, struct dis_discrete *dvar)
564 {
565 FPRINTF(out,"%d",dis_sindex(dvar));
566 }
567
568 void slv_print_logrel_sindex(FILE *out, struct logrel_relation *lrel)
569 {
570 FPRINTF(out,"%d",logrel_sindex(lrel));
571 }
572
573 #ifdef NEWSTUFF
574 void slv_print_obj_index(FILE *out, obj_objective_t obj)
575 {
576 FPRINTF(out,"%d",obj_index(obj));
577 }
578 #endif
579
580 #define destroy_array(p) \
581 if( (p) != NULL ) ascfree((p))
582
583 /**
584 *** Attempt to directly solve the given relation (equality constraint) for
585 *** the given variable, leaving the others fixed. Returns an integer
586 *** signifying the status as one of the following three:
587 ***
588 *** 0 ==> Unable to determine anything, possibly due to fp exception.
589 *** 1 ==> Solution found, variable value set to this solution.
590 *** -1 ==> No solution found.
591 ***
592 *** The variable bounds will be upheld, unless they are to be ignored.
593 **/
594 int slv_direct_solve(slv_system_t server, struct rel_relation *rel,
595 struct var_variable *var, FILE *fp,
596 real64 epsilon, int ignore_bounds, int scaled)
597 {
598 int32 able, status;
599 int nsolns, allsolns;
600 real64 *slist, save;
601
602 slist = relman_directly_solve_new(rel,var,&able,&nsolns,epsilon);
603 if( !able ) {
604 return(0);
605 }
606 allsolns = nsolns;
607 while( --nsolns >= 0 ) {
608 if( ignore_bounds ) {
609 var_set_value(var,slist[nsolns]);
610 break;
611 }
612 if( var_lower_bound(var) > slist[nsolns] ) {
613 save = var_value(var);
614 var_set_value(var,var_lower_bound(var));
615 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
616 (void)relman_eval(rel,&status,SAFE_FIX_ME);
617 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
618 if (scaled) {
619 if( relman_calc_satisfied_scaled(rel,epsilon) ) break;
620 } else {
621 if( relman_calc_satisfied(rel,epsilon) ) break;
622 }
623 var_set_value(var,save);
624 } else if( var_upper_bound(var) < slist[nsolns] ) {
625 save = var_value(var);
626 var_set_value(var,var_upper_bound(var));
627 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
628 (void)relman_eval(rel,&status,SAFE_FIX_ME);
629 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
630 if (scaled) {
631 if( relman_calc_satisfied_scaled(rel,epsilon) ) break;
632 } else {
633 if( relman_calc_satisfied(rel,epsilon) ) break;
634 }
635 var_set_value(var,save);
636 } else {
637 var_set_value(var,slist[nsolns]);
638 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
639 (void)relman_eval(rel,&status,SAFE_FIX_ME);
640 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
641 break;
642 }
643 }
644 if (nsolns<0 && allsolns>0 && fp !=NULL) {
645 /* dump the rejected solutions to give the user a clue */
646 error_reporter_start(ASC_PROG_WARNING,NULL,0);
647 FPRINTF(ASCERR,"Ignoring potential solutions for variable '");
648 var_write_name(server,var,ASCERR);
649 FPRINTF(ASCERR,"' in equation '");
650 rel_write_name(server,rel,ASCERR);
651 FPRINTF(ASCERR,"'. ");
652 for (--allsolns; allsolns >= 0; allsolns--) {
653 FPRINTF(ASCERR,"Rejected solution: %.18g\n",slist[allsolns]);
654 }
655 error_reporter_end_flush();
656 }
657 /* destroy_array(slist); do not do this */
658 return( nsolns >= 0 ? 1 : -1 );
659 }
660
661
662 /**
663 *** Attempt to directly solve the given relation (equality constraint) for
664 *** the given variable, leaving the others fixed. Returns an integer
665 *** signifying the status as one of the following three:
666 ***
667 *** 0 ==> Unable to determine anything. Bad logrelation or dvar
668 *** 1 ==> Solution found, discrete variable value set to this solution.
669 *** 2 ==> More than one solution, Do not modify the current value
670 *** of the variable and display an error message.
671 *** -1 ==> No solution found.
672 **/
673 int slv_direct_log_solve(slv_system_t server, struct logrel_relation *lrel,
674 struct dis_discrete *dvar, FILE *fp, int perturb,
675 struct gl_list_t *insts)
676 {
677 int32 able;
678 int32 nsolns, c;
679 int32 *slist;
680
681 (void)fp;
682
683 slist = logrelman_directly_solve(lrel,dvar,&able,&nsolns,perturb,insts);
684 if( !able ) return(0);
685 if(nsolns == -1) return (-1);
686
687 if (nsolns == 1) {
688 dis_set_boolean_value(dvar,slist[1]);
689 logrel_set_residual(lrel,TRUE);
690 destroy_array(slist);
691 return 1;
692 } else {
693 error_reporter_start(ASC_PROG_ERROR,NULL,0);
694 FPRINTF(ASCERR,"Ignoring potential solutions for discrete variable '");
695 dis_write_name(server,dvar,ASCERR);
696 FPRINTF(ASCERR,"' in equation '");
697 logrel_write_name(server,lrel,ASCERR);
698 FPRINTF(ASCERR,"'. ");
699 for (c = nsolns; c >= 1; c--) {
700 FPRINTF(ASCERR,"Rejected solution: %d \n",slist[c]);
701 }
702 error_reporter_end_flush();
703
704 destroy_array(slist); /* should we have to do this? */
705 return 2;
706 }
707 }
708
709 #endif /* SLV_INSTANCES */
710
711 int32 **slv_lnkmap_from_mtx(mtx_matrix_t mtx, mtx_region_t *clientregion)
712 {
713 int32 **map, *data, rl, order;
714 real64 val;
715 mtx_coord_t coord;
716 mtx_range_t range;
717 mtx_region_t region;
718
719 if (mtx == NULL) {
720 FPRINTF(stderr,"Warning: slv_lnkmap_from_mtx called with null mtx\n");
721 return NULL;
722 }
723
724 order = mtx_order(mtx);
725 if (clientregion != NULL) {
726 if ((clientregion->row.low < 0) ||
727 (clientregion->row.high > (order-1)) ||
728 (clientregion->col.low < 0) ||
729 (clientregion->col.high > (order-1))) {
730 FPRINTF(stderr,"Warning: slv_lnkmap_from_mtx called with null or invalid region\n");
731 return NULL;
732 }
733 mtx_region(&region, clientregion->row.low, clientregion->row.high,
734 clientregion->col.low, clientregion->col.high);
735 } else {
736 mtx_region(&region, 0, order-1, 0, order-1);
737 }
738
739 range.low = region.col.low;
740 range.high = region.col.high;
741
742 data = (int32 *)ascmalloc((order+2*mtx_nonzeros_in_region(mtx, &region))*sizeof(int32));
743 if (NULL == data) {
744 return NULL;
745 }
746 map = (int32 **)ascmalloc(order*sizeof(int32 *));
747 if (NULL == map) {
748 ascfree(data);
749 return NULL;
750 }
751 for (coord.row=0 ; coord.row<region.row.low ; ++coord.row) {
752 map[coord.row] = data;
753 data[0] = 0;
754 ++data;
755 }
756 for(coord.row=region.row.low ; coord.row <= region.row.high ; coord.row ++) {
757 rl = mtx_nonzeros_in_row(mtx, coord.row, &range);
758 map[coord.row] = data;
759 data[0] = rl;
760 data++;
761 coord.col = mtx_FIRST;
762 while( val = mtx_next_in_row(mtx, &coord, &range),
763 coord.col != mtx_LAST) {
764 data[0] = coord.col;
765 data[1] = (int32)ascnint(val);
766 data += 2;
767 }
768 }
769 for (coord.row=(region.row.high+1) ; coord.row<order ; ++coord.row) {
770 map[coord.row] = data;
771 data[0] = 0;
772 ++data;
773 }
774 return map;
775 }
776
777 int32 **slv_create_lnkmap(int m, int n, int len, int *hi, int *hj)
778 {
779 mtx_matrix_t mtx;
780 mtx_coord_t coord;
781 int32 i, **map;
782
783 mtx = mtx_create();
784 if (NULL == mtx) {
785 FPRINTF(stderr,"Warning: slv_create_lnkmap called with null mtx\n");
786 return NULL;
787 }
788 mtx_set_order(mtx,MAX(m,n));
789 for (i=0 ; i<len ; i++) {
790 if ((hi[i] >= m) || (hj[i] >= n)) {
791 FPRINTF(stderr,"Warning: index out of range in slv_create_lnkmap\n");
792 mtx_destroy(mtx);
793 return NULL;
794 }
795 coord.row = hi[i];
796 coord.col = hj[i];
797 mtx_fill_value(mtx,&coord,(real64)i);
798 }
799 map = slv_lnkmap_from_mtx(mtx,mtx_ENTIRE_MATRIX);
800 mtx_destroy(mtx);
801 return map;
802 }
803
804
805 void slv_destroy_lnkmap(int32 **map)
806 {
807 if (NOTNULL(map)) {
808 ascfree(map[0]);
809 ascfree(map);
810 }
811 }
812
813 void slv_write_lnkmap(FILE *fp, int32 m, int32 **map)
814 {
815 int32 i,j,nv, *v;
816 if (ISNULL(map)) {
817 FPRINTF(stderr,"slv_write_lnkmap: Cannot write NULL lnkmap\n");
818 return;
819 }
820 for(i=0; i<m; i++) {
821 if (NOTNULL(map[i])) {
822 v=map[i];
823 nv=v[0];
824 v++;
825 FPRINTF(fp,"Row %d len = %d\n",i,nv);
826 for (j=0;j<nv;j++) {
827 FPRINTF(fp," Col %d, lnk location %d\n",v[0],v[1]);
828 v+=2;
829 }
830 }
831 }
832 }
833

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22