/[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 731 - (show annotations) (download) (as text)
Tue Jul 4 07:42:06 2006 UTC (14 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 47489 byte(s)
Removed some debug messages.
Fixed up return values for Integrators functions to comply with integrator.c API.
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 CONSOLE_DEBUG("TEARDROP reordering...");
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_insure_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 FPRINTF(mif,"Bounds for variable ");
1185 var_write_name(sys,var,mif);
1186 FPRINTF(mif," are inconsistent [%g,%g].\n",low,high);
1187 FPRINTF(mif,"Bounds will be swapped.\n");
1188 }
1189 var_set_upper_bound(var, low);
1190 var_set_lower_bound(var, high);
1191 low = var_lower_bound(var);
1192 high = var_upper_bound(var);
1193 nchange++;
1194 }
1195
1196 if( low > val ) {
1197 if (mif!=NULL) {
1198 FPRINTF(mif,"Variable ");
1199 var_write_name(sys,var,mif);
1200 FPRINTF(mif," was initialized below its lower bound.\n");
1201 FPRINTF(mif,"It will be moved to its lower bound.\n");
1202 }
1203 var_set_value(var, low);
1204 } else {
1205 if( val > high ) {
1206 if (mif!=NULL) {
1207 FPRINTF(mif,"Variable ");
1208 var_write_name(sys,var,mif);
1209 FPRINTF(mif," was initialized above its upper bound.\n");
1210 FPRINTF(mif,"It will be moved to its upper bound.\n");
1211 }
1212 var_set_value(var, high);
1213 nchange++;
1214 }
1215 }
1216 }
1217 return nchange;
1218 }
1219
1220
1221 void slv_check_bounds(const slv_system_t sys,int32 lo,int32 hi,
1222 FILE *mif,const char *label)
1223 {
1224 real64 val,low,high;
1225 int32 c,len;
1226 struct var_variable *var, **vp;
1227 static char defaultlabel[] = "";
1228
1229 if (label==NULL) label = defaultlabel;
1230 if (sys==NULL || mif == NULL) return;
1231 vp = slv_get_solvers_var_list(sys);
1232 if (vp==NULL) return;
1233 len = slv_get_num_solvers_vars(sys);
1234 if (lo > len || hi > len) {
1235 FPRINTF(stderr,"slv_check_bounds miscalled\n");
1236 return;
1237 }
1238 for (c= lo; c <= hi; c++) {
1239 var = vp[c];
1240 low = var_lower_bound(var);
1241 high = var_upper_bound(var);
1242 val = var_value(var);
1243 if( low > high ) {
1244 FPRINTF(mif,"Bounds for %s variable ",label);
1245 var_write_name(sys,var,mif);
1246 FPRINTF(mif,"\nare inconsistent [%g,%g].\n",low,high);
1247 }
1248
1249 if( low > val ) {
1250 FPRINTF(mif,"%s variable ",label);
1251 var_write_name(sys,var,mif);
1252 FPRINTF(mif,"\nwas initialized below its lower bound.\n");
1253 } else {
1254 if( val > high ) {
1255 FPRINTF(mif,"%s variable ",label);
1256 var_write_name(sys,var,mif);
1257 FPRINTF(mif," was initialized above its upper bound.\n");
1258 }
1259 }
1260 }
1261 return;
1262 }
1263
1264 /*------------------------------------------------------------------------------
1265 SOLVER REGISTRATION
1266 */
1267
1268 /* rewrote this stuff to get rid of all the #ifdefs -- JP */
1269
1270 struct StaticSolverRegistration{
1271 int is_active;
1272 const char *name;
1273 SlvRegistration *regfunc;
1274 };
1275
1276 /*
1277 The names here are only used to provide information in the case where
1278 solver registration fails. The definitive solver names are in the slv*.c
1279 files.
1280 */
1281 static const struct StaticSolverRegistration slv_reg[]={
1282 /* {0,"SLV",&slv0_register} */
1283 /* ,{0,"MINOS",&slv1_register} */
1284 {HAVE_QRSLV,"QRSLV",&slv3_register}
1285 /* ,{0,"CSLV",&slv4_register} */
1286 /* ,{0,"LSSLV",&slv5_register} */
1287 /* ,{0,"MPS",&slv6_register} */
1288 /* ,{0,"NGSLV",&slv7_register} */
1289 /* ,{0,"OPTSQP",&slv2_register} */
1290 ,{HAVE_CONOPT,"CONOPT",&slv8_register}
1291 ,{HAVE_CMSLV,"CMSLV",&slv9_register}
1292 /* ,{0,"LRSLV",&slv9a_register} */
1293 ,{0,NULL,NULL}
1294 };
1295
1296 int SlvRegisterStandardClients(void){
1297 int nclients = 0;
1298 int newclient=0;
1299 int error;
1300 int i;
1301 /* CONSOLE_DEBUG("REGISTERING STANDARD SOLVER ENGINES"); */
1302 for(i=0;slv_reg[i].name!=NULL;++i){
1303 if(slv_reg[i].is_active && slv_reg[i].regfunc){
1304 error = slv_register_client(slv_reg[i].regfunc,NULL,NULL,&newclient);
1305 if(error){
1306 ERROR_REPORTER_HERE(ASC_PROG_ERR
1307 ,"Unable to register solver '%s' (error %d)."
1308 ,slv_reg[i].name,error
1309 );
1310 }else{
1311 CONSOLE_DEBUG("Solver '%s' registered OK",slv_solver_name(newclient));
1312 nclients++;
1313 }
1314 }else{
1315 CONSOLE_DEBUG("Solver '%s' was not compiled.",slv_reg[i].name);
1316 }
1317 }
1318 return nclients;
1319 }
1320
1321 /*------------------------------------------------------------------------------
1322 OUTPUT ASSIGNMENT AND PARTITIONG IN LOGICAL RELATIONS
1323 */
1324
1325 #define MLIMDEBUG 0
1326 /*
1327 * Populates a matrix according to the sys solvers_dvars, solvers_logrels
1328 * lists and the filters given. mtx given must be created (not null)
1329 * and with order >= MAX( slv_get_num_solvers_logrels(sys),
1330 * slv_get_num_solvers_dvars(sys));
1331 * returns 0 if went ok.
1332 */
1333 int slv_std_make_log_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
1334 dis_filter_t *dvf,logrel_filter_t *lrf)
1335 {
1336 #if MLIMDEBUG
1337 FILE *fp;
1338 #endif
1339 int32 lr, lrlen,ord;
1340 struct logrel_relation **lrp;
1341
1342 if (sys==NULL || mtx == NULL || dvf == NULL || lrf == NULL) {
1343 FPRINTF(stderr,"make_log_incidence called with null\n");
1344 return 1;
1345 }
1346 lrp = slv_get_solvers_logrel_list(sys);
1347 assert(lrp!=NULL);
1348 lrlen = slv_get_num_solvers_logrels(sys);
1349 ord = MAX(lrlen,slv_get_num_solvers_dvars(sys));
1350 if (ord > mtx_order(mtx)) {
1351 FPRINTF(stderr,"make_incidence called with undersized matrix\n");
1352 return 2;
1353 }
1354 for (lr=0; lr < lrlen; lr++) {
1355 if (logrel_apply_filter(lrp[lr],lrf)) {
1356 logrelman_get_incidence(lrp[lr],dvf,mtx);
1357 }
1358 }
1359 #if MLIMDEBUG
1360 fp = fopen("/tmp/mim.plot","w+");
1361 if (fp !=NULL) {
1362 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1363 fclose(fp);
1364 }
1365 #endif
1366 return 0;
1367 }
1368
1369 #if USEDCODE
1370
1371 /* returns 0 if successful, 1 if insufficient memory. Does not change
1372 * the data in mtx. Orders the solvers_dvar list of the system to
1373 * match the permutation on the given mtx. It is assumed that
1374 * the org cols of mtx == dvar list position. Only the
1375 * dvars in range lo to hi of the dvar list are affected and
1376 * only these dvars should appear in the same org column range of
1377 * the mtx. This should not be called on blocks less than 3x3.
1378 */
1379 static int reindex_dvars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
1380 const mtx_matrix_t mtx)
1381 {
1382 struct dis_discrete **dvtmp, **dvp;
1383 int32 c,v,vlen;
1384
1385 if (lo >= hi +1) return 0; /* job too small */
1386 dvp = slv_get_solvers_dvar_list(sys);
1387 vlen = slv_get_num_solvers_dvars(sys);
1388 /* on dvtmp we DONT have the terminating null */
1389 dvtmp = (struct dis_discrete **)
1390 ascmalloc(vlen*sizeof(struct dis_discrete *));
1391 if (dvtmp == NULL) {
1392 return 1;
1393 }
1394 /* copy pointers to dvtmp in order desired and copy back rather than sort */
1395 for (c=lo;c<=hi;c++) {
1396 v = mtx_col_to_org(mtx,c);
1397 dvtmp[c] = dvp[v];
1398 }
1399 /* copying back and re-sindexing */
1400 for (c=lo;c<=hi;c++) {
1401 dvp[c] = dvtmp[c];
1402 dis_set_sindex(dvp[c],c);
1403 }
1404 ascfree(dvtmp);
1405 return 0;
1406 }
1407
1408 /* returns 0 if successful, 1 if insufficient memory. Does not change
1409 * the data in mtx. Orders the solvers_logrel list of the system to
1410 * match the permutation on the given mtx. It is assumed that
1411 * the org rows of mtx == logrel list position. Only the
1412 *logrels in range lo to hi of the logrel list are affected and
1413 * only these logrels should appear in the same org row range of
1414 * the input mtx. This should not be called on blocks less than 3x3.
1415 */
1416 static int reindex_logrels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
1417 const mtx_matrix_t mtx)
1418 {
1419 struct logrel_relation **lrtmp, **lrp;
1420 int32 c,v,rlen;
1421
1422 lrp = slv_get_solvers_logrel_list(sys);
1423 rlen = slv_get_num_solvers_logrels(sys);
1424 /* on lrtmp we DONT have the terminating null */
1425 lrtmp = (struct logrel_relation **)
1426 ascmalloc(rlen*sizeof(struct logrel_relation *));
1427 if (lrtmp == NULL) {
1428 return 1;
1429 }
1430 /* copy pointers to lrtmp in order desired and copy back rather than sort.
1431 * do this only in the row range of interest. */
1432 for (c=lo;c<=hi;c++) {
1433 v = mtx_row_to_org(mtx,c);
1434 #if RIDEBUG
1435 FPRINTF(stderr,"Old logrel sindex(org) %d becoming sindex(cur) %d\n",v,c);
1436 #endif
1437 lrtmp[c] = lrp[v];
1438 }
1439 /* copying back and re-sindexing */
1440 for (c=lo;c<=hi;c++) {
1441 lrp[c] = lrtmp[c];
1442 logrel_set_sindex(lrp[c],c);
1443 }
1444 ascfree(lrtmp);
1445 return 0;
1446 }
1447
1448 #endif /* USEDCODE */
1449
1450 /*
1451 * returns 0 if ok, OTHERWISE if madness detected.
1452 */
1453 int slv_log_block_partition(slv_system_t sys)
1454 {
1455 #if SLBPDEBUG
1456 FILE *fp;
1457 #endif
1458 struct logrel_relation **lrp;
1459 struct dis_discrete **dvp;
1460 mtx_region_t *newblocks;
1461 dof_t *d;
1462 int32 nrow,ncol;
1463 int32 ncolpfix, nrowpun, vlenmnv;
1464 mtx_matrix_t mtx;
1465 int32 order, rank;
1466 int32 c,len,dvlen,r,lrlen;
1467 dis_filter_t dvf;
1468 logrel_filter_t lrf;
1469
1470 lrp = slv_get_solvers_logrel_list(sys);
1471 dvp = slv_get_solvers_dvar_list(sys);
1472 lrlen = slv_get_num_solvers_logrels(sys);
1473 dvlen = slv_get_num_solvers_dvars(sys);
1474 if (lrlen ==0 || dvlen == 0) return 1;
1475 order = MAX(lrlen,dvlen);
1476
1477 lrf.matchbits = (LOGREL_ACTIVE);
1478 lrf.matchvalue = (LOGREL_ACTIVE);
1479 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1480 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1481
1482 /* nrow plus unincluded logrels */
1483 nrowpun = slv_count_solvers_logrels(sys,&lrf);
1484 /* ncol plus fixed discrete vars */
1485 ncolpfix = slv_count_solvers_dvars(sys,&dvf);
1486
1487 dvf.matchbits = (DIS_BVAR);
1488 dvf.matchvalue = 0;
1489
1490 /* dvlen minus non boolean vars */
1491 vlenmnv = dvlen - slv_count_solvers_dvars(sys,&dvf);
1492
1493 lrf.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
1494 lrf.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
1495 dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_FIXED | DIS_ACTIVE);
1496 dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1497
1498 nrow = slv_count_solvers_logrels(sys,&lrf);
1499 ncol = slv_count_solvers_dvars(sys,&dvf);
1500
1501 mtx = slv_get_sys_mtx(sys);
1502 mtx_set_order(mtx,order);
1503
1504 if (slv_std_make_log_incidence_mtx(sys,mtx,&dvf,&lrf)) {
1505 FPRINTF(stderr,
1506 "slv_log_block_partition: failure in creating incidence matrix.\n");
1507 mtx_destroy(mtx);
1508 return 1;
1509 }
1510
1511 mtx_output_assign(mtx,lrlen,dvlen);
1512 rank = mtx_symbolic_rank(mtx);
1513 if (rank == 0 ) return 1;
1514 if (rank < nrow) {
1515 FPRINTF(stderr,"rank<nrow in slv_log_block_partition\n");
1516 FPRINTF(stderr,"dependent logical relations ? \n");
1517 }
1518 if (rank < ncol) {
1519 if ( nrow != rank) {
1520 FPRINTF(stderr,"rank<ncol and nrow!=rank in slv_log_block_partition.\n");
1521 FPRINTF(stderr,"Excess of columns ? \n");
1522 } else {
1523 FPRINTF(stderr,"rank<ncol but nrow==rank in slv_log_block_partition.\n");
1524 FPRINTF(stderr,"Degrees of freedom ? \n");
1525 }
1526 }
1527 if (ncol == nrow) {
1528 if (ncol == rank) {
1529 FPRINTF(stderr,"\n");
1530 FPRINTF(stderr,
1531 "System of logical relations does not need Logic Inference \n");
1532 }
1533 if (ncol != rank) {
1534 FPRINTF(stderr,"but ncol!=rank.\n");
1535 FPRINTF(stderr,"Rank deficient ? \n");
1536 } else {
1537 FPRINTF(stderr,"\n");
1538 }
1539 }
1540 mtx_partition(mtx);
1541 /* copy the block list. there is at least one if rank >=1 as above */
1542 len = mtx_number_of_blocks(mtx);
1543 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
1544 if (newblocks == NULL) {
1545 mtx_destroy(mtx);
1546 return 2;
1547 }
1548 for (c = 0 ; c < len; c++) {
1549 mtx_block(mtx,c,&(newblocks[c]));
1550 }
1551 slv_set_solvers_log_blocks(sys,len,newblocks); /* give it to the system */
1552 d = slv_get_log_dofdata(sys);
1553 d->structural_rank = rank;
1554 d->n_rows = nrow;
1555 d->n_cols = ncol;
1556 d->n_fixed = ncolpfix - ncol;
1557 d->n_unincluded = nrowpun - nrow;
1558 d->reorder.partition = 1; /* yes */
1559 d->reorder.basis_selection = 0; /* none yet */
1560 d->reorder.block_reordering = 0; /* none */
1561
1562 #if SLBPDEBUG
1563 fp = fopen("/tmp/sbp1.plot","w+");
1564 if (fp !=NULL) {
1565 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1566 fclose(fp);
1567 }
1568 #endif
1569
1570 len = vlenmnv - 1; /* row of last possible boolean variable */
1571 for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
1572 r = mtx_col_to_org(mtx,c);
1573 if ( dis_fixed(dvp[r]) && dis_active(dvp[r])) {
1574 c++;
1575 } else {
1576 mtx_swap_cols(mtx,c,len);
1577 len--;
1578 }
1579 }
1580
1581 /* Now, leftover columns are where we want them.
1582 * That is, unassigned incidence is next to assigned and fixed on right with
1583 * no columns which do not correspond to dis vars in the range 0..dvlen-1.
1584 */
1585 /* there aren't any nonincident dvars in the solvers dvarlist, so
1586 * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
1587
1588 /* rows 0->lrlen-1 should be actual logrels, though, and not phantoms.
1589 * At this point we should have
1590 * rows - assigned
1591 * - unassigned w/incidence
1592 * - unincluded and included w/out free incidence and inactive mixed.
1593 */
1594
1595 /* sort unincluded logrelations to bottom of good region */
1596 /* this is needed because the idiot user may have fixed everything */
1597 /* in a logrelation, but we want those rows not to be confused with */
1598 /* unincluded rows.*/
1599 len = lrlen - 1; /* row of last possible relation */
1600 /* sort the inactive out of the unassigned/unincluded */
1601 for (c=rank; len > c ; ) {
1602 r = mtx_row_to_org(mtx,c);
1603 if ( logrel_active(lrp[r])) {
1604 c++;
1605 } else {
1606 mtx_swap_rows(mtx,c,len);
1607 len--;
1608 }
1609 }
1610
1611 /* sort the unincluded out of the unassigned */
1612 len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
1613 for (c=rank; len > c ; ) {
1614 r = mtx_row_to_org(mtx,c);
1615 if ( logrel_included(lrp[r]) ) {
1616 c++;
1617 } else {
1618 mtx_swap_rows(mtx,c,len);
1619 len--;
1620 }
1621 }
1622
1623 #if SLBPDEBUG
1624 fp = fopen("/tmp/sbp2.plot","w+");
1625 if(fp !=NULL){
1626 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1627 fclose(fp);
1628 }
1629 #endif
1630
1631 /* These functions are called for slv_block_partitioning, but I do not
1632 * know if I need them here:
1633 *
1634 * now, at last we have cols in the order we want the lists to
1635 * be handed to the solver. So, we need to reset the dvar sindex values
1636 * and reorder the lists.
1637 */
1638
1639 #if USEDCODE
1640 if (reindex_dvars_from_mtx(sys,0,dvlen-1,mtx)) {
1641 mtx_destroy(mtx);
1642 return 2;
1643 }
1644 #endif /* USEDCODE */
1645
1646 /*
1647 * now, at last we have rows jacobian in the order we want the lists to
1648 * be handed to the solvers. So, we need to reset the rel sindex values.
1649 */
1650
1651 #if USEDCODE
1652 if (reindex_logrels_from_mtx(sys,0,lrlen-1,mtx)) {
1653 mtx_destroy(mtx);
1654 return 2;
1655 }
1656 #endif /* USEDCODE */
1657
1658 return 0;
1659 }
1660

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