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

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