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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1271 - (show annotations) (download) (as text)
Sun Feb 4 02:00:01 2007 UTC (15 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 20981 byte(s)
Added bounds checking to IDA.
Suppressed some debug output.
TestSteam.testpeturbida is broken now.
1 /* ASCEND modelling environment
2 Copyright (C) 1998, 2006, 2007 Carnegie Mellon University
3 Copyright (C) 1996 Benjamin Andrew Allan
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2, or (at your option)
8 any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330,
18 Boston, MA 02111-1307, USA.
19 *//*
20 by Benjamin Andrew Allan
21 5/19/96
22 Last in CVS: $Revision: 1.28 $ $Date: 1998/06/16 16:53:04 $ $Author: mthomas $
23 */
24
25 #include "slv_stdcalls.h"
26
27 #include <utilities/ascMalloc.h>
28 #include <general/mathmacros.h>
29
30 #include "relman.h"
31 #include "logrelman.h"
32
33 /* headers of registered clients */
34 #include "slv0.h"
35 #include "slv1.h"
36 #include "slv2.h"
37 #include "slv3.h"
38 #include "slv6.h"
39 #include "slv7.h"
40 #include "slv8.h"
41 #include "slv9.h"
42 #include "slv9a.h"
43
44 #define MIMDEBUG 0 /* slv_make_incidence_mtx debugging */
45 #define SLBPDEBUG 0 /* slv_log_block_partition debugging */
46
47 /*
48 Here we play some hefty Jacobian reordering games.
49
50 What we want to happen eventually is as follows:
51
52 - Get the free & incident pattern for include relations.
53 - Output-assign the Jacobian.
54 - BLT permute the Jacobian. If underspecified, fake rows
55 to make things appear square.
56 - For all leading square blocks apply reordering (kirk, etc),
57 - If trailing rectangular block, "apply clever kirkbased scheme not known."???
58
59 At present, we aren't quite so clever. We:
60
61 - Get the free & incident pattern for include relations.
62 - Output-assign the Jacobian.
63 - BLT permute the square output assigned region of the Jacobian.
64 - For all leading square blocks apply reordering (kirk, etc)
65 - Collect the block list as part of the master data structure.
66 - Set sindices for rels, vars as current rows/cols in matrix so
67 that the jacobian is 'naturally' preordered for all solvers.
68
69 Solvers are still free to reorder their own matrices any way they like.
70 It's probably a dumb idea, though.
71 */
72
73 /*------------------------------------------------------------------------------
74 MATRIX CREATION
75 */
76
77 /* see slv_stdcalls.h */
78 int slv_make_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
79 var_filter_t *vf,rel_filter_t *rf)
80 {
81 #if MIMDEBUG
82 FILE *fp;
83 #endif
84 int32 r, rlen,ord;
85 struct rel_relation **rp;
86
87 if (sys==NULL || mtx == NULL || vf == NULL || rf == NULL) {
88 ERROR_REPORTER_HERE(ASC_PROG_ERR,"called with null");
89 return 1;
90 }
91 rp = slv_get_solvers_rel_list(sys);
92 assert(rp!=NULL);
93 rlen = slv_get_num_solvers_rels(sys);
94 ord = MAX(rlen,slv_get_num_solvers_vars(sys));
95 if (ord > mtx_order(mtx)) {
96 ERROR_REPORTER_HERE(ASC_PROG_ERR,"undersized matrix");
97 return 2;
98 }
99 for (r=0; r < rlen; r++) {
100 if (rel_apply_filter(rp[r],rf)) {
101 relman_get_incidence(rp[r],vf,mtx);
102 }
103 }
104 #if MIMDEBUG
105 fp = fopen("/tmp/mim.plot","w+");
106 if (fp !=NULL) {
107 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
108 fclose(fp);
109 }
110 #endif
111 return 0;
112 }
113
114 void slv_sort_rels_and_vars(slv_system_t sys
115 ,int32 *rel_count, int32 *var_count
116 ){
117 struct rel_relation **rp, **rtmp, *rel;
118 struct var_variable **vp, **vtmp, *var;
119
120 int32 nrow,ncol,rlen,vlen,rel_tmp_end;
121 int32 r,c,var_tmp_end,len;
122 var_filter_t vf;
123 rel_filter_t rf;
124
125 /* get rel and var info */
126 rp = slv_get_solvers_rel_list(sys);
127 vp = slv_get_solvers_var_list(sys);
128 rlen = slv_get_num_solvers_rels(sys);
129 vlen = slv_get_num_solvers_vars(sys);
130 if (rlen ==0 || vlen == 0) return;
131
132 assert(var_count != NULL && rel_count != NULL);
133 *var_count = *rel_count = -1;
134
135 /* allocate temp arrays */
136 vtmp = ASC_NEW_ARRAY(struct var_variable *,vlen);
137 if (vtmp == NULL) {
138 return;
139 }
140 rtmp = ASC_NEW_ARRAY(struct rel_relation *,rlen);
141 if (rtmp == NULL) {
142 ascfree(vtmp);
143 return;
144 }
145
146 /* set up filters */
147 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
148 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
149 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
150 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
151
152 /* count rows and cols */
153 nrow = slv_count_solvers_rels(sys,&rf);
154 ncol = slv_count_solvers_vars(sys,&vf);
155
156
157 /* Sort out unincluded and inactive rels to
158 * the end of the temporary relation list.
159 */
160 *rel_count = 0;
161 rel_tmp_end = rlen -1;
162 for (r=0; r < rlen; r++) {
163 rel = rp[r];
164 if (rel_apply_filter(rel,&rf)) {
165 rtmp[*rel_count] = rel;
166 (*rel_count)++;
167 } else {
168 rtmp[rel_tmp_end] = rel;
169 rel_tmp_end--;
170 }
171 }
172
173 /* sort the inactive out of the unincluded and move to end */
174 len = rlen -1;
175 for(c = rel_tmp_end + 1; len > c; ) {
176 rel = rtmp[c];
177 if (rel_active(rel)) {
178 c++;
179 } else {
180 rtmp[c] = rtmp[len];
181 rtmp[len] = rel;
182 len--;
183 }
184 }
185
186 /* Sort variable list */
187 *var_count = 0;
188 var_tmp_end = vlen - 1;
189 for (c = 0; c < vlen; c++) {
190 var = vp[c];
191 if (var_apply_filter(var,&vf)) {
192 vtmp[*var_count] = var;
193 (*var_count)++;
194 } else {
195 vtmp[var_tmp_end] = var;
196 var_tmp_end--;
197 }
198 }
199 /* sort the inactive out of the unincluded and move to end */
200 len = vlen -1;
201 for(c = var_tmp_end + 1; len > c; ) {
202 var = vtmp[c];
203 if (var_active(var) && var_active(var)) {
204 c++;
205 } else {
206 vtmp[c] = vtmp[len];
207 vtmp[len] = var;
208 len--;
209 }
210 }
211
212 /* reset solver indicies */
213 for (c = 0; c < vlen; c++) {
214 vp[c] = vtmp[c];
215 var_set_sindex(vp[c],c);
216 }
217
218 for (c = 0; c < rlen; c++) {
219 rp[c] = rtmp[c];
220 rel_set_sindex(rp[c],c);
221 }
222 ascfree(vtmp);
223 ascfree(rtmp);
224 return;
225 }
226
227 /*------------------------------------------------------------------------------
228 QUERYING AND ENFORCING BOUNDS
229 */
230
231 int slv_ensure_bounds(slv_system_t sys,int32 lo,int32 hi, FILE *mif)
232 {
233 real64 val,low,high;
234 int32 c,nchange=0;
235 struct var_variable *var, **vp;
236
237 vp = slv_get_solvers_var_list(sys);
238 if (vp==NULL) return -1;
239 for (c= lo; c <= hi; c++) {
240 var = vp[c];
241 low = var_lower_bound(var);
242 high = var_upper_bound(var);
243 val = var_value(var);
244 if( low > high ) {
245 if (mif!=NULL) {
246 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
247 FPRINTF(ASCERR,"Bounds for variable '");
248 var_write_name(sys,var,ASCERR);
249 FPRINTF(ASCERR,"' are inconsistent [%g,%g].\n",low,high);
250 FPRINTF(ASCERR,"Bounds will be swapped.\n");
251 error_reporter_end_flush();
252 }
253 var_set_upper_bound(var, low);
254 var_set_lower_bound(var, high);
255 low = var_lower_bound(var);
256 high = var_upper_bound(var);
257 nchange++;
258 }
259
260 if( low > val ) {
261 if (mif!=NULL) {
262 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
263 FPRINTF(ASCERR,"Variable '");
264 var_write_name(sys,var,ASCERR);
265 FPRINTF(ASCERR,"' was set below its lower bound.\n");
266 FPRINTF(ASCERR,"It will be moved to its lower bound.");
267 error_reporter_end_flush();
268 }
269 var_set_value(var, low);
270 } else {
271 if( val > high ) {
272 if (mif!=NULL) {
273 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
274 FPRINTF(ASCERR,"Variable '");
275 var_write_name(sys,var,ASCERR);
276 FPRINTF(ASCERR,"' was set above its upper bound.\n");
277 FPRINTF(ASCERR,"It will be moved to its upper bound.");
278 error_reporter_end_flush();
279 }
280 var_set_value(var, high);
281 nchange++;
282 }
283 }
284 }
285 return nchange;
286 }
287
288 /* return 0 on success (ie bounds are met) */
289 int slv_check_bounds(const slv_system_t sys
290 , int32 lo,int32 hi
291 , const char *label
292 ){
293 real64 val,low,high;
294 int32 c,len;
295 struct var_variable *var, **vp;
296 int err = 0;
297
298 if(label==NULL) label = "";
299 if(sys==NULL) return -1;
300 vp = slv_get_solvers_var_list(sys);
301 if(vp==NULL) return -2;
302 len = slv_get_num_solvers_vars(sys);
303 if(hi < 0)hi+= len; /* so you can use -1 to mean 'the last' */
304 if(lo < 0)lo+= len;
305 if(lo > len || hi > len || lo < 0 || hi < 0 || lo > hi){
306 ERROR_REPORTER_HERE(ASC_PROG_ERR,"slv_check_bounds miscalled");
307 return -1;
308 }
309
310 for (c= lo; c <= hi; c++) {
311 var = vp[c];
312 low = var_lower_bound(var);
313 high = var_upper_bound(var);
314 val = var_value(var);
315 if( low > high ) {
316 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
317 FPRINTF(ASCERR,"Bounds for %s variable '",label);
318 var_write_name(sys,var,ASCERR);
319 FPRINTF(ASCERR,"' are inconsistent [%g,%g]",low,high);
320 error_reporter_end_flush();
321 err = err | 0x1;
322 }
323
324 if(low > val){
325 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
326 FPRINTF(ASCERR,"The %s variable '",label);
327 var_write_name(sys,var,ASCERR);
328 FPRINTF(ASCERR,"' was set below its lower bound.");
329 error_reporter_end_flush();
330 err = err | 0x2;
331 }else if( val > high ){
332 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
333 FPRINTF(ASCERR,"The %s variable '",label);
334 var_write_name(sys,var,ASCERR);
335 FPRINTF(ASCERR,"' was set above its upper bound.");
336 error_reporter_end_flush();
337 err = err | 0x4;
338 }
339 }
340 return err;
341 }
342
343 /*------------------------------------------------------------------------------
344 SOLVER REGISTRATION
345 */
346
347 /* rewrote this stuff to get rid of all the #ifdefs -- JP */
348
349 struct StaticSolverRegistration{
350 int is_active;
351 const char *name;
352 SlvRegistration *regfunc;
353 };
354
355 /*
356 The names here are only used to provide information in the case where
357 solver registration fails. The definitive solver names are in the slv*.c
358 files.
359 */
360 static const struct StaticSolverRegistration slv_reg[]={
361 /* {0,"SLV",&slv0_register} */
362 /* ,{0,"MINOS",&slv1_register} */
363 {HAVE_QRSLV,"QRSLV",&slv3_register}
364 /* ,{0,"CSLV",&slv4_register} */
365 /* ,{0,"LSSLV",&slv5_register} */
366 /* ,{0,"MPS",&slv6_register} */
367 /* ,{0,"NGSLV",&slv7_register} */
368 /* ,{0,"OPTSQP",&slv2_register} */
369 ,{HAVE_CONOPT,"CONOPT",&slv8_register}
370 ,{HAVE_CMSLV,"CMSLV",&slv9_register}
371 ,{HAVE_LRSLV,"LRSLV",&slv9a_register}
372 ,{0,NULL,NULL}
373 };
374
375 int SlvRegisterStandardClients(void){
376 int nclients = 0;
377 int newclient=0;
378 int error;
379 int i;
380 /* CONSOLE_DEBUG("REGISTERING STANDARD SOLVER ENGINES"); */
381 for(i=0;slv_reg[i].name!=NULL;++i){
382 if(slv_reg[i].is_active && slv_reg[i].regfunc){
383 error = slv_register_client(slv_reg[i].regfunc,NULL,NULL,&newclient);
384 if(error){
385 ERROR_REPORTER_HERE(ASC_PROG_ERR
386 ,"Unable to register solver '%s' (error %d)."
387 ,slv_reg[i].name,error
388 );
389 }else{
390 CONSOLE_DEBUG("Solver '%s' registered OK",slv_solver_name(newclient));
391 nclients++;
392 }
393 }else{
394 CONSOLE_DEBUG("Solver '%s' was not compiled.",slv_reg[i].name);
395 }
396 }
397 return nclients;
398 }
399
400 /*------------------------------------------------------------------------------
401 OUTPUT ASSIGNMENT AND PARTITIONG IN LOGICAL RELATIONS
402 */
403
404 #define MLIMDEBUG 0
405 /*
406 * Populates a matrix according to the sys solvers_dvars, solvers_logrels
407 * lists and the filters given. mtx given must be created (not null)
408 * and with order >= MAX( slv_get_num_solvers_logrels(sys),
409 * slv_get_num_solvers_dvars(sys));
410 * returns 0 if went ok.
411 */
412 int slv_make_log_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
413 dis_filter_t *dvf,logrel_filter_t *lrf)
414 {
415 #if MLIMDEBUG
416 FILE *fp;
417 #endif
418 int32 lr, lrlen,ord;
419 struct logrel_relation **lrp;
420
421 if (sys==NULL || mtx == NULL || dvf == NULL || lrf == NULL) {
422 FPRINTF(stderr,"make_log_incidence called with null\n");
423 return 1;
424 }
425 lrp = slv_get_solvers_logrel_list(sys);
426 assert(lrp!=NULL);
427 lrlen = slv_get_num_solvers_logrels(sys);
428 ord = MAX(lrlen,slv_get_num_solvers_dvars(sys));
429 if (ord > mtx_order(mtx)) {
430 FPRINTF(stderr,"make_incidence called with undersized matrix\n");
431 return 2;
432 }
433 for (lr=0; lr < lrlen; lr++) {
434 if (logrel_apply_filter(lrp[lr],lrf)) {
435 logrelman_get_incidence(lrp[lr],dvf,mtx);
436 }
437 }
438 #if MLIMDEBUG
439 fp = fopen("/tmp/mim.plot","w+");
440 if (fp !=NULL) {
441 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
442 fclose(fp);
443 }
444 #endif
445 return 0;
446 }
447
448 #if USEDCODE
449
450 /* returns 0 if successful, 1 if insufficient memory. Does not change
451 * the data in mtx. Orders the solvers_dvar list of the system to
452 * match the permutation on the given mtx. It is assumed that
453 * the org cols of mtx == dvar list position. Only the
454 * dvars in range lo to hi of the dvar list are affected and
455 * only these dvars should appear in the same org column range of
456 * the mtx. This should not be called on blocks less than 3x3.
457 */
458 static int reindex_dvars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
459 const mtx_matrix_t mtx)
460 {
461 struct dis_discrete **dvtmp, **dvp;
462 int32 c,v,vlen;
463
464 if (lo >= hi +1) return 0; /* job too small */
465 dvp = slv_get_solvers_dvar_list(sys);
466 vlen = slv_get_num_solvers_dvars(sys);
467 /* on dvtmp we DONT have the terminating null */
468 dvtmp = ASC_NEW_ARRAY(struct dis_discrete *,vlen);
469 if (dvtmp == NULL) {
470 return 1;
471 }
472 /* copy pointers to dvtmp in order desired and copy back rather than sort */
473 for (c=lo;c<=hi;c++) {
474 v = mtx_col_to_org(mtx,c);
475 dvtmp[c] = dvp[v];
476 }
477 /* copying back and re-sindexing */
478 for (c=lo;c<=hi;c++) {
479 dvp[c] = dvtmp[c];
480 dis_set_sindex(dvp[c],c);
481 }
482 ascfree(dvtmp);
483 return 0;
484 }
485
486 /* returns 0 if successful, 1 if insufficient memory. Does not change
487 * the data in mtx. Orders the solvers_logrel list of the system to
488 * match the permutation on the given mtx. It is assumed that
489 * the org rows of mtx == logrel list position. Only the
490 *logrels in range lo to hi of the logrel list are affected and
491 * only these logrels should appear in the same org row range of
492 * the input mtx. This should not be called on blocks less than 3x3.
493 */
494 static int reindex_logrels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
495 const mtx_matrix_t mtx)
496 {
497 struct logrel_relation **lrtmp, **lrp;
498 int32 c,v,rlen;
499
500 lrp = slv_get_solvers_logrel_list(sys);
501 rlen = slv_get_num_solvers_logrels(sys);
502 /* on lrtmp we DONT have the terminating null */
503 lrtmp = ASC_NEW_ARRAY(struct logrel_relation *,rlen);
504 if (lrtmp == NULL) {
505 return 1;
506 }
507 /* copy pointers to lrtmp in order desired and copy back rather than sort.
508 * do this only in the row range of interest. */
509 for (c=lo;c<=hi;c++) {
510 v = mtx_row_to_org(mtx,c);
511 #if RIDEBUG
512 FPRINTF(stderr,"Old logrel sindex(org) %d becoming sindex(cur) %d\n",v,c);
513 #endif
514 lrtmp[c] = lrp[v];
515 }
516 /* copying back and re-sindexing */
517 for (c=lo;c<=hi;c++) {
518 lrp[c] = lrtmp[c];
519 logrel_set_sindex(lrp[c],c);
520 }
521 ascfree(lrtmp);
522 return 0;
523 }
524
525 #endif /* USEDCODE */
526
527 /*
528 * returns 0 if ok, OTHERWISE if madness detected.
529 */
530 int slv_log_block_partition(slv_system_t sys)
531 {
532 #if SLBPDEBUG
533 FILE *fp;
534 #endif
535 struct logrel_relation **lrp;
536 struct dis_discrete **dvp;
537 mtx_region_t *newblocks;
538 dof_t *d;
539 int32 nrow,ncol;
540 int32 ncolpfix, nrowpun, vlenmnv;
541 mtx_matrix_t mtx;
542 int32 order, rank;
543 int32 c,len,dvlen,r,lrlen;
544 dis_filter_t dvf;
545 logrel_filter_t lrf;
546
547 lrp = slv_get_solvers_logrel_list(sys);
548 dvp = slv_get_solvers_dvar_list(sys);
549 lrlen = slv_get_num_solvers_logrels(sys);
550 dvlen = slv_get_num_solvers_dvars(sys);
551 if (lrlen ==0 || dvlen == 0) return 1;
552 order = MAX(lrlen,dvlen);
553
554 lrf.matchbits = (LOGREL_ACTIVE);
555 lrf.matchvalue = (LOGREL_ACTIVE);
556 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
557 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
558
559 /* nrow plus unincluded logrels */
560 nrowpun = slv_count_solvers_logrels(sys,&lrf);
561 /* ncol plus fixed discrete vars */
562 ncolpfix = slv_count_solvers_dvars(sys,&dvf);
563
564 dvf.matchbits = (DIS_BVAR);
565 dvf.matchvalue = 0;
566
567 /* dvlen minus non boolean vars */
568 vlenmnv = dvlen - slv_count_solvers_dvars(sys,&dvf);
569
570 lrf.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
571 lrf.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
572 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_FIXED | DIS_ACTIVE);
573 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
574
575 nrow = slv_count_solvers_logrels(sys,&lrf);
576 ncol = slv_count_solvers_dvars(sys,&dvf);
577
578 mtx = slv_get_sys_mtx(sys);
579 mtx_set_order(mtx,order);
580
581 if (slv_make_log_incidence_mtx(sys,mtx,&dvf,&lrf)) {
582 ERROR_REPORTER_HERE(ASC_PROG_ERR,"failure in creating incidence matrix.");
583 mtx_destroy(mtx);
584 return 1;
585 }
586
587 mtx_output_assign(mtx,lrlen,dvlen);
588 rank = mtx_symbolic_rank(mtx);
589 if (rank == 0 ) return 1;
590 if (rank < nrow) {
591 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank < nrow. Dependent logical relations?");
592 }
593 if (rank < ncol) {
594 if ( nrow != rank) {
595 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank < ncol and nrow != rank. Excess of columns?");
596 } else {
597 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank<ncol but nrow==rank. Degrees of freedom?");
598 }
599 }
600 if (ncol == nrow) {
601 if (ncol == rank) {
602 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"System of logical relations does not need Logic Inference.\n");
603 }
604 if (ncol != rank) {
605 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"but ncol!=rank. Rank deficient?");
606 }
607 }
608 mtx_partition(mtx);
609 /* copy the block list. there is at least one if rank >=1 as above */
610 len = mtx_number_of_blocks(mtx);
611 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
612 if (newblocks == NULL) {
613 mtx_destroy(mtx);
614 return 2;
615 }
616 for (c = 0 ; c < len; c++) {
617 mtx_block(mtx,c,&(newblocks[c]));
618 }
619 slv_set_solvers_log_blocks(sys,len,newblocks); /* give it to the system */
620 d = slv_get_log_dofdata(sys);
621 d->structural_rank = rank;
622 d->n_rows = nrow;
623 d->n_cols = ncol;
624 d->n_fixed = ncolpfix - ncol;
625 d->n_unincluded = nrowpun - nrow;
626 d->reorder.partition = 1; /* yes */
627 d->reorder.basis_selection = 0; /* none yet */
628 d->reorder.block_reordering = 0; /* none */
629
630 #if SLBPDEBUG
631 fp = fopen("/tmp/sbp1.plot","w+");
632 if (fp !=NULL) {
633 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
634 fclose(fp);
635 }
636 #endif
637
638 len = vlenmnv - 1; /* row of last possible boolean variable */
639 for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
640 r = mtx_col_to_org(mtx,c);
641 if ( dis_fixed(dvp[r]) && dis_active(dvp[r])) {
642 c++;
643 } else {
644 mtx_swap_cols(mtx,c,len);
645 len--;
646 }
647 }
648
649 /* Now, leftover columns are where we want them.
650 * That is, unassigned incidence is next to assigned and fixed on right with
651 * no columns which do not correspond to dis vars in the range 0..dvlen-1.
652 */
653 /* there aren't any nonincident dvars in the solvers dvarlist, so
654 * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
655
656 /* rows 0->lrlen-1 should be actual logrels, though, and not phantoms.
657 * At this point we should have
658 * rows - assigned
659 * - unassigned w/incidence
660 * - unincluded and included w/out free incidence and inactive mixed.
661 */
662
663 /* sort unincluded logrelations to bottom of good region */
664 /* this is needed because the idiot user may have fixed everything */
665 /* in a logrelation, but we want those rows not to be confused with */
666 /* unincluded rows.*/
667 len = lrlen - 1; /* row of last possible relation */
668 /* sort the inactive out of the unassigned/unincluded */
669 for (c=rank; len > c ; ) {
670 r = mtx_row_to_org(mtx,c);
671 if ( logrel_active(lrp[r])) {
672 c++;
673 } else {
674 mtx_swap_rows(mtx,c,len);
675 len--;
676 }
677 }
678
679 /* sort the unincluded out of the unassigned */
680 len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
681 for (c=rank; len > c ; ) {
682 r = mtx_row_to_org(mtx,c);
683 if ( logrel_included(lrp[r]) ) {
684 c++;
685 } else {
686 mtx_swap_rows(mtx,c,len);
687 len--;
688 }
689 }
690
691 #if SLBPDEBUG
692 fp = fopen("/tmp/sbp2.plot","w+");
693 if(fp !=NULL){
694 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
695 fclose(fp);
696 }
697 #endif
698
699 /* These functions are called for slv_block_partitioning, but I do not
700 * know if I need them here:
701 *
702 * now, at last we have cols in the order we want the lists to
703 * be handed to the solver. So, we need to reset the dvar sindex values
704 * and reorder the lists.
705 */
706
707 #if USEDCODE
708 if (reindex_dvars_from_mtx(sys,0,dvlen-1,mtx)) {
709 mtx_destroy(mtx);
710 return 2;
711 }
712 #endif /* USEDCODE */
713
714 /*
715 * now, at last we have rows jacobian in the order we want the lists to
716 * be handed to the solvers. So, we need to reset the rel sindex values.
717 */
718
719 #if USEDCODE
720 if (reindex_logrels_from_mtx(sys,0,lrlen-1,mtx)) {
721 mtx_destroy(mtx);
722 return 2;
723 }
724 #endif /* USEDCODE */
725
726 return 0;
727 }
728

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