/[ascend]/trunk/ascend4/solver/slv_stdcalls.c
ViewVC logotype

Contents of /trunk/ascend4/solver/slv_stdcalls.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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