/[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 1207 - (show annotations) (download) (as text)
Tue Jan 23 03:17:50 2007 UTC (17 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 21272 byte(s)
Cleaned up some #define and #includes
1 /* ASCEND modelling environment
2 Copyright (C) 1998, 2006 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_std_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_std_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,
116 int32 *var_count)
117 {
118 struct rel_relation **rp;
119 struct rel_relation **rtmp;
120 struct rel_relation *rel;
121 struct var_variable **vp;
122 struct var_variable **vtmp;
123 struct var_variable *var;
124 int32 nrow,ncol,rlen,vlen,order,rel_tmp_end;
125 int32 r,c,var_tmp_end,len;
126 var_filter_t vf;
127 rel_filter_t rf;
128
129 /* get rel and var info */
130 rp = slv_get_solvers_rel_list(sys);
131 vp = slv_get_solvers_var_list(sys);
132 rlen = slv_get_num_solvers_rels(sys);
133 vlen = slv_get_num_solvers_vars(sys);
134 if (rlen ==0 || vlen == 0) return;
135 order = MAX(rlen,vlen);
136
137 assert(var_count != NULL && rel_count != NULL);
138 *var_count = *rel_count = -1;
139
140 /* allocate temp arrays */
141 vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
142 if (vtmp == NULL) {
143 return;
144 }
145 rtmp = ASC_NEW_ARRAY(struct rel_relation *,rlen);
146 if (rtmp == NULL) {
147 ascfree(vtmp);
148 return;
149 }
150
151 /* set up filters */
152 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
153 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
154 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
155 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
156
157 /* count rows and cols */
158 nrow = slv_count_solvers_rels(sys,&rf);
159 ncol = slv_count_solvers_vars(sys,&vf);
160
161
162 /* Sort out unincluded and inactive rels to
163 * the end of the temporary relation list.
164 */
165 *rel_count = 0;
166 rel_tmp_end = rlen -1;
167 for (r=0; r < rlen; r++) {
168 rel = rp[r];
169 if (rel_apply_filter(rel,&rf)) {
170 rtmp[*rel_count] = rel;
171 (*rel_count)++;
172 } else {
173 rtmp[rel_tmp_end] = rel;
174 rel_tmp_end--;
175 }
176 }
177
178 /* sort the inactive out of the unincluded and move to end */
179 len = rlen -1;
180 for(c = rel_tmp_end + 1; len > c; ) {
181 rel = rtmp[c];
182 if (rel_active(rel)) {
183 c++;
184 } else {
185 rtmp[c] = rtmp[len];
186 rtmp[len] = rel;
187 len--;
188 }
189 }
190
191 /* Sort variable list */
192 *var_count = 0;
193 var_tmp_end = vlen - 1;
194 for (c = 0; c < vlen; c++) {
195 var = vp[c];
196 if (var_apply_filter(var,&vf)) {
197 vtmp[*var_count] = var;
198 (*var_count)++;
199 } else {
200 vtmp[var_tmp_end] = var;
201 var_tmp_end--;
202 }
203 }
204 /* sort the inactive out of the unincluded and move to end */
205 len = vlen -1;
206 for(c = var_tmp_end + 1; len > c; ) {
207 var = vtmp[c];
208 if (var_active(var) && var_active(var)) {
209 c++;
210 } else {
211 vtmp[c] = vtmp[len];
212 vtmp[len] = var;
213 len--;
214 }
215 }
216
217 /* reset solver indicies */
218 for (c = 0; c < vlen; c++) {
219 vp[c] = vtmp[c];
220 var_set_sindex(vp[c],c);
221 }
222
223 for (c = 0; c < rlen; c++) {
224 rp[c] = rtmp[c];
225 rel_set_sindex(rp[c],c);
226 }
227 ascfree(vtmp);
228 ascfree(rtmp);
229 return;
230 }
231
232 /*------------------------------------------------------------------------------
233 QUERYING AND ENFORCING BOUNDS
234 */
235
236 int slv_ensure_bounds(slv_system_t sys,int32 lo,int32 hi, FILE *mif)
237 {
238 real64 val,low,high;
239 int32 c,nchange=0;
240 struct var_variable *var, **vp;
241
242 vp = slv_get_solvers_var_list(sys);
243 if (vp==NULL) return -1;
244 for (c= lo; c <= hi; c++) {
245 var = vp[c];
246 low = var_lower_bound(var);
247 high = var_upper_bound(var);
248 val = var_value(var);
249 if( low > high ) {
250 if (mif!=NULL) {
251 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
252 FPRINTF(ASCERR,"Bounds for variable '");
253 var_write_name(sys,var,ASCERR);
254 FPRINTF(ASCERR,"' are inconsistent [%g,%g].\n",low,high);
255 FPRINTF(ASCERR,"Bounds will be swapped.\n");
256 error_reporter_end_flush();
257 }
258 var_set_upper_bound(var, low);
259 var_set_lower_bound(var, high);
260 low = var_lower_bound(var);
261 high = var_upper_bound(var);
262 nchange++;
263 }
264
265 if( low > val ) {
266 if (mif!=NULL) {
267 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
268 FPRINTF(ASCERR,"Variable '");
269 var_write_name(sys,var,ASCERR);
270 FPRINTF(ASCERR,"' was set below its lower bound.\n");
271 FPRINTF(ASCERR,"It will be moved to its lower bound.");
272 error_reporter_end_flush();
273 }
274 var_set_value(var, low);
275 } else {
276 if( val > high ) {
277 if (mif!=NULL) {
278 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
279 FPRINTF(ASCERR,"Variable '");
280 var_write_name(sys,var,ASCERR);
281 FPRINTF(ASCERR,"' was set above its upper bound.\n");
282 FPRINTF(ASCERR,"It will be moved to its upper bound.");
283 error_reporter_end_flush();
284 }
285 var_set_value(var, high);
286 nchange++;
287 }
288 }
289 }
290 return nchange;
291 }
292
293 /* return 0 on success (ie bounds are met) */
294 int slv_check_bounds(const slv_system_t sys
295 , int32 lo,int32 hi
296 , const char *label
297 ){
298 real64 val,low,high;
299 int32 c,len;
300 struct var_variable *var, **vp;
301 int err = 0;
302
303 if(label==NULL) label = "";
304 if(sys==NULL) return -1;
305 vp = slv_get_solvers_var_list(sys);
306 if(vp==NULL) return -2;
307 len = slv_get_num_solvers_vars(sys);
308 if(hi < 0)hi+= len; /* so you can use -1 to mean 'the last' */
309 if(lo < 0)lo+= len;
310 if (lo > len || hi > len || lo < 0 || hi < 0 || lo > hi) {
311 ERROR_REPORTER_HERE(ASC_PROG_ERR,"slv_check_bounds miscalled");
312 return -1;
313 }
314
315 for (c= lo; c <= hi; c++) {
316 var = vp[c];
317 low = var_lower_bound(var);
318 high = var_upper_bound(var);
319 val = var_value(var);
320 if( low > high ) {
321 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
322 FPRINTF(ASCERR,"Bounds for %s variable '",label);
323 var_write_name(sys,var,ASCERR);
324 FPRINTF(ASCERR,"' are inconsistent [%g,%g]",low,high);
325 error_reporter_end_flush();
326 err = err | 0x1;
327 }
328
329 if(low > val){
330 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
331 FPRINTF(ASCERR,"The %s variable '",label);
332 var_write_name(sys,var,ASCERR);
333 FPRINTF(ASCERR,"' was set below its lower bound.");
334 error_reporter_end_flush();
335 err = err | 0x2;
336 }else if( val > high ){
337 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
338 FPRINTF(ASCERR,"The %s variable '",label);
339 var_write_name(sys,var,ASCERR);
340 FPRINTF(ASCERR,"' was set above its upper bound.");
341 error_reporter_end_flush();
342 err = err | 0x4;
343 }
344 }
345 return err;
346 }
347
348 /*------------------------------------------------------------------------------
349 SOLVER REGISTRATION
350 */
351
352 /* rewrote this stuff to get rid of all the #ifdefs -- JP */
353
354 struct StaticSolverRegistration{
355 int is_active;
356 const char *name;
357 SlvRegistration *regfunc;
358 };
359
360 /*
361 The names here are only used to provide information in the case where
362 solver registration fails. The definitive solver names are in the slv*.c
363 files.
364 */
365 static const struct StaticSolverRegistration slv_reg[]={
366 /* {0,"SLV",&slv0_register} */
367 /* ,{0,"MINOS",&slv1_register} */
368 {HAVE_QRSLV,"QRSLV",&slv3_register}
369 /* ,{0,"CSLV",&slv4_register} */
370 /* ,{0,"LSSLV",&slv5_register} */
371 /* ,{0,"MPS",&slv6_register} */
372 /* ,{0,"NGSLV",&slv7_register} */
373 /* ,{0,"OPTSQP",&slv2_register} */
374 ,{HAVE_CONOPT,"CONOPT",&slv8_register}
375 ,{HAVE_CMSLV,"CMSLV",&slv9_register}
376 ,{HAVE_LRSLV,"LRSLV",&slv9a_register}
377 ,{0,NULL,NULL}
378 };
379
380 int SlvRegisterStandardClients(void){
381 int nclients = 0;
382 int newclient=0;
383 int error;
384 int i;
385 /* CONSOLE_DEBUG("REGISTERING STANDARD SOLVER ENGINES"); */
386 for(i=0;slv_reg[i].name!=NULL;++i){
387 if(slv_reg[i].is_active && slv_reg[i].regfunc){
388 error = slv_register_client(slv_reg[i].regfunc,NULL,NULL,&newclient);
389 if(error){
390 ERROR_REPORTER_HERE(ASC_PROG_ERR
391 ,"Unable to register solver '%s' (error %d)."
392 ,slv_reg[i].name,error
393 );
394 }else{
395 CONSOLE_DEBUG("Solver '%s' registered OK",slv_solver_name(newclient));
396 nclients++;
397 }
398 }else{
399 CONSOLE_DEBUG("Solver '%s' was not compiled.",slv_reg[i].name);
400 }
401 }
402 return nclients;
403 }
404
405 /*------------------------------------------------------------------------------
406 OUTPUT ASSIGNMENT AND PARTITIONG IN LOGICAL RELATIONS
407 */
408
409 #define MLIMDEBUG 0
410 /*
411 * Populates a matrix according to the sys solvers_dvars, solvers_logrels
412 * lists and the filters given. mtx given must be created (not null)
413 * and with order >= MAX( slv_get_num_solvers_logrels(sys),
414 * slv_get_num_solvers_dvars(sys));
415 * returns 0 if went ok.
416 */
417 int slv_std_make_log_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
418 dis_filter_t *dvf,logrel_filter_t *lrf)
419 {
420 #if MLIMDEBUG
421 FILE *fp;
422 #endif
423 int32 lr, lrlen,ord;
424 struct logrel_relation **lrp;
425
426 if (sys==NULL || mtx == NULL || dvf == NULL || lrf == NULL) {
427 FPRINTF(stderr,"make_log_incidence called with null\n");
428 return 1;
429 }
430 lrp = slv_get_solvers_logrel_list(sys);
431 assert(lrp!=NULL);
432 lrlen = slv_get_num_solvers_logrels(sys);
433 ord = MAX(lrlen,slv_get_num_solvers_dvars(sys));
434 if (ord > mtx_order(mtx)) {
435 FPRINTF(stderr,"make_incidence called with undersized matrix\n");
436 return 2;
437 }
438 for (lr=0; lr < lrlen; lr++) {
439 if (logrel_apply_filter(lrp[lr],lrf)) {
440 logrelman_get_incidence(lrp[lr],dvf,mtx);
441 }
442 }
443 #if MLIMDEBUG
444 fp = fopen("/tmp/mim.plot","w+");
445 if (fp !=NULL) {
446 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
447 fclose(fp);
448 }
449 #endif
450 return 0;
451 }
452
453 #if USEDCODE
454
455 /* returns 0 if successful, 1 if insufficient memory. Does not change
456 * the data in mtx. Orders the solvers_dvar list of the system to
457 * match the permutation on the given mtx. It is assumed that
458 * the org cols of mtx == dvar list position. Only the
459 * dvars in range lo to hi of the dvar list are affected and
460 * only these dvars should appear in the same org column range of
461 * the mtx. This should not be called on blocks less than 3x3.
462 */
463 static int reindex_dvars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
464 const mtx_matrix_t mtx)
465 {
466 struct dis_discrete **dvtmp, **dvp;
467 int32 c,v,vlen;
468
469 if (lo >= hi +1) return 0; /* job too small */
470 dvp = slv_get_solvers_dvar_list(sys);
471 vlen = slv_get_num_solvers_dvars(sys);
472 /* on dvtmp we DONT have the terminating null */
473 dvtmp = (struct dis_discrete **)
474 ascmalloc(vlen*sizeof(struct dis_discrete *));
475 if (dvtmp == NULL) {
476 return 1;
477 }
478 /* copy pointers to dvtmp in order desired and copy back rather than sort */
479 for (c=lo;c<=hi;c++) {
480 v = mtx_col_to_org(mtx,c);
481 dvtmp[c] = dvp[v];
482 }
483 /* copying back and re-sindexing */
484 for (c=lo;c<=hi;c++) {
485 dvp[c] = dvtmp[c];
486 dis_set_sindex(dvp[c],c);
487 }
488 ascfree(dvtmp);
489 return 0;
490 }
491
492 /* returns 0 if successful, 1 if insufficient memory. Does not change
493 * the data in mtx. Orders the solvers_logrel list of the system to
494 * match the permutation on the given mtx. It is assumed that
495 * the org rows of mtx == logrel list position. Only the
496 *logrels in range lo to hi of the logrel list are affected and
497 * only these logrels should appear in the same org row range of
498 * the input mtx. This should not be called on blocks less than 3x3.
499 */
500 static int reindex_logrels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
501 const mtx_matrix_t mtx)
502 {
503 struct logrel_relation **lrtmp, **lrp;
504 int32 c,v,rlen;
505
506 lrp = slv_get_solvers_logrel_list(sys);
507 rlen = slv_get_num_solvers_logrels(sys);
508 /* on lrtmp we DONT have the terminating null */
509 lrtmp = (struct logrel_relation **)
510 ascmalloc(rlen*sizeof(struct logrel_relation *));
511 if (lrtmp == NULL) {
512 return 1;
513 }
514 /* copy pointers to lrtmp in order desired and copy back rather than sort.
515 * do this only in the row range of interest. */
516 for (c=lo;c<=hi;c++) {
517 v = mtx_row_to_org(mtx,c);
518 #if RIDEBUG
519 FPRINTF(stderr,"Old logrel sindex(org) %d becoming sindex(cur) %d\n",v,c);
520 #endif
521 lrtmp[c] = lrp[v];
522 }
523 /* copying back and re-sindexing */
524 for (c=lo;c<=hi;c++) {
525 lrp[c] = lrtmp[c];
526 logrel_set_sindex(lrp[c],c);
527 }
528 ascfree(lrtmp);
529 return 0;
530 }
531
532 #endif /* USEDCODE */
533
534 /*
535 * returns 0 if ok, OTHERWISE if madness detected.
536 */
537 int slv_log_block_partition(slv_system_t sys)
538 {
539 #if SLBPDEBUG
540 FILE *fp;
541 #endif
542 struct logrel_relation **lrp;
543 struct dis_discrete **dvp;
544 mtx_region_t *newblocks;
545 dof_t *d;
546 int32 nrow,ncol;
547 int32 ncolpfix, nrowpun, vlenmnv;
548 mtx_matrix_t mtx;
549 int32 order, rank;
550 int32 c,len,dvlen,r,lrlen;
551 dis_filter_t dvf;
552 logrel_filter_t lrf;
553
554 lrp = slv_get_solvers_logrel_list(sys);
555 dvp = slv_get_solvers_dvar_list(sys);
556 lrlen = slv_get_num_solvers_logrels(sys);
557 dvlen = slv_get_num_solvers_dvars(sys);
558 if (lrlen ==0 || dvlen == 0) return 1;
559 order = MAX(lrlen,dvlen);
560
561 lrf.matchbits = (LOGREL_ACTIVE);
562 lrf.matchvalue = (LOGREL_ACTIVE);
563 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
564 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
565
566 /* nrow plus unincluded logrels */
567 nrowpun = slv_count_solvers_logrels(sys,&lrf);
568 /* ncol plus fixed discrete vars */
569 ncolpfix = slv_count_solvers_dvars(sys,&dvf);
570
571 dvf.matchbits = (DIS_BVAR);
572 dvf.matchvalue = 0;
573
574 /* dvlen minus non boolean vars */
575 vlenmnv = dvlen - slv_count_solvers_dvars(sys,&dvf);
576
577 lrf.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
578 lrf.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
579 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_FIXED | DIS_ACTIVE);
580 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
581
582 nrow = slv_count_solvers_logrels(sys,&lrf);
583 ncol = slv_count_solvers_dvars(sys,&dvf);
584
585 mtx = slv_get_sys_mtx(sys);
586 mtx_set_order(mtx,order);
587
588 if (slv_std_make_log_incidence_mtx(sys,mtx,&dvf,&lrf)) {
589 ERROR_REPORTER_HERE(ASC_PROG_ERR,"failure in creating incidence matrix.");
590 mtx_destroy(mtx);
591 return 1;
592 }
593
594 mtx_output_assign(mtx,lrlen,dvlen);
595 rank = mtx_symbolic_rank(mtx);
596 if (rank == 0 ) return 1;
597 if (rank < nrow) {
598 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank < nrow. Dependent logical relations?");
599 }
600 if (rank < ncol) {
601 if ( nrow != rank) {
602 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank < ncol and nrow != rank. Excess of columns?");
603 } else {
604 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank<ncol but nrow==rank. Degrees of freedom?");
605 }
606 }
607 if (ncol == nrow) {
608 if (ncol == rank) {
609 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"System of logical relations does not need Logic Inference.");
610 }
611 if (ncol != rank) {
612 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"but ncol!=rank. Rank deficient?");
613 }
614 }
615 mtx_partition(mtx);
616 /* copy the block list. there is at least one if rank >=1 as above */
617 len = mtx_number_of_blocks(mtx);
618 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
619 if (newblocks == NULL) {
620 mtx_destroy(mtx);
621 return 2;
622 }
623 for (c = 0 ; c < len; c++) {
624 mtx_block(mtx,c,&(newblocks[c]));
625 }
626 slv_set_solvers_log_blocks(sys,len,newblocks); /* give it to the system */
627 d = slv_get_log_dofdata(sys);
628 d->structural_rank = rank;
629 d->n_rows = nrow;
630 d->n_cols = ncol;
631 d->n_fixed = ncolpfix - ncol;
632 d->n_unincluded = nrowpun - nrow;
633 d->reorder.partition = 1; /* yes */
634 d->reorder.basis_selection = 0; /* none yet */
635 d->reorder.block_reordering = 0; /* none */
636
637 #if SLBPDEBUG
638 fp = fopen("/tmp/sbp1.plot","w+");
639 if (fp !=NULL) {
640 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
641 fclose(fp);
642 }
643 #endif
644
645 len = vlenmnv - 1; /* row of last possible boolean variable */
646 for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
647 r = mtx_col_to_org(mtx,c);
648 if ( dis_fixed(dvp[r]) && dis_active(dvp[r])) {
649 c++;
650 } else {
651 mtx_swap_cols(mtx,c,len);
652 len--;
653 }
654 }
655
656 /* Now, leftover columns are where we want them.
657 * That is, unassigned incidence is next to assigned and fixed on right with
658 * no columns which do not correspond to dis vars in the range 0..dvlen-1.
659 */
660 /* there aren't any nonincident dvars in the solvers dvarlist, so
661 * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
662
663 /* rows 0->lrlen-1 should be actual logrels, though, and not phantoms.
664 * At this point we should have
665 * rows - assigned
666 * - unassigned w/incidence
667 * - unincluded and included w/out free incidence and inactive mixed.
668 */
669
670 /* sort unincluded logrelations to bottom of good region */
671 /* this is needed because the idiot user may have fixed everything */
672 /* in a logrelation, but we want those rows not to be confused with */
673 /* unincluded rows.*/
674 len = lrlen - 1; /* row of last possible relation */
675 /* sort the inactive out of the unassigned/unincluded */
676 for (c=rank; len > c ; ) {
677 r = mtx_row_to_org(mtx,c);
678 if ( logrel_active(lrp[r])) {
679 c++;
680 } else {
681 mtx_swap_rows(mtx,c,len);
682 len--;
683 }
684 }
685
686 /* sort the unincluded out of the unassigned */
687 len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
688 for (c=rank; len > c ; ) {
689 r = mtx_row_to_org(mtx,c);
690 if ( logrel_included(lrp[r]) ) {
691 c++;
692 } else {
693 mtx_swap_rows(mtx,c,len);
694 len--;
695 }
696 }
697
698 #if SLBPDEBUG
699 fp = fopen("/tmp/sbp2.plot","w+");
700 if(fp !=NULL){
701 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
702 fclose(fp);
703 }
704 #endif
705
706 /* These functions are called for slv_block_partitioning, but I do not
707 * know if I need them here:
708 *
709 * now, at last we have cols in the order we want the lists to
710 * be handed to the solver. So, we need to reset the dvar sindex values
711 * and reorder the lists.
712 */
713
714 #if USEDCODE
715 if (reindex_dvars_from_mtx(sys,0,dvlen-1,mtx)) {
716 mtx_destroy(mtx);
717 return 2;
718 }
719 #endif /* USEDCODE */
720
721 /*
722 * now, at last we have rows jacobian in the order we want the lists to
723 * be handed to the solvers. So, we need to reset the rel sindex values.
724 */
725
726 #if USEDCODE
727 if (reindex_logrels_from_mtx(sys,0,lrlen-1,mtx)) {
728 mtx_destroy(mtx);
729 return 2;
730 }
731 #endif /* USEDCODE */
732
733 return 0;
734 }
735

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