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

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