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