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

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