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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 108 - (hide annotations) (download) (as text)
Tue Dec 13 05:53:20 2005 UTC (14 years, 10 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 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 jds 61 #include "utilities/ascPanic.h"
42 aw0a 1 #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 jds 61 *** slv_create_vector(low,high)
72     *** slv_init_vector(vec,low,high)
73     *** slv_destroy_vector(vec)
74 aw0a 1 *** 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 jds 61 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 aw0a 1 void slv_zero_vector( struct vector_data *vec)
145     {
146     real64 *p;
147     int32 len;
148    
149 jds 61 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 aw0a 1 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 jds 61 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 aw0a 1 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 jds 61 struct vector_data *vec2)
186 aw0a 1 /**
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 jds 61 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 aw0a 1 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 jds 61 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 aw0a 1 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 jds 61 void slv_write_vector(FILE *fp, struct vector_data *vec)
286 aw0a 1 /**
287     *** Outputs a vector.
288     **/
289     {
290     int32 ndx,hi;
291     real64 *vvec;
292 jds 61
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 aw0a 1 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 jds 61
353     asc_assert((NULL != p1) && (NULL != p2) && (len >= 0));
354    
355 aw0a 1 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 jds 61 *** NOT yet implemented correctly:
467 aw0a 1 *** 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 jds 61 FILE *slv_get_output_file(FILE *fp)
480 aw0a 1 {
481     static FILE *nuldev = NULL;
482 jds 61 #ifndef __WIN32__
483     const char fname[] = "/dev/null";
484     #else
485     const char fname[] = "nul";
486     #endif
487 aw0a 1
488 jds 61 if( fp == NULL ) {
489     if(nuldev == NULL)
490     if( (nuldev = fopen(fname,"w")) == NULL ) {
491     FPRINTF(stderr,"ERROR: slv_get_output_file\n");
492 aw0a 1 FPRINTF(stderr," Unable to open %s.\n",fname);
493     }
494 jds 61 fp = nuldev;
495 aw0a 1 }
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 jds 61 char *name = NULL;
506 aw0a 1 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 jds 61 void slv_print_dis_name( FILE *out, slv_system_t sys, struct dis_discrete *dvar)
523 aw0a 1 {
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 jds 61 void slv_print_obj_name(FILE *out, obj_objective_t obj)
544 aw0a 1 {
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 jds 61 void slv_print_var_sindex(FILE *out, struct var_variable *var)
554 aw0a 1 {
555     FPRINTF(out,"%d",var_sindex(var));
556     }
557    
558 jds 61 void slv_print_rel_sindex(FILE *out, struct rel_relation *rel)
559 aw0a 1 {
560     FPRINTF(out,"%d",rel_sindex(rel));
561     }
562    
563 jds 61 void slv_print_dis_sindex(FILE *out, struct dis_discrete *dvar)
564 aw0a 1 {
565     FPRINTF(out,"%d",dis_sindex(dvar));
566     }
567    
568 jds 61 void slv_print_logrel_sindex(FILE *out, struct logrel_relation *lrel)
569 aw0a 1 {
570     FPRINTF(out,"%d",logrel_sindex(lrel));
571     }
572    
573     #ifdef NEWSTUFF
574 jds 61 void slv_print_obj_index(FILE *out, obj_objective_t obj)
575 aw0a 1 {
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 jds 61 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 aw0a 1 {
598 jds 61 int32 able, status;
599     int nsolns, allsolns;
600     real64 *slist, save;
601 aw0a 1
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 jds 61 (void)relman_eval(rel,&status,SAFE_FIX_ME);
617 aw0a 1 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 jds 61 (void)relman_eval(rel,&status,SAFE_FIX_ME);
629 aw0a 1 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 jds 61 (void)relman_eval(rel,&status,SAFE_FIX_ME);
640 aw0a 1 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 johnpye 86 error_reporter_start(ASC_PROG_WARNING,NULL,0);
647 johnpye 108 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 aw0a 1 for (--allsolns; allsolns >= 0; allsolns--) {
653 johnpye 108 FPRINTF(ASCERR,"Rejected solution: %.18g\n",slist[allsolns]);
654 aw0a 1 }
655 johnpye 108 error_reporter_end_flush();
656 aw0a 1 }
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 jds 61 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 aw0a 1 {
677     int32 able;
678     int32 nsolns, c;
679     int32 *slist;
680    
681 johnpye 108 (void)fp;
682    
683 aw0a 1 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 johnpye 108 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 aw0a 1 for (c = nsolns; c >= 1; c--) {
700 johnpye 108 FPRINTF(ASCERR,"Rejected solution: %d \n",slist[c]);
701 aw0a 1 }
702 johnpye 108 error_reporter_end_flush();
703    
704 aw0a 1 destroy_array(slist); /* should we have to do this? */
705     return 2;
706     }
707     }
708    
709 jds 61 #endif /* SLV_INSTANCES */
710 aw0a 1
711 jds 61 int32 **slv_lnkmap_from_mtx(mtx_matrix_t mtx, mtx_region_t *clientregion)
712 aw0a 1 {
713 jds 61 int32 **map, *data, rl, order;
714 aw0a 1 real64 val;
715     mtx_coord_t coord;
716 jds 61 mtx_range_t range;
717     mtx_region_t region;
718 aw0a 1
719 jds 61 if (mtx == NULL) {
720 aw0a 1 FPRINTF(stderr,"Warning: slv_lnkmap_from_mtx called with null mtx\n");
721     return NULL;
722     }
723 jds 61
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 aw0a 1 data++;
761 jds 61 coord.col = mtx_FIRST;
762     while( val = mtx_next_in_row(mtx, &coord, &range),
763 aw0a 1 coord.col != mtx_LAST) {
764 jds 61 data[0] = coord.col;
765     data[1] = (int32)ascnint(val);
766     data += 2;
767 aw0a 1 }
768     }
769 jds 61 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 aw0a 1 return map;
775     }
776 jds 61
777     int32 **slv_create_lnkmap(int m, int n, int len, int *hi, int *hj)
778     {
779 aw0a 1 mtx_matrix_t mtx;
780     mtx_coord_t coord;
781     int32 i, **map;
782    
783 jds 61 mtx = mtx_create();
784     if (NULL == mtx) {
785     FPRINTF(stderr,"Warning: slv_create_lnkmap called with null mtx\n");
786     return NULL;
787     }
788 aw0a 1 mtx_set_order(mtx,MAX(m,n));
789 jds 61 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 aw0a 1 mtx_fill_value(mtx,&coord,(real64)i);
798     }
799 jds 61 map = slv_lnkmap_from_mtx(mtx,mtx_ENTIRE_MATRIX);
800 aw0a 1 mtx_destroy(mtx);
801     return map;
802     }
803    
804 jds 61
805     void slv_destroy_lnkmap(int32 **map)
806     {
807 aw0a 1 if (NOTNULL(map)) {
808     ascfree(map[0]);
809     ascfree(map);
810     }
811     }
812    
813 jds 61 void slv_write_lnkmap(FILE *fp, int32 m, int32 **map)
814     {
815 aw0a 1 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