/[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 126 - (hide annotations) (download) (as text)
Tue Dec 20 13:27:11 2005 UTC (14 years, 10 months ago) by johnpye
File MIME type: text/x-csrc
File size: 22613 byte(s)
Reformating comments in the SLV files
1 aw0a 1 /*
2 johnpye 126 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 aw0a 1
9 johnpye 126 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 aw0a 1
14 johnpye 126 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 aw0a 1 #include <math.h>
31     #include "utilities/ascConfig.h"
32     #include "utilities/ascSignal.h"
33     #include "compiler/compiler.h"
34     #include "utilities/ascMalloc.h"
35 jds 61 #include "utilities/ascPanic.h"
36 aw0a 1 #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 jds 61 *** slv_create_vector(low,high)
66     *** slv_init_vector(vec,low,high)
67     *** slv_destroy_vector(vec)
68 aw0a 1 *** 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 jds 61 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 aw0a 1 void slv_zero_vector( struct vector_data *vec)
139     {
140     real64 *p;
141     int32 len;
142    
143 jds 61 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 aw0a 1 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 jds 61 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 aw0a 1 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 jds 61 struct vector_data *vec2)
180 aw0a 1 /**
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 jds 61 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 aw0a 1 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 jds 61 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 aw0a 1 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 jds 61 void slv_write_vector(FILE *fp, struct vector_data *vec)
280 aw0a 1 /**
281     *** Outputs a vector.
282     **/
283     {
284     int32 ndx,hi;
285     real64 *vvec;
286 jds 61
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 aw0a 1 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 jds 61
347     asc_assert((NULL != p1) && (NULL != p2) && (len >= 0));
348    
349 aw0a 1 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 jds 61 *** NOT yet implemented correctly:
461 aw0a 1 *** 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 jds 61 FILE *slv_get_output_file(FILE *fp)
474 aw0a 1 {
475     static FILE *nuldev = NULL;
476 jds 61 #ifndef __WIN32__
477     const char fname[] = "/dev/null";
478     #else
479     const char fname[] = "nul";
480     #endif
481 aw0a 1
482 jds 61 if( fp == NULL ) {
483     if(nuldev == NULL)
484     if( (nuldev = fopen(fname,"w")) == NULL ) {
485     FPRINTF(stderr,"ERROR: slv_get_output_file\n");
486 aw0a 1 FPRINTF(stderr," Unable to open %s.\n",fname);
487     }
488 jds 61 fp = nuldev;
489 aw0a 1 }
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 jds 61 char *name = NULL;
500 aw0a 1 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 jds 61 void slv_print_dis_name( FILE *out, slv_system_t sys, struct dis_discrete *dvar)
517 aw0a 1 {
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 jds 61 void slv_print_obj_name(FILE *out, obj_objective_t obj)
538 aw0a 1 {
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 jds 61 void slv_print_var_sindex(FILE *out, struct var_variable *var)
548 aw0a 1 {
549     FPRINTF(out,"%d",var_sindex(var));
550     }
551    
552 jds 61 void slv_print_rel_sindex(FILE *out, struct rel_relation *rel)
553 aw0a 1 {
554     FPRINTF(out,"%d",rel_sindex(rel));
555     }
556    
557 jds 61 void slv_print_dis_sindex(FILE *out, struct dis_discrete *dvar)
558 aw0a 1 {
559     FPRINTF(out,"%d",dis_sindex(dvar));
560     }
561    
562 jds 61 void slv_print_logrel_sindex(FILE *out, struct logrel_relation *lrel)
563 aw0a 1 {
564     FPRINTF(out,"%d",logrel_sindex(lrel));
565     }
566    
567     #ifdef NEWSTUFF
568 jds 61 void slv_print_obj_index(FILE *out, obj_objective_t obj)
569 aw0a 1 {
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 jds 61 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 aw0a 1 {
592 jds 61 int32 able, status;
593     int nsolns, allsolns;
594     real64 *slist, save;
595 aw0a 1
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 jds 61 (void)relman_eval(rel,&status,SAFE_FIX_ME);
611 aw0a 1 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 jds 61 (void)relman_eval(rel,&status,SAFE_FIX_ME);
623 aw0a 1 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 jds 61 (void)relman_eval(rel,&status,SAFE_FIX_ME);
634 aw0a 1 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 johnpye 86 error_reporter_start(ASC_PROG_WARNING,NULL,0);
641 johnpye 108 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 aw0a 1 for (--allsolns; allsolns >= 0; allsolns--) {
647 johnpye 108 FPRINTF(ASCERR,"Rejected solution: %.18g\n",slist[allsolns]);
648 aw0a 1 }
649 johnpye 108 error_reporter_end_flush();
650 aw0a 1 }
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 jds 61 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 aw0a 1 {
671     int32 able;
672     int32 nsolns, c;
673     int32 *slist;
674    
675 johnpye 108 (void)fp;
676    
677 aw0a 1 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 johnpye 108 error_reporter_start(ASC_PROG_ERROR,NULL,0);
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 aw0a 1 for (c = nsolns; c >= 1; c--) {
694 johnpye 108 FPRINTF(ASCERR,"Rejected solution: %d \n",slist[c]);
695 aw0a 1 }
696 johnpye 108 error_reporter_end_flush();
697    
698 aw0a 1 destroy_array(slist); /* should we have to do this? */
699     return 2;
700     }
701     }
702    
703 jds 61 #endif /* SLV_INSTANCES */
704 aw0a 1
705 jds 61 int32 **slv_lnkmap_from_mtx(mtx_matrix_t mtx, mtx_region_t *clientregion)
706 aw0a 1 {
707 jds 61 int32 **map, *data, rl, order;
708 aw0a 1 real64 val;
709     mtx_coord_t coord;
710 jds 61 mtx_range_t range;
711     mtx_region_t region;
712 aw0a 1
713 jds 61 if (mtx == NULL) {
714 aw0a 1 FPRINTF(stderr,"Warning: slv_lnkmap_from_mtx called with null mtx\n");
715     return NULL;
716     }
717 jds 61
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 aw0a 1 data++;
755 jds 61 coord.col = mtx_FIRST;
756     while( val = mtx_next_in_row(mtx, &coord, &range),
757 aw0a 1 coord.col != mtx_LAST) {
758 jds 61 data[0] = coord.col;
759     data[1] = (int32)ascnint(val);
760     data += 2;
761 aw0a 1 }
762     }
763 jds 61 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 aw0a 1 return map;
769     }
770 jds 61
771     int32 **slv_create_lnkmap(int m, int n, int len, int *hi, int *hj)
772     {
773 aw0a 1 mtx_matrix_t mtx;
774     mtx_coord_t coord;
775     int32 i, **map;
776    
777 jds 61 mtx = mtx_create();
778     if (NULL == mtx) {
779     FPRINTF(stderr,"Warning: slv_create_lnkmap called with null mtx\n");
780     return NULL;
781     }
782 aw0a 1 mtx_set_order(mtx,MAX(m,n));
783 jds 61 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 aw0a 1 mtx_fill_value(mtx,&coord,(real64)i);
792     }
793 jds 61 map = slv_lnkmap_from_mtx(mtx,mtx_ENTIRE_MATRIX);
794 aw0a 1 mtx_destroy(mtx);
795     return map;
796     }
797    
798 jds 61
799     void slv_destroy_lnkmap(int32 **map)
800     {
801 aw0a 1 if (NOTNULL(map)) {
802     ascfree(map[0]);
803     ascfree(map);
804     }
805     }
806    
807 jds 61 void slv_write_lnkmap(FILE *fp, int32 m, int32 **map)
808     {
809 aw0a 1 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