/[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 1055 - (show annotations) (download) (as text)
Sat Jan 6 14:47:13 2007 UTC (12 years, 11 months ago) by johnpye
File MIME type: text/x-csrc
File size: 47931 byte(s)
simplified slv_check_bounds call.
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 <utilities/ascConfig.h>
26 #include <compiler/compiler.h>
27 #include <utilities/ascMalloc.h>
28 #include <general/list.h>
29 #include "mtx.h"
30 #include "slv_types.h"
31 #include "var.h"
32 #include "rel.h"
33 #include "discrete.h"
34 #include "conditional.h"
35 #include "logrel.h"
36 #include "relman.h"
37 #include "logrelman.h"
38 #include "bnd.h"
39 #include "slv_common.h"
40 #include "linsol.h"
41 #include "linsolqr.h"
42 #include "slv_client.h"
43
44 /* headers of registered clients */
45 #include "slv0.h"
46 #include "slv1.h"
47 #include "slv2.h"
48 #include "slv3.h"
49 #include "slv6.h"
50 #include "slv7.h"
51 #include "slv8.h"
52 #include "slv9.h"
53 #include "slv9a.h"
54
55 #include "model_reorder.h"
56 #include "slv_stdcalls.h"
57 #include <general/mathmacros.h>
58
59 #define KILL 0 /* deleteme code. compile old code if kill = 1 */
60 #define NEEDSTOBEDONE 0 /* indicates code/comments that are not yet ready */
61 #define USECODE 0 /* Code in good shape but not used currently */
62
63 #define MIMDEBUG 0 /* slv_std_make_incidence_mtx debugging */
64 #define RIDEBUG 0 /* reindex debugging */
65 #define SBPDEBUG 0 /* slv_block_partition_real debugging */
66 #define SLBPDEBUG 0 /* slv_log_block_partition debugging */
67
68 /* global to get around the mr header (for tear_subreorder) */
69 static
70 enum mtx_reorder_method g_blockmethod = mtx_UNKNOWN;
71
72 /*
73 Here we play some hefty Jacobian reordering games.
74
75 What we want to happen eventually is as follows:
76
77 - Get the free & incident pattern for include relations.
78 - Output-assign the Jacobian.
79 - BLT permute the Jacobian. If underspecified, fake rows
80 to make things appear square.
81 - For all leading square blocks apply reordering (kirk, etc),
82 - If trailing rectangular block, "apply clever kirkbased scheme not known."???
83
84 At present, we aren't quite so clever. We:
85
86 - Get the free & incident pattern for include relations.
87 - Output-assign the Jacobian.
88 - BLT permute the square output assigned region of the Jacobian.
89 - For all leading square blocks apply reordering (kirk, etc)
90 - Collect the block list as part of the master data structure.
91 - Set sindices for rels, vars as current rows/cols in matrix so
92 that the jacobian is 'naturally' preordered for all solvers.
93
94 Solvers are still free to reorder their own matrices any way they like.
95 It's probably a dumb idea, though.
96 */
97
98 /* see slv_stdcalls.h */
99 int slv_std_make_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
100 var_filter_t *vf,rel_filter_t *rf)
101 {
102 #if MIMDEBUG
103 FILE *fp;
104 #endif
105 int32 r, rlen,ord;
106 struct rel_relation **rp;
107
108 if (sys==NULL || mtx == NULL || vf == NULL || rf == NULL) {
109 ERROR_REPORTER_HERE(ASC_PROG_ERR,"called with null");
110 return 1;
111 }
112 rp = slv_get_solvers_rel_list(sys);
113 assert(rp!=NULL);
114 rlen = slv_get_num_solvers_rels(sys);
115 ord = MAX(rlen,slv_get_num_solvers_vars(sys));
116 if (ord > mtx_order(mtx)) {
117 ERROR_REPORTER_HERE(ASC_PROG_ERR,"undersized matrix");
118 return 2;
119 }
120 for (r=0; r < rlen; r++) {
121 if (rel_apply_filter(rp[r],rf)) {
122 relman_get_incidence(rp[r],vf,mtx);
123 }
124 }
125 #if MIMDEBUG
126 fp = fopen("/tmp/mim.plot","w+");
127 if (fp !=NULL) {
128 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
129 fclose(fp);
130 }
131 #endif
132 return 0;
133 }
134
135
136 /**
137 Orders the solvers_var list of the system to match the permutation
138 on the given mtx. Does not change the data in mtx.
139
140 @return 0 on success, 1 on out-of-memory
141
142 It is assumed that the org cols of mtx == var list position. Only the
143 vars in range lo to hi of the var list are affected and
144 only these vars should appear in the same org column range of
145 the mtx. This should not be called on blocks less than 3x3.
146 */
147 static int reindex_vars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
148 const mtx_matrix_t mtx)
149 {
150 struct var_variable **vtmp, **vp;
151 int32 c,v,vlen;
152
153 if (lo >= hi +1) return 0; /* job too small */
154 vp = slv_get_solvers_var_list(sys);
155 vlen = slv_get_num_solvers_vars(sys);
156 /* on vtmp we DONT have the terminating null */
157 vtmp = ASC_NEW_ARRAY(struct var_variable *,vlen);
158 if (vtmp == NULL) {
159 return 1;
160 }
161 /* copy pointers to vtmp in order desired and copy back rather than sort */
162 for (c=lo;c<=hi;c++) {
163 v = mtx_col_to_org(mtx,c);
164 vtmp[c] = vp[v];
165 }
166 /* copying back and re-sindexing */
167 for (c=lo;c<=hi;c++) {
168 vp[c] = vtmp[c];
169 var_set_sindex(vp[c],c);
170 }
171 ascfree(vtmp);
172 return 0;
173 }
174 /**
175 Orders the solvers_rel list of the system to match the permutation
176 on the given mtx. Does not change the data in mtx.
177
178 @return 0 on success, 1 on out-of-memory
179
180 It is assumed that the org rows of mtx == rel list position. Only the
181 rels in range lo to hi of the rel list are affected and
182 only these rels should appear in the same org row range of
183 the input mtx. This should not be called on blocks less than 3x3.
184 */
185 static int reindex_rels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
186 const mtx_matrix_t mtx)
187 {
188 struct rel_relation **rtmp, **rp;
189 int32 c,v,rlen;
190
191 rp = slv_get_solvers_rel_list(sys);
192 rlen = slv_get_num_solvers_rels(sys);
193 /* on rtmp we DONT have the terminating null */
194 rtmp = ASC_NEW_ARRAY(struct rel_relation*,rlen);
195 if (rtmp == NULL) {
196 return 1;
197 }
198 /* copy pointers to rtmp in order desired and copy back rather than sort.
199 * do this only in the row range of interest. */
200 for (c=lo;c<=hi;c++) {
201 v = mtx_row_to_org(mtx,c);
202 #if RIDEBUG
203 if(c!=v){
204 CONSOLE_DEBUG("Old rel sindex (org) %d becoming sindex (cur) %d\n",v,c);
205 }
206 #endif
207 rtmp[c] = rp[v];
208 }
209 /* copying back and re-sindexing */
210 for (c=lo;c<=hi;c++) {
211 rp[c] = rtmp[c];
212 rel_set_sindex(rp[c],c);
213 }
214 ascfree(rtmp);
215 return 0;
216 }
217
218 /**
219 Perform var and rel reordering to achieve block form.
220
221 @param uppertrianguler if non-zero, BUT form. if 0, make BLT form.
222
223 @callergraph
224 */
225 int slv_block_partition_real(slv_system_t sys,int uppertriangular){
226 #if SBPDEBUG
227 FILE *fp;
228 #endif
229 struct rel_relation **rp;
230 struct var_variable **vp;
231 mtx_region_t *newblocks;
232 dof_t *d;
233 int32 nrow,ncol;
234 int32 ncolpfix, nrowpun, vlenmnv;
235 mtx_matrix_t mtx;
236 int32 order, rank;
237 int32 c,len,vlen,r,rlen;
238 var_filter_t vf;
239 rel_filter_t rf;
240
241 /* CONSOLE_DEBUG("..."); */
242
243 rp = slv_get_solvers_rel_list(sys);
244 vp = slv_get_solvers_var_list(sys);
245 rlen = slv_get_num_solvers_rels(sys);
246 vlen = slv_get_num_solvers_vars(sys);
247 if (rlen ==0 || vlen == 0) return 1;
248 order = MAX(rlen,vlen);
249
250 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
251 rf.matchvalue = (REL_ACTIVE);
252 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
253 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
254
255 /* nrow plus unincluded rels and ineqs */
256 nrowpun = slv_count_solvers_rels(sys,&rf);
257 ncolpfix = slv_count_solvers_vars(sys,&vf); /* ncol plus fixed vars */
258
259 vf.matchbits = (VAR_SVAR);
260 vf.matchvalue = 0;
261 vlenmnv = vlen - slv_count_solvers_vars(sys,&vf); /* vlen minus non vars */
262
263 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
264 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
265 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
266 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
267
268 nrow = slv_count_solvers_rels(sys,&rf);
269 ncol = slv_count_solvers_vars(sys,&vf);
270
271 mtx = mtx_create();
272 mtx_set_order(mtx,order);
273
274 if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
275 ERROR_REPORTER_HERE(ASC_PROG_ERR,"failure in creating incidence matrix.");
276 mtx_destroy(mtx);
277 return 1;
278 }
279
280 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
281
282 mtx_output_assign(mtx,rlen,vlen);
283 rank = mtx_symbolic_rank(mtx);
284 if (rank == 0 ) return 1; /* nothing to do, eh? */
285
286 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
287
288 /* lot of whining about dof */
289 if (rank < nrow) {
290 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is row rank deficient (%d dependent equations)",
291 nrow - rank);
292 }
293 if (rank < ncol) {
294 if ( nrow != rank) {
295 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is row rank deficient with %d excess columns.",
296 ncol - rank);
297 } else {
298 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System has %d degrees of freedom.", ncol - rank);
299 }
300 }
301 if (ncol == nrow) {
302 if (ncol != rank) {
303 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is (%d) square but rank deficient.",ncol);
304 } else {
305 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"System is (%d) square.",ncol);
306 }
307 }
308 if (uppertriangular) {
309 mtx_ut_partition(mtx);
310 } else {
311 mtx_partition(mtx);
312 }
313 /* copy the block list. there is at least one if rank >=1 as above */
314 len = mtx_number_of_blocks(mtx);
315 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
316 if (newblocks == NULL) {
317 mtx_destroy(mtx);
318 return 2;
319 }
320 for (c = 0 ; c < len; c++) {
321 mtx_block(mtx,c,&(newblocks[c]));
322 }
323 slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
324 d = slv_get_dofdata(sys);
325 d->structural_rank = rank;
326 d->n_rows = nrow;
327 d->n_cols = ncol;
328
329 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
330
331 /* the next two lines assume inequalities are un-included. */
332 #define MINE 1
333 #if MINE
334 d->n_fixed = ncolpfix - ncol;
335 d->n_unincluded = nrowpun - nrow;
336 #else
337 d->n_fixed = vlen - ncol;
338 d->n_unincluded = rlen - nrow;
339 #endif
340 d->reorder.partition = 1; /* yes */
341 d->reorder.basis_selection = 0; /* none yet */
342 d->reorder.block_reordering = 0; /* none */
343
344 #if SBPDEBUG
345 fp = fopen("/tmp/sbp1.plot","w+");
346 if (fp !=NULL) {
347 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
348 fclose(fp);
349 }
350 #endif
351
352 #if MINE
353 len = vlenmnv - 1; /* row of last possible solver variable */
354 for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
355 r = mtx_col_to_org(mtx,c);
356 if ( var_fixed(vp[r]) && var_active(vp[r])) {
357 c++;
358 } else {
359 mtx_swap_cols(mtx,c,len);
360 len--;
361 }
362 }
363 #endif
364
365 /* Now, leftover columns are where we want them.
366 * That is, unassigned incidence is next to assigned and fixed on right with
367 * no columns which do not correspond to variables in the range 0..vlen-1.
368 */
369 /* there aren't any nonincident vars in the solvers varlist, so
370 * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
371
372 /* rows 0->rlen-1 should be actual rels, though, and not phantoms.
373 * At this point we should have
374 * rows - assigned
375 * - unassigned w/incidence
376 * - unincluded and included w/out free incidence and inactive mixed.
377 */
378 /* sort unincluded relations to bottom of good region */
379 /* this is needed because the idiot user may have fixed everything */
380 /* in a relation, but we want those rows not to be confused with */
381 /* unincluded rows. we're sort of dropping inequalities on the floor.*/
382 len = rlen - 1; /* row of last possible relation */
383 /* sort the inactive out of the unassigned/unincluded */
384 for (c=rank; len > c ; ) {
385 r = mtx_row_to_org(mtx,c);
386 if ( rel_active(rp[r])) {
387 c++;
388 } else {
389 mtx_swap_rows(mtx,c,len);
390 len--;
391 }
392 }
393
394 #if MINE
395 /* sort the unincluded out of the unassigned */
396 len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
397 for (c=rank; len > c ; ) {
398 r = mtx_row_to_org(mtx,c);
399 if ( rel_included(rp[r]) ) {
400 c++;
401 } else {
402 mtx_swap_rows(mtx,c,len);
403 len--;
404 }
405 }
406 #endif
407
408
409 #if SBPDEBUG
410 fp = fopen("/tmp/sbp2.plot","w+");
411 if (fp !=NULL) {
412 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
413 fclose(fp);
414 }
415 #endif
416 /* now, at last we have cols jacobian in the order we want the lists to
417 * be handed to the solvers. So, we need to reset the var sindex values
418 * and reorder the lists.
419 */
420 if (reindex_vars_from_mtx(sys,0,vlen-1,mtx)) {
421 mtx_destroy(mtx);
422 return 2;
423 }
424 /* now, at last we have rows jacobian in the order we want the lists to
425 * be handed to the solvers. So, we need to reset the rel sindex values.
426 */
427 if (reindex_rels_from_mtx(sys,0,rlen-1,mtx)) {
428 mtx_destroy(mtx);
429 return 2;
430 }
431
432 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
433
434 mtx_destroy(mtx);
435 return 0;
436 }
437
438
439 #if 0 /* code not currently used */
440 #\ifdef STATIC_HARWELL /* added the backslash so that syntax highlighting behaves */
441 /* note that you can't statically link to Harwell routines under the terms of the GPL */
442 extern void mc21b();
443 extern void mc13emod();
444 /* returns 0 if ok, OTHERWISE if madness detected.
445 * doesn't grok inequalities.
446 * CURRENTLY DOESN'T DO WELL WHEN NCOL<NROW
447 */
448 #define SBPDEBUG 0
449 int slv_block_partition_harwell(slv_system_t sys)
450 {
451 #\if SBPDEBUG
452 FILE *fp;
453 #\endif
454 struct rel_relation **rp;
455 struct rel_relation **rtmp;
456 struct rel_relation *rel;
457 struct var_variable **vp;
458 struct var_variable **vtmp;
459 struct var_variable *var;
460 const struct var_variable **list;
461 mtx_region_t *newblocks;
462 dof_t *d;
463 int32 nrow,ncol,size;
464 int32 order, rank;
465 int32 c,len,vlen,r,v,rlen,col,row;
466 int32 licn,rel_tmp_end,rel_count,var_perm_count,var_tmp_end,row_perm_count;
467 int32 col_incidence,numnz,row_high,col_high,dummy_ptr,row_ptr,num;
468 int32 *dummy_array_ptr;
469 var_filter_t vf;
470 rel_filter_t rf;
471 int32 *icn, *ip, *lenr, *iperm, *iw1, *iw2, *iw3, *arp, *ipnew, *ib;
472 int32 *row_perm;
473
474 /* get rel and var info */
475 rp = slv_get_solvers_rel_list(sys);
476 vp = slv_get_solvers_var_list(sys);
477 rlen = slv_get_num_solvers_rels(sys);
478 vlen = slv_get_num_solvers_vars(sys);
479 if (rlen ==0 || vlen == 0) return 1;
480 order = MAX(rlen,vlen);
481
482 /* allocate temp arrays */
483 vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
484 var = (struct var_variable *)ascmalloc(sizeof(struct var_variable *));
485 if (vtmp == NULL || var ==NULL) {
486 return 1;
487 }
488 rtmp = ASC_NEW_ARRAY(struct rel_relation *,rlen);
489 rel = ASC_NEW(struct rel_relation);
490 if (rtmp == NULL || rel ==NULL) {
491 return 1;
492 }
493
494 /* set up filters */
495 rf.matchbits = (REL_INCLUDED | REL_EQUALITY);
496 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY);
497 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
498 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR);
499
500 /* count rows and cols */
501 nrow = slv_count_solvers_rels(sys,&rf);
502 ncol = slv_count_solvers_vars(sys,&vf);
503
504
505 /* count incidence for equality relations that are included
506 and active. Sorts out unincluded and inactive rels to
507 the end of the temporary relation list */
508 licn = 0;
509 rel_tmp_end = rlen -1;
510 for (r=0; r < rlen; r++) {
511 rel = rp[r];
512 if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
513 len = rel_n_incidences(rel);
514 list = rel_incidence_list(rel);
515 for (c=0; c < len; c++) {
516 if( var_apply_filter(list[c],&vf) ) {
517 licn++;
518 }
519 }
520 } else {
521 rtmp[rel_tmp_end] = rp[r];
522 rel_tmp_end--;
523 }
524 }
525
526 /* sort the inactive out of the unincluded and move to end */
527 len = rlen -1;
528 for(c = rel_tmp_end + 1; len > c; ) {
529 rel = rtmp[c];
530 if (rel_active(rel)) {
531 c++;
532 } else {
533 rtmp[c] = rtmp[len];
534 rtmp[len] = rel;
535 len--;
536 }
537 }
538
539 /* Sort variable list */
540 var_perm_count = 0;
541 var_tmp_end = vlen - 1;
542 for (c = 0; c < vlen; c++) {
543 var = vp[c];
544 if (var_apply_filter(var,&vf)) {
545 vtmp[var_perm_count] = var;
546 var_perm_count++;
547 } else {
548 vtmp[var_tmp_end] = var;
549 var_tmp_end--;
550 }
551 }
552 /* sort the inactive out of the unincluded and move to end */
553 len = vlen -1;
554 for(c = var_tmp_end + 1; len > c; ) {
555 var = vtmp[c];
556 if (var_active(var) && var_active(var)) {
557 c++;
558 } else {
559 vtmp[c] = vtmp[len];
560 vtmp[len] = var;
561 len--;
562 }
563 }
564
565 /* reset solver indicies to make life easier */
566 for (c = 0; c < vlen; c++) {
567 vp[c] = vtmp[c];
568 var_set_sindex(vp[c],c);
569 }
570
571 size = MAX(nrow,ncol);
572 /* Create vectors for fortran calls */
573 icn = ASC_NEW_ARRAY(int32,licn);
574 ip = ASC_NEW_ARRAY(int32,size);
575 lenr = ASC_NEW_ARRAY(int32,size);
576 iperm = ASC_NEW_ARRAY(int32,size);
577 iw1 = ASC_NEW_ARRAY(int32,size);
578 iw2 = ASC_NEW_ARRAY(int32,size);
579 iw3 = ASC_NEW_ARRAY(int32,size);
580 arp = ASC_NEW_ARRAY(int32,size);
581 ipnew = ASC_NEW_ARRAY(int32,size);
582 ib = ASC_NEW_ARRAY(int32,size);
583 row_perm = ASC_NEW_ARRAY(int32,size);
584
585
586 /* Fill incidence vectors and place included relations w/o
587 * incidence in temporary relation list before the unincluded
588 * and inactive relations
589 */
590 col_incidence = 0;
591 row_perm_count = 0;
592 for (r=0; r < rlen; r++) {
593 rel = rp[r];
594 if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
595 len = rel_n_incidences(rel);
596 if (len > 0) {
597 ip[row_perm_count] = col_incidence + 1; /*FORTRAN*/
598 lenr[row_perm_count] = len;
599 row_perm[row_perm_count] = r;
600 row_perm_count++;
601 list = rel_incidence_list(rel);
602 for (c=0; c < len; c++) {
603 if( var_apply_filter(list[c],&vf) ) {
604 col = var_sindex(list[c]);
605 icn[col_incidence] = col + 1; /*FORTRAN*/
606 col_incidence++;
607 }
608 }
609 } else {
610 rtmp[rel_tmp_end] = rel;
611 rel_tmp_end--;
612 }
613 }
614 }
615
616 licn = col_incidence;
617 order = row_perm_count;
618 mc21b(&order,icn,&licn,ip,lenr,iperm,&numnz,iw1,arp,iw2,iw3);
619 if (order == numnz){
620 FPRINTF(stderr,"they are equal\n");
621 } else if(order == numnz + 1){
622 FPRINTF(stderr,"ADD ONE\n");
623 } else {
624 FPRINTF(stderr,"no relationship\n");
625 }
626
627 if (rank == 0 ) {
628 return 1; /* nothing to do, eh? */
629 }
630 row_high = nrow - 1;
631 col_high = ncol -1;
632
633 /* lot of whining about dof */
634 if (rank < nrow) {
635 FPRINTF(stderr,"System is row rank deficient (dependent equations)\n");
636 row_high = rank - 1;
637 /* KHACK: have found that this is executed erroneously */
638 }
639 if (rank < ncol) {
640 if ( nrow != rank) {
641 FPRINTF(stderr,"System is row rank deficient with excess columns.\n");
642 } else {
643 FPRINTF(stderr,"System has degrees of freedom.\n");
644 }
645 }
646 if (ncol == nrow) {
647 FPRINTF(stderr,"System is (%d) square",ncol);
648 if (ncol != rank) {
649 FPRINTF(stderr,"but rank deficient.\n");
650 } else {
651 FPRINTF(stderr,".\n");
652 }
653 }
654
655 dummy_array_ptr = arp;
656 arp = lenr;
657 lenr = dummy_array_ptr;
658
659 /* Fix row pointers for assigned relations. Put unassigned
660 onto temporary relation list (before included w/o incidence)
661 and put dummy dense rows into ipnew */
662 dummy_ptr = licn + 1; /* FORTRAN */
663 for (row = 0; row < order; row++) {
664 row_ptr = iperm[row] - 1; /* FORTRAN */
665 if (row_ptr != -1) { /* indicates unassigned row */
666 ipnew[row_ptr] = ip[row];
667 lenr[row_ptr] = arp[row];
668 } else {
669 ipnew[row_ptr] = dummy_ptr;
670 dummy_ptr += ncol;
671 lenr[row_ptr] = ncol;
672 /* unassigned before fixed*/
673 rtmp[rel_tmp_end] = rp[row_perm[row_ptr]];
674 rel_tmp_end--;
675 }
676 }
677 for (row = order; row < ncol; row++) {
678 ipnew[row] = dummy_ptr;
679 dummy_ptr += ncol;
680 lenr[row] = ncol;
681 }
682
683 mc13emod(&size,icn,&licn,ipnew,lenr,arp,ib,&num,iw1,iw2,iw3);
684
685
686 /* copy the block list. there is at least one if rank >=1 as above */
687 len = num;
688 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
689 if (newblocks == NULL) {
690 return 2;
691 }
692 /* set up locations of block corners with the row and column
693 orderings which will be set soon */
694 for (c = 0 ; c < len; c++) {
695 newblocks[c].row.low = ib[c]-1;
696 newblocks[c].col.low = newblocks[c].row.low;
697 }
698 for (c = 0 ; c < len -1; c++) {
699 newblocks[c].row.high = newblocks[c+1].row.low -1;
700 newblocks[c].col.high = newblocks[c].row.high;
701 }
702 newblocks[len - 1].row.high = row_high;
703 newblocks[len - 1].col.high = col_high;
704
705 slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
706 d = slv_get_dofdata(sys);
707 d->structural_rank = rank;
708 d->n_rows = nrow;
709 d->n_cols = ncol;
710 /* the next two lines assume only free and fixed incident solvervar
711 * appear in the solvers var list and inequalities are unincluded.
712 */
713 d->n_fixed = vlen - ncol;
714 d->n_unincluded = rlen - nrow;
715 d->reorder.partition = 1; /* yes */
716 d->reorder.basis_selection = 0; /* none yet */
717 d->reorder.block_reordering = 0; /* none */
718
719 for (c = 0; c < ncol; c++) {
720 v = arp[c] -1; /* FORTRAN */
721 vtmp[c] = vp[v];
722 }
723 for (c = 0; c < vlen; c++) {
724 vp[c] = vtmp[c];
725 var_set_sindex(vp[c],c);
726 }
727
728 rel_count = 0;
729 for (c = 0; c < size; c++) {
730 r = arp[c] - 1;
731 if (r < order){ /* weed out fake rows */
732 r = iperm[r] - 1;
733 if (r != -1) { /* unassigned already in rtmp */
734 rtmp[rel_count] = rp[row_perm[r]];
735 rel_count++;
736 }
737 }
738 }
739 for (c = 0; c < rlen; c++) {
740 rp[c] = rtmp[c];
741 rel_set_sindex(rp[c],c);
742 }
743 ascfree(vtmp);
744 ascfree(rtmp);
745 ascfree(icn);
746 ascfree(ip);
747 ascfree(lenr);
748 ascfree(iperm);
749 ascfree(iw1);
750 ascfree(iw2);
751 ascfree(iw3);
752 ascfree(arp);
753 ascfree(ipnew);
754 ascfree(ib);
755 return 0;
756 }
757 #\endif /*STATIC_HARWELL*/
758 #endif /* 0 */
759
760
761 void slv_sort_rels_and_vars(slv_system_t sys,
762 int32 *rel_count,
763 int32 *var_count)
764 {
765 struct rel_relation **rp;
766 struct rel_relation **rtmp;
767 struct rel_relation *rel;
768 struct var_variable **vp;
769 struct var_variable **vtmp;
770 struct var_variable *var;
771 int32 nrow,ncol,rlen,vlen,order,rel_tmp_end;
772 int32 r,c,var_tmp_end,len;
773 var_filter_t vf;
774 rel_filter_t rf;
775
776 /* get rel and var info */
777 rp = slv_get_solvers_rel_list(sys);
778 vp = slv_get_solvers_var_list(sys);
779 rlen = slv_get_num_solvers_rels(sys);
780 vlen = slv_get_num_solvers_vars(sys);
781 if (rlen ==0 || vlen == 0) return;
782 order = MAX(rlen,vlen);
783
784 assert(var_count != NULL && rel_count != NULL);
785 *var_count = *rel_count = -1;
786
787 /* allocate temp arrays */
788 vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
789 if (vtmp == NULL) {
790 return;
791 }
792 rtmp = ASC_NEW_ARRAY(struct rel_relation *,rlen);
793 if (rtmp == NULL) {
794 ascfree(vtmp);
795 return;
796 }
797
798 /* set up filters */
799 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
800 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
801 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
802 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
803
804 /* count rows and cols */
805 nrow = slv_count_solvers_rels(sys,&rf);
806 ncol = slv_count_solvers_vars(sys,&vf);
807
808
809 /* Sort out unincluded and inactive rels to
810 * the end of the temporary relation list.
811 */
812 *rel_count = 0;
813 rel_tmp_end = rlen -1;
814 for (r=0; r < rlen; r++) {
815 rel = rp[r];
816 if (rel_apply_filter(rel,&rf)) {
817 rtmp[*rel_count] = rel;
818 (*rel_count)++;
819 } else {
820 rtmp[rel_tmp_end] = rel;
821 rel_tmp_end--;
822 }
823 }
824
825 /* sort the inactive out of the unincluded and move to end */
826 len = rlen -1;
827 for(c = rel_tmp_end + 1; len > c; ) {
828 rel = rtmp[c];
829 if (rel_active(rel)) {
830 c++;
831 } else {
832 rtmp[c] = rtmp[len];
833 rtmp[len] = rel;
834 len--;
835 }
836 }
837
838 /* Sort variable list */
839 *var_count = 0;
840 var_tmp_end = vlen - 1;
841 for (c = 0; c < vlen; c++) {
842 var = vp[c];
843 if (var_apply_filter(var,&vf)) {
844 vtmp[*var_count] = var;
845 (*var_count)++;
846 } else {
847 vtmp[var_tmp_end] = var;
848 var_tmp_end--;
849 }
850 }
851 /* sort the inactive out of the unincluded and move to end */
852 len = vlen -1;
853 for(c = var_tmp_end + 1; len > c; ) {
854 var = vtmp[c];
855 if (var_active(var) && var_active(var)) {
856 c++;
857 } else {
858 vtmp[c] = vtmp[len];
859 vtmp[len] = var;
860 len--;
861 }
862 }
863
864 /* reset solver indicies */
865 for (c = 0; c < vlen; c++) {
866 vp[c] = vtmp[c];
867 var_set_sindex(vp[c],c);
868 }
869
870 for (c = 0; c < rlen; c++) {
871 rp[c] = rtmp[c];
872 rel_set_sindex(rp[c],c);
873 }
874 ascfree(vtmp);
875 ascfree(rtmp);
876 return;
877 }
878
879
880
881
882
883 int slv_block_unify(slv_system_t sys)
884 {
885 dof_t *d;
886 const mtx_block_t *mbt;
887 mtx_region_t *newblocks;
888
889 if (sys==NULL) return 1;
890
891 d = slv_get_dofdata(sys);
892 mbt = slv_get_solvers_blocks(sys);
893 assert(d!=NULL && mbt!=NULL);
894 if (d->structural_rank && mbt->nblocks >1) {
895 newblocks = ASC_NEW(mtx_region_t);
896 if (newblocks == NULL) return 2;
897 newblocks->row.low = newblocks->col.low = 0;
898 newblocks->row.high = d->structural_rank - 1;
899 newblocks->col.high = d->n_cols - 1;
900 if ( newblocks->col.high < 0 || newblocks->row.high < 0) {
901 ascfree(newblocks);
902 return 1;
903 }
904 slv_set_solvers_blocks(sys,1,newblocks);
905 }
906 return 0;
907 }
908
909 int slv_set_up_block(slv_system_t sys,int bnum)
910 {
911 struct rel_relation **rp;
912 struct var_variable **vp;
913 mtx_region_t reg;
914 const mtx_block_t *b;
915 int32 c,vlen,rlen;
916 var_filter_t vf;
917 rel_filter_t rf;
918
919 if (sys==NULL) return 1;
920 rlen = slv_get_num_solvers_rels(sys);
921 vlen = slv_get_num_solvers_vars(sys);
922 if (rlen ==0 || vlen == 0) return 1;
923
924 rp = slv_get_solvers_rel_list(sys);
925 vp = slv_get_solvers_var_list(sys);
926 assert(rp!=NULL);
927 assert(vp!=NULL);
928
929 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
930 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
931 vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
932 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
933
934 b = slv_get_solvers_blocks(sys);
935 assert(b!=NULL);
936 if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
937 reg = b->block[bnum];
938 for (c=reg.col.low; c<=reg.col.high; c++) {
939 var_set_in_block(vp[c],1);
940 }
941 for (c=reg.row.low; c<=reg.row.high; c++) {
942 rel_set_in_block(rp[c],1);
943 }
944 return 0;
945 }
946
947 /* returns 0 if ok, OTHERWISE if madness detected.
948 */
949 int slv_spk1_reorder_block(slv_system_t sys,int bnum,int transpose)
950 {
951 struct rel_relation **rp;
952 struct var_variable **vp;
953 mtx_region_t reg;
954 const mtx_block_t *b;
955 mtx_matrix_t mtx;
956 mtx_coord_t coord;
957 int32 c,vlen,rlen;
958 var_filter_t vf;
959 rel_filter_t rf;
960 dof_t *d;
961
962 if (sys==NULL) return 1;
963 rlen = slv_get_num_solvers_rels(sys);
964 vlen = slv_get_num_solvers_vars(sys);
965 if (rlen ==0 || vlen == 0) return 1;
966
967 rp = slv_get_solvers_rel_list(sys);
968 vp = slv_get_solvers_var_list(sys);
969 assert(rp!=NULL);
970 assert(vp!=NULL);
971 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
972 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
973 vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
974 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
975
976 mtx = mtx_create();
977 mtx_set_order(mtx,MAX(rlen,vlen));
978 b = slv_get_solvers_blocks(sys);
979 assert(b!=NULL);
980 if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
981 reg = b->block[bnum];
982 for (c=reg.col.low; c<=reg.col.high; c++) {
983 var_set_in_block(vp[c],1);
984 }
985 for (c=reg.row.low; c<=reg.row.high; c++) {
986 rel_set_in_block(rp[c],1);
987 }
988 if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
989 return 1; /* must be square */
990 /* could also enforce minsize 3x3, but someone my call with a
991 * partitionable region, so don't want to.
992 */
993 }
994
995 if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
996 FPRINTF(stderr,
997 "slv_spk1_reorder_block: failure in creating incidence matrix.\n");
998 mtx_destroy(mtx);
999 return 1;
1000 }
1001 /* verify that block has no empty columns, though not checking diagonal */
1002 for (c = reg.row.low; c <= reg.row.high; c++) {
1003 coord.col = mtx_FIRST;
1004 coord.row = c;
1005 if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
1006 mtx_destroy(mtx);
1007 FPRINTF(stderr, "slv_spk1_reorder_block: empty row (%d) found.\n",c);
1008 return 1;
1009 }
1010 coord.row = mtx_FIRST;
1011 coord.col = c;
1012 if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
1013 FPRINTF(stderr, "slv_spk1_reorder_block: empty col (%d) found.\n",c);
1014 mtx_destroy(mtx);
1015 return 1;
1016 }
1017 }
1018 if (transpose) {
1019 mtx_reorder(mtx,&reg,mtx_TSPK1);
1020 } else {
1021 mtx_reorder(mtx,&reg,mtx_SPK1);
1022 }
1023 if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
1024 mtx_destroy(mtx);
1025 return 2;
1026 }
1027 if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
1028 mtx_destroy(mtx);
1029 return 2;
1030 }
1031 d = slv_get_dofdata(sys);
1032 d->reorder.block_reordering = 1; /* spk1 */
1033
1034 mtx_destroy(mtx);
1035 return 0;
1036 }
1037
1038 /*
1039 A function to be called on the leftover diagonal blocks of mr_bisect.
1040 This function should be thoroughly investigated, which it has not been.
1041 */
1042 static int tear_subreorder(slv_system_t server,
1043 mtx_matrix_t mtx, mtx_region_t *reg)
1044 {
1045 assert(mtx != NULL && reg !=NULL);
1046 (void)server;
1047 if (g_blockmethod==mtx_UNKNOWN) {
1048 mtx_reorder(mtx,reg,mtx_SPK1);
1049 } else {
1050 mtx_reorder(mtx,reg,g_blockmethod);
1051 }
1052 return 0;
1053 }
1054
1055 int slv_tear_drop_reorder_block(slv_system_t sys, int32 bnum,
1056 int32 cutoff, int two,
1057 enum mtx_reorder_method blockmethod
1058 ){
1059 struct rel_relation **rp;
1060 struct var_variable **vp;
1061 const mtx_block_t *b;
1062 dof_t *d;
1063 mr_reorder_t *mrsys;
1064 mtx_matrix_t mtx;
1065 mtx_region_t reg;
1066 mtx_coord_t coord;
1067 int32 c,vlen,rlen,modcount;
1068 var_filter_t vf;
1069 rel_filter_t rf;
1070
1071 if (sys==NULL) return 1;
1072 rlen = slv_get_num_solvers_rels(sys);
1073 vlen = slv_get_num_solvers_vars(sys);
1074 if (rlen ==0 || vlen == 0) return 1;
1075
1076 rp = slv_get_solvers_rel_list(sys);
1077 vp = slv_get_solvers_var_list(sys);
1078 assert(rp!=NULL);
1079 assert(vp!=NULL);
1080 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
1081 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK| REL_ACTIVE);
1082 vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
1083 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
1084
1085 mtx = mtx_create();
1086 mtx_set_order(mtx,MAX(rlen,vlen));
1087 b = slv_get_solvers_blocks(sys);
1088 assert(b!=NULL); /* probably shouldn't be an assert here ... */
1089 if (bnum <0 || bnum >= b->nblocks || b->block == NULL) {
1090 mtx_destroy(mtx);
1091 return 1;
1092 }
1093 reg = b->block[bnum];
1094 /* do the flag setting even if we don't do the reorder */
1095 for (c=reg.col.low; c<=reg.col.high; c++) {
1096 var_set_in_block(vp[c],1);
1097 }
1098 for (c=reg.row.low; c<=reg.row.high; c++) {
1099 rel_set_in_block(rp[c],1);
1100 }
1101 if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
1102 mtx_destroy(mtx);
1103 return 1; /* must be square */
1104 }
1105 if (reg.row.high - reg.row.low < 3) {
1106 mtx_destroy(mtx);
1107 return 0; /* must be 3x3 or bigger to have any effect */
1108 }
1109
1110 if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
1111 FPRINTF(stderr,
1112 "slv_tear_drop_reorder_block: failure in creating incidence matrix.\n");
1113 mtx_destroy(mtx);
1114 return 1;
1115 }
1116 /* verify that block has no empty columns, though not checking diagonal */
1117 for (c = reg.row.low; c <= reg.row.high; c++) {
1118 coord.col = mtx_FIRST;
1119 coord.row = c;
1120 if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
1121 mtx_destroy(mtx);
1122 FPRINTF(stderr, "slv_tear_drop_reorder_block: empty row (%d) found.\n",c);
1123 return 1;
1124 }
1125 coord.row = mtx_FIRST;
1126 coord.col = c;
1127 if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
1128 FPRINTF(stderr, "slv_tear_drop_reorder_block: empty col (%d) found.\n",c);
1129 mtx_destroy(mtx);
1130 return 1;
1131 }
1132 }
1133
1134 modcount = slv_get_num_models(sys);
1135 mrsys = mr_reorder_create(sys,mtx,modcount);
1136 if (mrsys==NULL) return 1;
1137 mrsys->cutoff = cutoff;
1138 g_blockmethod = blockmethod;
1139 if (!two) {
1140 mr_bisect_partition(mrsys,&reg,1,tear_subreorder);
1141 } else {
1142 mr_bisect_partition2(mrsys,&reg,1,tear_subreorder);
1143 }
1144 #if 1
1145 FPRINTF(stderr,"Model-based reordering: (tear-style = %d)\n",two);
1146 FPRINTF(stderr,"Min. tearable size:\t%d\n",mrsys->cutoff);
1147 FPRINTF(stderr,"Tears:\t\t%d\n",mrsys->ntears);
1148 FPRINTF(stderr,"Partitionings:\t%d\n",mrsys->nblts);
1149 #endif
1150 reg = b->block[bnum]; /* bisect likely munged reg */
1151 if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
1152 mtx_destroy(mtx);
1153 mr_reorder_destroy(mrsys);
1154 return 2;
1155 }
1156 if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
1157 mtx_destroy(mtx);
1158 mr_reorder_destroy(mrsys);
1159 return 2;
1160 }
1161 d = slv_get_dofdata(sys);
1162 d->reorder.block_reordering = 2; /* tear_drop_baa */
1163
1164 mtx_destroy(mtx);
1165 mr_reorder_destroy(mrsys);
1166 return 0;
1167 }
1168
1169 int slv_ensure_bounds(slv_system_t sys,int32 lo,int32 hi, FILE *mif)
1170 {
1171 real64 val,low,high;
1172 int32 c,nchange=0;
1173 struct var_variable *var, **vp;
1174
1175 vp = slv_get_solvers_var_list(sys);
1176 if (vp==NULL) return -1;
1177 for (c= lo; c <= hi; c++) {
1178 var = vp[c];
1179 low = var_lower_bound(var);
1180 high = var_upper_bound(var);
1181 val = var_value(var);
1182 if( low > high ) {
1183 if (mif!=NULL) {
1184 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
1185 FPRINTF(ASCERR,"Bounds for variable '");
1186 var_write_name(sys,var,ASCERR);
1187 FPRINTF(ASCERR,"' are inconsistent [%g,%g].\n",low,high);
1188 FPRINTF(ASCERR,"Bounds will be swapped.\n");
1189 error_reporter_end_flush();
1190 }
1191 var_set_upper_bound(var, low);
1192 var_set_lower_bound(var, high);
1193 low = var_lower_bound(var);
1194 high = var_upper_bound(var);
1195 nchange++;
1196 }
1197
1198 if( low > val ) {
1199 if (mif!=NULL) {
1200 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
1201 FPRINTF(ASCERR,"Variable '");
1202 var_write_name(sys,var,ASCERR);
1203 FPRINTF(ASCERR,"' was set below its lower bound.\n");
1204 FPRINTF(ASCERR,"It will be moved to its lower bound.");
1205 error_reporter_end_flush();
1206 }
1207 var_set_value(var, low);
1208 } else {
1209 if( val > high ) {
1210 if (mif!=NULL) {
1211 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
1212 FPRINTF(ASCERR,"Variable '");
1213 var_write_name(sys,var,ASCERR);
1214 FPRINTF(ASCERR,"' was set above its upper bound.\n");
1215 FPRINTF(ASCERR,"It will be moved to its upper bound.");
1216 error_reporter_end_flush();
1217 }
1218 var_set_value(var, high);
1219 nchange++;
1220 }
1221 }
1222 }
1223 return nchange;
1224 }
1225
1226 /* return 0 on success (ie bounds are met) */
1227 int slv_check_bounds(const slv_system_t sys
1228 , int32 lo,int32 hi
1229 , const char *label
1230 ){
1231 real64 val,low,high;
1232 int32 c,len;
1233 struct var_variable *var, **vp;
1234 int err = 0;
1235
1236 if(label==NULL) label = "";
1237 if(sys==NULL) return -1;
1238 vp = slv_get_solvers_var_list(sys);
1239 if(vp==NULL) return -2;
1240 len = slv_get_num_solvers_vars(sys);
1241 if(hi < 0)hi+= len; /* so you can use -1 to mean 'the last' */
1242 if(lo < 0)lo+= len;
1243 if (lo > len || hi > len || lo < 0 || hi < 0 || lo > hi) {
1244 ERROR_REPORTER_HERE(ASC_PROG_ERR,"slv_check_bounds miscalled");
1245 return -1;
1246 }
1247
1248 for (c= lo; c <= hi; c++) {
1249 var = vp[c];
1250 low = var_lower_bound(var);
1251 high = var_upper_bound(var);
1252 val = var_value(var);
1253 if( low > high ) {
1254 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
1255 FPRINTF(ASCERR,"Bounds for %s variable '",label);
1256 var_write_name(sys,var,ASCERR);
1257 FPRINTF(ASCERR,"' are inconsistent [%g,%g]",low,high);
1258 error_reporter_end_flush();
1259 err = err | 0x1;
1260 }
1261
1262 if(low > val){
1263 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
1264 FPRINTF(ASCERR,"The %s variable '",label);
1265 var_write_name(sys,var,ASCERR);
1266 FPRINTF(ASCERR,"' was set below its lower bound.");
1267 error_reporter_end_flush();
1268 err = err | 0x2;
1269 }else if( val > high ){
1270 ERROR_REPORTER_START_NOLINE(ASC_USER_ERROR);
1271 FPRINTF(ASCERR,"The %s variable '",label);
1272 var_write_name(sys,var,ASCERR);
1273 FPRINTF(ASCERR,"' was set above its upper bound.");
1274 error_reporter_end_flush();
1275 err = err | 0x4;
1276 }
1277 }
1278 return err;
1279 }
1280
1281 /*------------------------------------------------------------------------------
1282 SOLVER REGISTRATION
1283 */
1284
1285 /* rewrote this stuff to get rid of all the #ifdefs -- JP */
1286
1287 struct StaticSolverRegistration{
1288 int is_active;
1289 const char *name;
1290 SlvRegistration *regfunc;
1291 };
1292
1293 /*
1294 The names here are only used to provide information in the case where
1295 solver registration fails. The definitive solver names are in the slv*.c
1296 files.
1297 */
1298 static const struct StaticSolverRegistration slv_reg[]={
1299 /* {0,"SLV",&slv0_register} */
1300 /* ,{0,"MINOS",&slv1_register} */
1301 {HAVE_QRSLV,"QRSLV",&slv3_register}
1302 /* ,{0,"CSLV",&slv4_register} */
1303 /* ,{0,"LSSLV",&slv5_register} */
1304 /* ,{0,"MPS",&slv6_register} */
1305 /* ,{0,"NGSLV",&slv7_register} */
1306 /* ,{0,"OPTSQP",&slv2_register} */
1307 ,{HAVE_CONOPT,"CONOPT",&slv8_register}
1308 ,{HAVE_CMSLV,"CMSLV",&slv9_register}
1309 ,{HAVE_LRSLV,"LRSLV",&slv9a_register}
1310 ,{0,NULL,NULL}
1311 };
1312
1313 int SlvRegisterStandardClients(void){
1314 int nclients = 0;
1315 int newclient=0;
1316 int error;
1317 int i;
1318 /* CONSOLE_DEBUG("REGISTERING STANDARD SOLVER ENGINES"); */
1319 for(i=0;slv_reg[i].name!=NULL;++i){
1320 if(slv_reg[i].is_active && slv_reg[i].regfunc){
1321 error = slv_register_client(slv_reg[i].regfunc,NULL,NULL,&newclient);
1322 if(error){
1323 ERROR_REPORTER_HERE(ASC_PROG_ERR
1324 ,"Unable to register solver '%s' (error %d)."
1325 ,slv_reg[i].name,error
1326 );
1327 }else{
1328 CONSOLE_DEBUG("Solver '%s' registered OK",slv_solver_name(newclient));
1329 nclients++;
1330 }
1331 }else{
1332 CONSOLE_DEBUG("Solver '%s' was not compiled.",slv_reg[i].name);
1333 }
1334 }
1335 return nclients;
1336 }
1337
1338 /*------------------------------------------------------------------------------
1339 OUTPUT ASSIGNMENT AND PARTITIONG IN LOGICAL RELATIONS
1340 */
1341
1342 #define MLIMDEBUG 0
1343 /*
1344 * Populates a matrix according to the sys solvers_dvars, solvers_logrels
1345 * lists and the filters given. mtx given must be created (not null)
1346 * and with order >= MAX( slv_get_num_solvers_logrels(sys),
1347 * slv_get_num_solvers_dvars(sys));
1348 * returns 0 if went ok.
1349 */
1350 int slv_std_make_log_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
1351 dis_filter_t *dvf,logrel_filter_t *lrf)
1352 {
1353 #if MLIMDEBUG
1354 FILE *fp;
1355 #endif
1356 int32 lr, lrlen,ord;
1357 struct logrel_relation **lrp;
1358
1359 if (sys==NULL || mtx == NULL || dvf == NULL || lrf == NULL) {
1360 FPRINTF(stderr,"make_log_incidence called with null\n");
1361 return 1;
1362 }
1363 lrp = slv_get_solvers_logrel_list(sys);
1364 assert(lrp!=NULL);
1365 lrlen = slv_get_num_solvers_logrels(sys);
1366 ord = MAX(lrlen,slv_get_num_solvers_dvars(sys));
1367 if (ord > mtx_order(mtx)) {
1368 FPRINTF(stderr,"make_incidence called with undersized matrix\n");
1369 return 2;
1370 }
1371 for (lr=0; lr < lrlen; lr++) {
1372 if (logrel_apply_filter(lrp[lr],lrf)) {
1373 logrelman_get_incidence(lrp[lr],dvf,mtx);
1374 }
1375 }
1376 #if MLIMDEBUG
1377 fp = fopen("/tmp/mim.plot","w+");
1378 if (fp !=NULL) {
1379 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1380 fclose(fp);
1381 }
1382 #endif
1383 return 0;
1384 }
1385
1386 #if USEDCODE
1387
1388 /* returns 0 if successful, 1 if insufficient memory. Does not change
1389 * the data in mtx. Orders the solvers_dvar list of the system to
1390 * match the permutation on the given mtx. It is assumed that
1391 * the org cols of mtx == dvar list position. Only the
1392 * dvars in range lo to hi of the dvar list are affected and
1393 * only these dvars should appear in the same org column range of
1394 * the mtx. This should not be called on blocks less than 3x3.
1395 */
1396 static int reindex_dvars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
1397 const mtx_matrix_t mtx)
1398 {
1399 struct dis_discrete **dvtmp, **dvp;
1400 int32 c,v,vlen;
1401
1402 if (lo >= hi +1) return 0; /* job too small */
1403 dvp = slv_get_solvers_dvar_list(sys);
1404 vlen = slv_get_num_solvers_dvars(sys);
1405 /* on dvtmp we DONT have the terminating null */
1406 dvtmp = (struct dis_discrete **)
1407 ascmalloc(vlen*sizeof(struct dis_discrete *));
1408 if (dvtmp == NULL) {
1409 return 1;
1410 }
1411 /* copy pointers to dvtmp in order desired and copy back rather than sort */
1412 for (c=lo;c<=hi;c++) {
1413 v = mtx_col_to_org(mtx,c);
1414 dvtmp[c] = dvp[v];
1415 }
1416 /* copying back and re-sindexing */
1417 for (c=lo;c<=hi;c++) {
1418 dvp[c] = dvtmp[c];
1419 dis_set_sindex(dvp[c],c);
1420 }
1421 ascfree(dvtmp);
1422 return 0;
1423 }
1424
1425 /* returns 0 if successful, 1 if insufficient memory. Does not change
1426 * the data in mtx. Orders the solvers_logrel list of the system to
1427 * match the permutation on the given mtx. It is assumed that
1428 * the org rows of mtx == logrel list position. Only the
1429 *logrels in range lo to hi of the logrel list are affected and
1430 * only these logrels should appear in the same org row range of
1431 * the input mtx. This should not be called on blocks less than 3x3.
1432 */
1433 static int reindex_logrels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
1434 const mtx_matrix_t mtx)
1435 {
1436 struct logrel_relation **lrtmp, **lrp;
1437 int32 c,v,rlen;
1438
1439 lrp = slv_get_solvers_logrel_list(sys);
1440 rlen = slv_get_num_solvers_logrels(sys);
1441 /* on lrtmp we DONT have the terminating null */
1442 lrtmp = (struct logrel_relation **)
1443 ascmalloc(rlen*sizeof(struct logrel_relation *));
1444 if (lrtmp == NULL) {
1445 return 1;
1446 }
1447 /* copy pointers to lrtmp in order desired and copy back rather than sort.
1448 * do this only in the row range of interest. */
1449 for (c=lo;c<=hi;c++) {
1450 v = mtx_row_to_org(mtx,c);
1451 #if RIDEBUG
1452 FPRINTF(stderr,"Old logrel sindex(org) %d becoming sindex(cur) %d\n",v,c);
1453 #endif
1454 lrtmp[c] = lrp[v];
1455 }
1456 /* copying back and re-sindexing */
1457 for (c=lo;c<=hi;c++) {
1458 lrp[c] = lrtmp[c];
1459 logrel_set_sindex(lrp[c],c);
1460 }
1461 ascfree(lrtmp);
1462 return 0;
1463 }
1464
1465 #endif /* USEDCODE */
1466
1467 /*
1468 * returns 0 if ok, OTHERWISE if madness detected.
1469 */
1470 int slv_log_block_partition(slv_system_t sys)
1471 {
1472 #if SLBPDEBUG
1473 FILE *fp;
1474 #endif
1475 struct logrel_relation **lrp;
1476 struct dis_discrete **dvp;
1477 mtx_region_t *newblocks;
1478 dof_t *d;
1479 int32 nrow,ncol;
1480 int32 ncolpfix, nrowpun, vlenmnv;
1481 mtx_matrix_t mtx;
1482 int32 order, rank;
1483 int32 c,len,dvlen,r,lrlen;
1484 dis_filter_t dvf;
1485 logrel_filter_t lrf;
1486
1487 lrp = slv_get_solvers_logrel_list(sys);
1488 dvp = slv_get_solvers_dvar_list(sys);
1489 lrlen = slv_get_num_solvers_logrels(sys);
1490 dvlen = slv_get_num_solvers_dvars(sys);
1491 if (lrlen ==0 || dvlen == 0) return 1;
1492 order = MAX(lrlen,dvlen);
1493
1494 lrf.matchbits = (LOGREL_ACTIVE);
1495 lrf.matchvalue = (LOGREL_ACTIVE);
1496 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1497 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1498
1499 /* nrow plus unincluded logrels */
1500 nrowpun = slv_count_solvers_logrels(sys,&lrf);
1501 /* ncol plus fixed discrete vars */
1502 ncolpfix = slv_count_solvers_dvars(sys,&dvf);
1503
1504 dvf.matchbits = (DIS_BVAR);
1505 dvf.matchvalue = 0;
1506
1507 /* dvlen minus non boolean vars */
1508 vlenmnv = dvlen - slv_count_solvers_dvars(sys,&dvf);
1509
1510 lrf.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
1511 lrf.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
1512 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_FIXED | DIS_ACTIVE);
1513 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1514
1515 nrow = slv_count_solvers_logrels(sys,&lrf);
1516 ncol = slv_count_solvers_dvars(sys,&dvf);
1517
1518 mtx = slv_get_sys_mtx(sys);
1519 mtx_set_order(mtx,order);
1520
1521 if (slv_std_make_log_incidence_mtx(sys,mtx,&dvf,&lrf)) {
1522 ERROR_REPORTER_HERE(ASC_PROG_ERR,"failure in creating incidence matrix.");
1523 mtx_destroy(mtx);
1524 return 1;
1525 }
1526
1527 mtx_output_assign(mtx,lrlen,dvlen);
1528 rank = mtx_symbolic_rank(mtx);
1529 if (rank == 0 ) return 1;
1530 if (rank < nrow) {
1531 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank < nrow. Dependent logical relations?");
1532 }
1533 if (rank < ncol) {
1534 if ( nrow != rank) {
1535 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank < ncol and nrow != rank. Excess of columns?");
1536 } else {
1537 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"rank<ncol but nrow==rank. Degrees of freedom?");
1538 }
1539 }
1540 if (ncol == nrow) {
1541 if (ncol == rank) {
1542 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"System of logical relations does not need Logic Inference.");
1543 }
1544 if (ncol != rank) {
1545 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"but ncol!=rank. Rank deficient?");
1546 }
1547 }
1548 mtx_partition(mtx);
1549 /* copy the block list. there is at least one if rank >=1 as above */
1550 len = mtx_number_of_blocks(mtx);
1551 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
1552 if (newblocks == NULL) {
1553 mtx_destroy(mtx);
1554 return 2;
1555 }
1556 for (c = 0 ; c < len; c++) {
1557 mtx_block(mtx,c,&(newblocks[c]));
1558 }
1559 slv_set_solvers_log_blocks(sys,len,newblocks); /* give it to the system */
1560 d = slv_get_log_dofdata(sys);
1561 d->structural_rank = rank;
1562 d->n_rows = nrow;
1563 d->n_cols = ncol;
1564 d->n_fixed = ncolpfix - ncol;
1565 d->n_unincluded = nrowpun - nrow;
1566 d->reorder.partition = 1; /* yes */
1567 d->reorder.basis_selection = 0; /* none yet */
1568 d->reorder.block_reordering = 0; /* none */
1569
1570 #if SLBPDEBUG
1571 fp = fopen("/tmp/sbp1.plot","w+");
1572 if (fp !=NULL) {
1573 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1574 fclose(fp);
1575 }
1576 #endif
1577
1578 len = vlenmnv - 1; /* row of last possible boolean variable */
1579 for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
1580 r = mtx_col_to_org(mtx,c);
1581 if ( dis_fixed(dvp[r]) && dis_active(dvp[r])) {
1582 c++;
1583 } else {
1584 mtx_swap_cols(mtx,c,len);
1585 len--;
1586 }
1587 }
1588
1589 /* Now, leftover columns are where we want them.
1590 * That is, unassigned incidence is next to assigned and fixed on right with
1591 * no columns which do not correspond to dis vars in the range 0..dvlen-1.
1592 */
1593 /* there aren't any nonincident dvars in the solvers dvarlist, so
1594 * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
1595
1596 /* rows 0->lrlen-1 should be actual logrels, though, and not phantoms.
1597 * At this point we should have
1598 * rows - assigned
1599 * - unassigned w/incidence
1600 * - unincluded and included w/out free incidence and inactive mixed.
1601 */
1602
1603 /* sort unincluded logrelations to bottom of good region */
1604 /* this is needed because the idiot user may have fixed everything */
1605 /* in a logrelation, but we want those rows not to be confused with */
1606 /* unincluded rows.*/
1607 len = lrlen - 1; /* row of last possible relation */
1608 /* sort the inactive out of the unassigned/unincluded */
1609 for (c=rank; len > c ; ) {
1610 r = mtx_row_to_org(mtx,c);
1611 if ( logrel_active(lrp[r])) {
1612 c++;
1613 } else {
1614 mtx_swap_rows(mtx,c,len);
1615 len--;
1616 }
1617 }
1618
1619 /* sort the unincluded out of the unassigned */
1620 len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
1621 for (c=rank; len > c ; ) {
1622 r = mtx_row_to_org(mtx,c);
1623 if ( logrel_included(lrp[r]) ) {
1624 c++;
1625 } else {
1626 mtx_swap_rows(mtx,c,len);
1627 len--;
1628 }
1629 }
1630
1631 #if SLBPDEBUG
1632 fp = fopen("/tmp/sbp2.plot","w+");
1633 if(fp !=NULL){
1634 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1635 fclose(fp);
1636 }
1637 #endif
1638
1639 /* These functions are called for slv_block_partitioning, but I do not
1640 * know if I need them here:
1641 *
1642 * now, at last we have cols in the order we want the lists to
1643 * be handed to the solver. So, we need to reset the dvar sindex values
1644 * and reorder the lists.
1645 */
1646
1647 #if USEDCODE
1648 if (reindex_dvars_from_mtx(sys,0,dvlen-1,mtx)) {
1649 mtx_destroy(mtx);
1650 return 2;
1651 }
1652 #endif /* USEDCODE */
1653
1654 /*
1655 * now, at last we have rows jacobian in the order we want the lists to
1656 * be handed to the solvers. So, we need to reset the rel sindex values.
1657 */
1658
1659 #if USEDCODE
1660 if (reindex_logrels_from_mtx(sys,0,lrlen-1,mtx)) {
1661 mtx_destroy(mtx);
1662 return 2;
1663 }
1664 #endif /* USEDCODE */
1665
1666 return 0;
1667 }
1668

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