/* * Solver Standard Clients * by Benjamin Andrew Allan * 5/19/96 * Version: $Revision: 1.28 $ * Version control file: $RCSfile: slv_stdcalls.c,v $ * Date last modified: $Date: 1998/06/16 16:53:04 $ * Last modified by: $Author: mthomas $ * * Copyright(C) 1996 Benjamin Andrew Allan * Copyright(C) 1998 Carnegie Mellon University * * This file is part of the ASCEND IV math programming system. * Here is where we register our normal solvers and keep things * like reordering clients and so forth. This file should probably * be broken up into one file-one client layout for the unregistered * clients so that we can load them selectively. * * The Ascend Math Programming System is free software; you can * redistribute it and/or modify it under the terms of the GNU * General Public License as published by the Free Software * Foundation; either version 2 of the License, or (at your option) * any later version. * * The Ascend Math Programming System is distributed in hope that it * will be useful, but WITHOUT ANY WARRANTY; without even the implied * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. * See the GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with the program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check * the file named COPYING. */ #include "utilities/ascConfig.h" #include "compiler/compiler.h" #include "utilities/ascMalloc.h" #include "general/list.h" #include "solver/mtx.h" #include "solver/slv_types.h" #include "solver/var.h" #include "solver/rel.h" #include "solver/discrete.h" #include "solver/conditional.h" #include "solver/logrel.h" #include "solver/relman.h" #include "solver/logrelman.h" #include "solver/bnd.h" #include "solver/slv_common.h" #include "solver/linsol.h" #include "solver/linsolqr.h" #include "solver/slv_client.h" /* header of registered clients */ #include "solver/slv0.h" #include "solver/slv1.h" #include "solver/slv2.h" #include "solver/slv3.h" #include "solver/slv6.h" #include "solver/slv7.h" #include "solver/slv8.h" #include "solver/slv9.h" #include "solver/slv9a.h" #include "solver/model_reorder.h" #include "solver/slv_stdcalls.h" #define KILL 0 /* deleteme code. compile old code if kill = 1 */ #define NEEDSTOBEDONE 0 /* indicates code/comments that are not yet ready */ #define USECODE 0 /* Code in good shape but not used currently */ /* * Here we play some hefty jacobian reordering games. * What we want to happen eventually is as follows: * * Get the free & incident pattern for include relations. * Output assign jacobian. * BLT permute the jacobian. if underspecified, fake rows * to make things appear square. * For all leading square blocks apply reordering (kirk, etc) * If trailing rectangular block, apply clever kirkbased scheme not known. * * At present, we aren't quite so clever. We: * * Get the free & incident pattern for include relations. * Output assign jacobian. * BLT permute the square output assigned region of the jacobian. * For all leading square blocks apply reordering (kirk, etc) * * Collect the block list as part of the master data structure. * Set sindices for rels, vars as current rows/cols in matrix so * that the jacobian is 'naturally' preordered for all solvers. * * Solvers are still free to reorder their own matrices any way they like; * it's probably a dumb idea, though. * * returns 0 if ok, 1 if out of memory. */ #define MIMDEBUG 0 /* * Populates a matrix according to the sys solvers_vars, solvers_rels * lists and the filters given. The filter should have at least * SVAR = 1 bit on. mtx given must be created (not null) and with * order >= MAX( slv_get_num_solvers_rels(sys),slv_get_num_solvers_vars(sys)); * returns 0 if went ok. */ int slv_std_make_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx, var_filter_t *vf,rel_filter_t *rf) { #if MIMDEBUG FILE *fp; #endif int32 r, rlen,ord; struct rel_relation **rp; if (sys==NULL || mtx == NULL || vf == NULL || rf == NULL) { FPRINTF(stderr,"make_incidence called with null\n"); return 1; } rp = slv_get_solvers_rel_list(sys); assert(rp!=NULL); rlen = slv_get_num_solvers_rels(sys); ord = MAX(rlen,slv_get_num_solvers_vars(sys)); if (ord > mtx_order(mtx)) { FPRINTF(stderr,"make_incidence called with undersized matrix\n"); return 2; } for (r=0; r < rlen; r++) { if (rel_apply_filter(rp[r],rf)) { relman_get_incidence(rp[r],vf,mtx); } } #if MIMDEBUG fp = fopen("/tmp/mim.plot","w+"); if (fp !=NULL) { mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX); fclose(fp); } #endif return 0; } #define RIDEBUG 0 /* returns 0 if successful, 1 if insufficient memory. Does not change * the data in mtx. Orders the solvers_var list of the system to * match the permutation on the given mtx. It is assumed that * the org cols of mtx == var list position. Only the * vars in range lo to hi of the var list are affected and * only these vars should appear in the same org column range of * the mtx. This should not be called on blocks less than 3x3. */ static int reindex_vars_from_mtx(slv_system_t sys, int32 lo, int32 hi, const mtx_matrix_t mtx) { struct var_variable **vtmp, **vp; int32 c,v,vlen; if (lo >= hi +1) return 0; /* job too small */ vp = slv_get_solvers_var_list(sys); vlen = slv_get_num_solvers_vars(sys); /* on vtmp we DONT have the terminating null */ vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *)); if (vtmp == NULL) { return 1; } /* copy pointers to vtmp in order desired and copy back rather than sort */ for (c=lo;c<=hi;c++) { v = mtx_col_to_org(mtx,c); vtmp[c] = vp[v]; } /* copying back and re-sindexing */ for (c=lo;c<=hi;c++) { vp[c] = vtmp[c]; var_set_sindex(vp[c],c); } ascfree(vtmp); return 0; } /* returns 0 if successful, 1 if insufficient memory. Does not change * the data in mtx. Orders the solvers_rel list of the system to * match the permutation on the given mtx. It is assumed that * the org rows of mtx == rel list position. Only the * rels in range lo to hi of the rel list are affected and * only these rels should appear in the same org row range of * the input mtx. This should not be called on blocks less than 3x3. */ static int reindex_rels_from_mtx(slv_system_t sys, int32 lo, int32 hi, const mtx_matrix_t mtx) { struct rel_relation **rtmp, **rp; int32 c,v,rlen; rp = slv_get_solvers_rel_list(sys); rlen = slv_get_num_solvers_rels(sys); /* on rtmp we DONT have the terminating null */ rtmp = (struct rel_relation **)ascmalloc(rlen*sizeof(struct rel_relation *)); if (rtmp == NULL) { return 1; } /* copy pointers to rtmp in order desired and copy back rather than sort. * do this only in the row range of interest. */ for (c=lo;c<=hi;c++) { v = mtx_row_to_org(mtx,c); #if RIDEBUG FPRINTF(stderr,"Old rel sindex (org) %d becoming sindex (cur) %d\n",v,c); #endif rtmp[c] = rp[v]; } /* copying back and re-sindexing */ for (c=lo;c<=hi;c++) { rp[c] = rtmp[c]; rel_set_sindex(rp[c],c); } ascfree(rtmp); return 0; } /* returns 0 if ok, OTHERWISE if madness detected. * doesn't grok inequalities. */ #define SBPDEBUG 0 int slv_block_partition_real(slv_system_t sys,int uppertriangular) { #if SBPDEBUG FILE *fp; #endif struct rel_relation **rp; struct var_variable **vp; mtx_region_t *newblocks; dof_t *d; int32 nrow,ncol; int32 ncolpfix, nrowpun, vlenmnv; mtx_matrix_t mtx; int32 order, rank; int32 c,len,vlen,r,rlen; var_filter_t vf; rel_filter_t rf; rp = slv_get_solvers_rel_list(sys); vp = slv_get_solvers_var_list(sys); rlen = slv_get_num_solvers_rels(sys); vlen = slv_get_num_solvers_vars(sys); if (rlen ==0 || vlen == 0) return 1; order = MAX(rlen,vlen); rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE); rf.matchvalue = (REL_ACTIVE); vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE); vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE); /* nrow plus unincluded rels and ineqs */ nrowpun = slv_count_solvers_rels(sys,&rf); ncolpfix = slv_count_solvers_vars(sys,&vf); /* ncol plus fixed vars */ vf.matchbits = (VAR_SVAR); vf.matchvalue = 0; vlenmnv = vlen - slv_count_solvers_vars(sys,&vf); /* vlen minus non vars */ rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE); rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE); vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE); vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE); nrow = slv_count_solvers_rels(sys,&rf); ncol = slv_count_solvers_vars(sys,&vf); mtx = mtx_create(); mtx_set_order(mtx,order); if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) { FPRINTF(stderr, "slv_block_partition: failure in creating incidence matrix.\n"); mtx_destroy(mtx); return 1; } mtx_output_assign(mtx,rlen,vlen); rank = mtx_symbolic_rank(mtx); if (rank == 0 ) return 1; /* nothing to do, eh? */ /* lot of whining about dof */ if (rank < nrow) { FPRINTF(stderr,"System is row rank deficient (%d dependent equations)\n", nrow - rank); } if (rank < ncol) { if ( nrow != rank) { FPRINTF(stderr,"System is row rank deficient with %d excess columns.\n", ncol - rank); } else { FPRINTF(stderr,"System has %d degrees of freedom.\n", ncol - rank); } } if (ncol == nrow) { FPRINTF(stderr,"System is (%d) square",ncol); if (ncol != rank) { FPRINTF(stderr,"but rank deficient.\n"); } else { FPRINTF(stderr,".\n"); } } if (uppertriangular) { mtx_ut_partition(mtx); } else { mtx_partition(mtx); } /* copy the block list. there is at least one if rank >=1 as above */ len = mtx_number_of_blocks(mtx); newblocks = (mtx_region_t *)ascmalloc(len*sizeof(mtx_region_t)); if (newblocks == NULL) { mtx_destroy(mtx); return 2; } for (c = 0 ; c < len; c++) { mtx_block(mtx,c,&(newblocks[c])); } slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */ d = slv_get_dofdata(sys); d->structural_rank = rank; d->n_rows = nrow; d->n_cols = ncol; /* the next two lines assume inequalities are unincluded. */ #define MINE 1 #if MINE d->n_fixed = ncolpfix - ncol; d->n_unincluded = nrowpun - nrow; #else d->n_fixed = vlen - ncol; d->n_unincluded = rlen - nrow; #endif d->reorder.partition = 1; /* yes */ d->reorder.basis_selection = 0; /* none yet */ d->reorder.block_reordering = 0; /* none */ #if SBPDEBUG fp = fopen("/tmp/sbp1.plot","w+"); if (fp !=NULL) { mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX); fclose(fp); } #endif #if MINE len = vlenmnv - 1; /* row of last possible solver variable */ for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */ r = mtx_col_to_org(mtx,c); if ( var_fixed(vp[r]) && var_active(vp[r])) { c++; } else { mtx_swap_cols(mtx,c,len); len--; } } #endif /* Now, leftover columns are where we want them. * That is, unassigned incidence is next to assigned and fixed on right with * no columns which do not correspond to variables in the range 0..vlen-1. */ /* there aren't any nonincident vars in the solvers varlist, so * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */ /* rows 0->rlen-1 should be actual rels, though, and not phantoms. * At this point we should have * rows - assigned * - unassigned w/incidence * - unincluded and included w/out free incidence and inactive mixed. */ /* sort unincluded relations to bottom of good region */ /* this is needed because the idiot user may have fixed everything */ /* in a relation, but we want those rows not to be confused with */ /* unincluded rows. we're sort of dropping inequalities on the floor.*/ len = rlen - 1; /* row of last possible relation */ /* sort the inactive out of the unassigned/unincluded */ for (c=rank; len > c ; ) { r = mtx_row_to_org(mtx,c); if ( rel_active(rp[r])) { c++; } else { mtx_swap_rows(mtx,c,len); len--; } } #if MINE /* sort the unincluded out of the unassigned */ len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */ for (c=rank; len > c ; ) { r = mtx_row_to_org(mtx,c); if ( rel_included(rp[r]) ) { c++; } else { mtx_swap_rows(mtx,c,len); len--; } } #endif #if SBPDEBUG fp = fopen("/tmp/sbp2.plot","w+"); if (fp !=NULL) { mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX); fclose(fp); } #endif /* now, at last we have cols jacobian in the order we want the lists to * be handed to the solvers. So, we need to reset the var sindex values * and reorder the lists. */ if (reindex_vars_from_mtx(sys,0,vlen-1,mtx)) { mtx_destroy(mtx); return 2; } /* now, at last we have rows jacobian in the order we want the lists to * be handed to the solvers. So, we need to reset the rel sindex values. */ if (reindex_rels_from_mtx(sys,0,rlen-1,mtx)) { mtx_destroy(mtx); return 2; } mtx_destroy(mtx); return 0; } #if 0 /* code not currently used */ #ifdef STATIC_HARWELL extern void mc21b(); extern void mc13emod(); /* returns 0 if ok, OTHERWISE if madness detected. * doesn't grok inequalities. * CURRENTLY DOESN'T DO WELL WHEN NCOL c; ) { rel = rtmp[c]; if (rel_active(rel)) { c++; } else { rtmp[c] = rtmp[len]; rtmp[len] = rel; len--; } } /* Sort variable list */ var_perm_count = 0; var_tmp_end = vlen - 1; for (c = 0; c < vlen; c++) { var = vp[c]; if (var_apply_filter(var,&vf)) { vtmp[var_perm_count] = var; var_perm_count++; } else { vtmp[var_tmp_end] = var; var_tmp_end--; } } /* sort the inactive out of the unincluded and move to end */ len = vlen -1; for(c = var_tmp_end + 1; len > c; ) { var = vtmp[c]; if (var_active(var) && var_active(var)) { c++; } else { vtmp[c] = vtmp[len]; vtmp[len] = var; len--; } } /* reset solver indicies to make life easier */ for (c = 0; c < vlen; c++) { vp[c] = vtmp[c]; var_set_sindex(vp[c],c); } size = MAX(nrow,ncol); /* Create vectors for fortran calls */ icn = (int32 *)ascmalloc(licn*sizeof(int32)); ip = (int32 *)ascmalloc(size*sizeof(int32)); lenr = (int32 *)ascmalloc(size*sizeof(int32)); iperm = (int32 *)ascmalloc(size*sizeof(int32)); iw1 = (int32 *)ascmalloc(size*sizeof(int32)); iw2 = (int32 *)ascmalloc(size*sizeof(int32)); iw3 = (int32 *)ascmalloc(size*sizeof(int32)); arp = (int32 *)ascmalloc(size*sizeof(int32)); ipnew = (int32 *)ascmalloc(size*sizeof(int32)); ib = (int32 *)ascmalloc(size*sizeof(int32)); row_perm = (int32 *)ascmalloc(size*sizeof(int32)); /* Fill incidence vectors and place included relations w/o * incidence in temporary relation list before the unincluded * and inactive relations */ col_incidence = 0; row_perm_count = 0; for (r=0; r < rlen; r++) { rel = rp[r]; if (rel_apply_filter(rel,&rf) && rel_active(rel)) { len = rel_n_incidences(rel); if (len > 0) { ip[row_perm_count] = col_incidence + 1; /*FORTRAN*/ lenr[row_perm_count] = len; row_perm[row_perm_count] = r; row_perm_count++; list = rel_incidence_list(rel); for (c=0; c < len; c++) { if( var_apply_filter(list[c],&vf) ) { col = var_sindex(list[c]); icn[col_incidence] = col + 1; /*FORTRAN*/ col_incidence++; } } } else { rtmp[rel_tmp_end] = rel; rel_tmp_end--; } } } licn = col_incidence; order = row_perm_count; mc21b(&order,icn,&licn,ip,lenr,iperm,&numnz,iw1,arp,iw2,iw3); if (order == numnz){ FPRINTF(stderr,"they are equal\n"); } else if(order == numnz + 1){ FPRINTF(stderr,"ADD ONE\n"); } else { FPRINTF(stderr,"no relationship\n"); } if (rank == 0 ) { return 1; /* nothing to do, eh? */ } row_high = nrow - 1; col_high = ncol -1; /* lot of whining about dof */ if (rank < nrow) { FPRINTF(stderr,"System is row rank deficient (dependent equations)\n"); row_high = rank - 1; /* KHACK: have found that this is executed erroneously */ } if (rank < ncol) { if ( nrow != rank) { FPRINTF(stderr,"System is row rank deficient with excess columns.\n"); } else { FPRINTF(stderr,"System has degrees of freedom.\n"); } } if (ncol == nrow) { FPRINTF(stderr,"System is (%d) square",ncol); if (ncol != rank) { FPRINTF(stderr,"but rank deficient.\n"); } else { FPRINTF(stderr,".\n"); } } dummy_array_ptr = arp; arp = lenr; lenr = dummy_array_ptr; /* Fix row pointers for assigned relations. Put unassigned onto temporary relation list (before included w/o incidence) and put dummy dense rows into ipnew */ dummy_ptr = licn + 1; /* FORTRAN */ for (row = 0; row < order; row++) { row_ptr = iperm[row] - 1; /* FORTRAN */ if (row_ptr != -1) { /* indicates unassigned row */ ipnew[row_ptr] = ip[row]; lenr[row_ptr] = arp[row]; } else { ipnew[row_ptr] = dummy_ptr; dummy_ptr += ncol; lenr[row_ptr] = ncol; /* unassigned before fixed*/ rtmp[rel_tmp_end] = rp[row_perm[row_ptr]]; rel_tmp_end--; } } for (row = order; row < ncol; row++) { ipnew[row] = dummy_ptr; dummy_ptr += ncol; lenr[row] = ncol; } mc13emod(&size,icn,&licn,ipnew,lenr,arp,ib,&num,iw1,iw2,iw3); /* copy the block list. there is at least one if rank >=1 as above */ len = num; newblocks = (mtx_region_t *)ascmalloc(len*sizeof(mtx_region_t)); if (newblocks == NULL) { return 2; } /* set up locations of block corners with the row and column orderings which will be set soon */ for (c = 0 ; c < len; c++) { newblocks[c].row.low = ib[c]-1; newblocks[c].col.low = newblocks[c].row.low; } for (c = 0 ; c < len -1; c++) { newblocks[c].row.high = newblocks[c+1].row.low -1; newblocks[c].col.high = newblocks[c].row.high; } newblocks[len - 1].row.high = row_high; newblocks[len - 1].col.high = col_high; slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */ d = slv_get_dofdata(sys); d->structural_rank = rank; d->n_rows = nrow; d->n_cols = ncol; /* the next two lines assume only free and fixed incident solvervar * appear in the solvers var list and inequalities are unincluded. */ d->n_fixed = vlen - ncol; d->n_unincluded = rlen - nrow; d->reorder.partition = 1; /* yes */ d->reorder.basis_selection = 0; /* none yet */ d->reorder.block_reordering = 0; /* none */ for (c = 0; c < ncol; c++) { v = arp[c] -1; /* FORTRAN */ vtmp[c] = vp[v]; } for (c = 0; c < vlen; c++) { vp[c] = vtmp[c]; var_set_sindex(vp[c],c); } rel_count = 0; for (c = 0; c < size; c++) { r = arp[c] - 1; if (r < order){ /* weed out fake rows */ r = iperm[r] - 1; if (r != -1) { /* unassigned already in rtmp */ rtmp[rel_count] = rp[row_perm[r]]; rel_count++; } } } for (c = 0; c < rlen; c++) { rp[c] = rtmp[c]; rel_set_sindex(rp[c],c); } ascfree(vtmp); ascfree(rtmp); ascfree(icn); ascfree(ip); ascfree(lenr); ascfree(iperm); ascfree(iw1); ascfree(iw2); ascfree(iw3); ascfree(arp); ascfree(ipnew); ascfree(ib); return 0; } #endif /*STATIC_HARWELL*/ #endif /* 0 */ void slv_sort_rels_and_vars(slv_system_t sys, int32 *rel_count, int32 *var_count) { struct rel_relation **rp; struct rel_relation **rtmp; struct rel_relation *rel; struct var_variable **vp; struct var_variable **vtmp; struct var_variable *var; int32 nrow,ncol,rlen,vlen,order,rel_tmp_end; int32 r,c,var_tmp_end,len; var_filter_t vf; rel_filter_t rf; /* get rel and var info */ rp = slv_get_solvers_rel_list(sys); vp = slv_get_solvers_var_list(sys); rlen = slv_get_num_solvers_rels(sys); vlen = slv_get_num_solvers_vars(sys); if (rlen ==0 || vlen == 0) return; order = MAX(rlen,vlen); assert(var_count != NULL && rel_count != NULL); *var_count = *rel_count = -1; /* allocate temp arrays */ vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *)); if (vtmp == NULL) { return; } rtmp = (struct rel_relation **)ascmalloc(rlen*sizeof(struct rel_relation *)); if (rtmp == NULL) { ascfree(vtmp); return; } /* set up filters */ rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE); rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE); vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE); vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE); /* count rows and cols */ nrow = slv_count_solvers_rels(sys,&rf); ncol = slv_count_solvers_vars(sys,&vf); /* Sort out unincluded and inactive rels to * the end of the temporary relation list. */ *rel_count = 0; rel_tmp_end = rlen -1; for (r=0; r < rlen; r++) { rel = rp[r]; if (rel_apply_filter(rel,&rf)) { rtmp[*rel_count] = rel; (*rel_count)++; } else { rtmp[rel_tmp_end] = rel; rel_tmp_end--; } } /* sort the inactive out of the unincluded and move to end */ len = rlen -1; for(c = rel_tmp_end + 1; len > c; ) { rel = rtmp[c]; if (rel_active(rel)) { c++; } else { rtmp[c] = rtmp[len]; rtmp[len] = rel; len--; } } /* Sort variable list */ *var_count = 0; var_tmp_end = vlen - 1; for (c = 0; c < vlen; c++) { var = vp[c]; if (var_apply_filter(var,&vf)) { vtmp[*var_count] = var; (*var_count)++; } else { vtmp[var_tmp_end] = var; var_tmp_end--; } } /* sort the inactive out of the unincluded and move to end */ len = vlen -1; for(c = var_tmp_end + 1; len > c; ) { var = vtmp[c]; if (var_active(var) && var_active(var)) { c++; } else { vtmp[c] = vtmp[len]; vtmp[len] = var; len--; } } /* reset solver indicies */ for (c = 0; c < vlen; c++) { vp[c] = vtmp[c]; var_set_sindex(vp[c],c); } for (c = 0; c < rlen; c++) { rp[c] = rtmp[c]; rel_set_sindex(rp[c],c); } ascfree(vtmp); ascfree(rtmp); return; } int slv_block_unify(slv_system_t sys) { dof_t *d; const mtx_block_t *mbt; mtx_region_t *newblocks; if (sys==NULL) return 1; d = slv_get_dofdata(sys); mbt = slv_get_solvers_blocks(sys); assert(d!=NULL && mbt!=NULL); if (d->structural_rank && mbt->nblocks >1) { newblocks = (mtx_region_t *)ascmalloc(sizeof(mtx_region_t)); if (newblocks == NULL) return 2; newblocks->row.low = newblocks->col.low = 0; newblocks->row.high = d->structural_rank - 1; newblocks->col.high = d->n_cols - 1; if ( newblocks->col.high < 0 || newblocks->row.high < 0) { ascfree(newblocks); return 1; } slv_set_solvers_blocks(sys,1,newblocks); } return 0; } int slv_set_up_block(slv_system_t sys,int bnum) { struct rel_relation **rp; struct var_variable **vp; mtx_region_t reg; const mtx_block_t *b; int32 c,vlen,rlen; var_filter_t vf; rel_filter_t rf; if (sys==NULL) return 1; rlen = slv_get_num_solvers_rels(sys); vlen = slv_get_num_solvers_vars(sys); if (rlen ==0 || vlen == 0) return 1; rp = slv_get_solvers_rel_list(sys); vp = slv_get_solvers_var_list(sys); assert(rp!=NULL); assert(vp!=NULL); rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE); rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE); vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE); vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE); b = slv_get_solvers_blocks(sys); assert(b!=NULL); if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1; reg = b->block[bnum]; for (c=reg.col.low; c<=reg.col.high; c++) { var_set_in_block(vp[c],1); } for (c=reg.row.low; c<=reg.row.high; c++) { rel_set_in_block(rp[c],1); } return 0; } /* returns 0 if ok, OTHERWISE if madness detected. */ int slv_spk1_reorder_block(slv_system_t sys,int bnum,int transpose) { struct rel_relation **rp; struct var_variable **vp; mtx_region_t reg; const mtx_block_t *b; mtx_matrix_t mtx; mtx_coord_t coord; int32 c,vlen,rlen; var_filter_t vf; rel_filter_t rf; dof_t *d; if (sys==NULL) return 1; rlen = slv_get_num_solvers_rels(sys); vlen = slv_get_num_solvers_vars(sys); if (rlen ==0 || vlen == 0) return 1; rp = slv_get_solvers_rel_list(sys); vp = slv_get_solvers_var_list(sys); assert(rp!=NULL); assert(vp!=NULL); rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE); rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE); vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE); vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE); mtx = mtx_create(); mtx_set_order(mtx,MAX(rlen,vlen)); b = slv_get_solvers_blocks(sys); assert(b!=NULL); if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1; reg = b->block[bnum]; for (c=reg.col.low; c<=reg.col.high; c++) { var_set_in_block(vp[c],1); } for (c=reg.row.low; c<=reg.row.high; c++) { rel_set_in_block(rp[c],1); } if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) { return 1; /* must be square */ /* could also enforce minsize 3x3, but someone my call with a * partitionable region, so don't want to. */ } if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) { FPRINTF(stderr, "slv_spk1_reorder_block: failure in creating incidence matrix.\n"); mtx_destroy(mtx); return 1; } /* verify that block has no empty columns, though not checking diagonal */ for (c = reg.row.low; c <= reg.row.high; c++) { coord.col = mtx_FIRST; coord.row = c; if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) { mtx_destroy(mtx); FPRINTF(stderr, "slv_spk1_reorder_block: empty row (%d) found.\n",c); return 1; } coord.row = mtx_FIRST; coord.col = c; if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) { FPRINTF(stderr, "slv_spk1_reorder_block: empty col (%d) found.\n",c); mtx_destroy(mtx); return 1; } } if (transpose) { mtx_reorder(mtx,®,mtx_TSPK1); } else { mtx_reorder(mtx,®,mtx_SPK1); } if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) { mtx_destroy(mtx); return 2; } if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) { mtx_destroy(mtx); return 2; } d = slv_get_dofdata(sys); d->reorder.block_reordering = 1; /* spk1 */ mtx_destroy(mtx); return 0; } /* global to get around the mr header */ static enum mtx_reorder_method g_blockmethod = mtx_UNKNOWN; /* A function to be called on the leftover diagonal blocks of mr_bisect. This function should be thoroughly investigated, which it has not been. */ static int tear_subreorder(slv_system_t server, mtx_matrix_t mtx, mtx_region_t *reg) { assert(mtx != NULL && reg !=NULL); (void)server; if (g_blockmethod==mtx_UNKNOWN) { mtx_reorder(mtx,reg,mtx_SPK1); } else { mtx_reorder(mtx,reg,g_blockmethod); } return 0; } int slv_tear_drop_reorder_block(slv_system_t sys, int32 bnum, int32 cutoff, int two, enum mtx_reorder_method blockmethod) { struct rel_relation **rp; struct var_variable **vp; const mtx_block_t *b; dof_t *d; mr_reorder_t *mrsys; mtx_matrix_t mtx; mtx_region_t reg; mtx_coord_t coord; int32 c,vlen,rlen,modcount; var_filter_t vf; rel_filter_t rf; if (sys==NULL) return 1; rlen = slv_get_num_solvers_rels(sys); vlen = slv_get_num_solvers_vars(sys); if (rlen ==0 || vlen == 0) return 1; rp = slv_get_solvers_rel_list(sys); vp = slv_get_solvers_var_list(sys); assert(rp!=NULL); assert(vp!=NULL); rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE); rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK| REL_ACTIVE); vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE); vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE); mtx = mtx_create(); mtx_set_order(mtx,MAX(rlen,vlen)); b = slv_get_solvers_blocks(sys); assert(b!=NULL); /* probably shouldn't be an assert here ... */ if (bnum <0 || bnum >= b->nblocks || b->block == NULL) { mtx_destroy(mtx); return 1; } reg = b->block[bnum]; /* do the flag setting even if we don't do the reorder */ for (c=reg.col.low; c<=reg.col.high; c++) { var_set_in_block(vp[c],1); } for (c=reg.row.low; c<=reg.row.high; c++) { rel_set_in_block(rp[c],1); } if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) { mtx_destroy(mtx); return 1; /* must be square */ } if (reg.row.high - reg.row.low < 3) { mtx_destroy(mtx); return 0; /* must be 3x3 or bigger to have any effect */ } if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) { FPRINTF(stderr, "slv_tear_drop_reorder_block: failure in creating incidence matrix.\n"); mtx_destroy(mtx); return 1; } /* verify that block has no empty columns, though not checking diagonal */ for (c = reg.row.low; c <= reg.row.high; c++) { coord.col = mtx_FIRST; coord.row = c; if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) { mtx_destroy(mtx); FPRINTF(stderr, "slv_tear_drop_reorder_block: empty row (%d) found.\n",c); return 1; } coord.row = mtx_FIRST; coord.col = c; if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) { FPRINTF(stderr, "slv_tear_drop_reorder_block: empty col (%d) found.\n",c); mtx_destroy(mtx); return 1; } } modcount = slv_get_num_models(sys); mrsys = mr_reorder_create(sys,mtx,modcount); if (mrsys==NULL) return 1; mrsys->cutoff = cutoff; g_blockmethod = blockmethod; if (!two) { mr_bisect_partition(mrsys,®,1,tear_subreorder); } else { mr_bisect_partition2(mrsys,®,1,tear_subreorder); } #if 1 FPRINTF(stderr,"Model-based reordering: (tear-style = %d)\n",two); FPRINTF(stderr,"Min. tearable size:\t%d\n",mrsys->cutoff); FPRINTF(stderr,"Tears:\t\t%d\n",mrsys->ntears); FPRINTF(stderr,"Partitionings:\t%d\n",mrsys->nblts); #endif reg = b->block[bnum]; /* bisect likely munged reg */ if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) { mtx_destroy(mtx); mr_reorder_destroy(mrsys); return 2; } if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) { mtx_destroy(mtx); mr_reorder_destroy(mrsys); return 2; } d = slv_get_dofdata(sys); d->reorder.block_reordering = 2; /* tear_drop_baa */ mtx_destroy(mtx); mr_reorder_destroy(mrsys); return 0; } int slv_insure_bounds(slv_system_t sys,int32 lo,int32 hi, FILE *mif) { real64 val,low,high; int32 c,nchange=0; struct var_variable *var, **vp; vp = slv_get_solvers_var_list(sys); if (vp==NULL) return -1; for (c= lo; c <= hi; c++) { var = vp[c]; low = var_lower_bound(var); high = var_upper_bound(var); val = var_value(var); if( low > high ) { if (mif!=NULL) { FPRINTF(mif,"Bounds for variable "); var_write_name(sys,var,mif); FPRINTF(mif," are inconsistent [%g,%g].\n",low,high); FPRINTF(mif,"Bounds will be swapped.\n"); } var_set_upper_bound(var, low); var_set_lower_bound(var, high); low = var_lower_bound(var); high = var_upper_bound(var); nchange++; } if( low > val ) { if (mif!=NULL) { FPRINTF(mif,"Variable "); var_write_name(sys,var,mif); FPRINTF(mif," was initialized below its lower bound.\n"); FPRINTF(mif,"It will be moved to its lower bound.\n"); } var_set_value(var, low); } else { if( val > high ) { if (mif!=NULL) { FPRINTF(mif,"Variable "); var_write_name(sys,var,mif); FPRINTF(mif," was initialized above its upper bound.\n"); FPRINTF(mif,"It will be moved to its upper bound.\n"); } var_set_value(var, high); nchange++; } } } return nchange; } void slv_check_bounds(const slv_system_t sys,int32 lo,int32 hi, FILE *mif,const char *label) { real64 val,low,high; int32 c,len; struct var_variable *var, **vp; static char defaultlabel[] = ""; if (label==NULL) label = defaultlabel; if (sys==NULL || mif == NULL) return; vp = slv_get_solvers_var_list(sys); if (vp==NULL) return; len = slv_get_num_solvers_vars(sys); if (lo > len || hi > len) { FPRINTF(stderr,"slv_check_bounds miscalled\n"); return; } for (c= lo; c <= hi; c++) { var = vp[c]; low = var_lower_bound(var); high = var_upper_bound(var); val = var_value(var); if( low > high ) { FPRINTF(mif,"Bounds for %s variable ",label); var_write_name(sys,var,mif); FPRINTF(mif,"\nare inconsistent [%g,%g].\n",low,high); } if( low > val ) { FPRINTF(mif,"%s variable ",label); var_write_name(sys,var,mif); FPRINTF(mif,"\nwas initialized below its lower bound.\n"); } else { if( val > high ) { FPRINTF(mif,"%s variable ",label); var_write_name(sys,var,mif); FPRINTF(mif," was initialized above its upper bound.\n"); } } } return; } int SlvRegisterStandardClients(void) { int nclients = 0; int status; #ifdef STATIC_SLV status = slv_register_client(slv0_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv0.\n"); } else { nclients++; } #endif #ifdef STATIC_MINOS status = slv_register_client(slv1_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv1.\n"); } else { nclients++; } #endif #ifdef STATIC_QRSLV status = slv_register_client(slv3_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv3.\n"); } else { nclients++; } #endif #ifdef STATIC_CSLV status = slv_register_client(slv4_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv4.\n"); } else { nclients++; } #endif #ifdef STATIC_LSSLV status = slv_register_client(slv5_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv5.\n"); } else { nclients++; } #endif #ifdef STATIC_MPS status = slv_register_client(slv6_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv6.\n"); } else { nclients++; } #endif #ifdef STATIC_NGSLV status = slv_register_client(slv7_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv7.\n"); } else { nclients++; } #endif #ifdef STATIC_OPTSQP status = slv_register_client(slv2_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv2.\n"); } else { nclients++; } #endif #if (defined(STATIC_CONOPT) || defined(DYNAMIC_CONOPT)) status = slv_register_client(slv8_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv8.\n"); } else { nclients++; } #endif #ifdef STATIC_CMSLV status = slv_register_client(slv9_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv9.\n"); } else { nclients++; } #endif #ifdef STATIC_LRSLV status = slv_register_client(slv9a_register,NULL,NULL); if (status) { FPRINTF(stderr,"Unable to register slv9a.\n"); } else { nclients++; } #endif return nclients; } /************************************************************************\ Output Assignment and partitiong in Logical Relations \************************************************************************/ #define MLIMDEBUG 0 /* * Populates a matrix according to the sys solvers_dvars, solvers_logrels * lists and the filters given. mtx given must be created (not null) * and with order >= MAX( slv_get_num_solvers_logrels(sys), * slv_get_num_solvers_dvars(sys)); * returns 0 if went ok. */ int slv_std_make_log_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx, dis_filter_t *dvf,logrel_filter_t *lrf) { #if MLIMDEBUG FILE *fp; #endif int32 lr, lrlen,ord; struct logrel_relation **lrp; if (sys==NULL || mtx == NULL || dvf == NULL || lrf == NULL) { FPRINTF(stderr,"make_log_incidence called with null\n"); return 1; } lrp = slv_get_solvers_logrel_list(sys); assert(lrp!=NULL); lrlen = slv_get_num_solvers_logrels(sys); ord = MAX(lrlen,slv_get_num_solvers_dvars(sys)); if (ord > mtx_order(mtx)) { FPRINTF(stderr,"make_incidence called with undersized matrix\n"); return 2; } for (lr=0; lr < lrlen; lr++) { if (logrel_apply_filter(lrp[lr],lrf)) { logrelman_get_incidence(lrp[lr],dvf,mtx); } } #if MLIMDEBUG fp = fopen("/tmp/mim.plot","w+"); if (fp !=NULL) { mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX); fclose(fp); } #endif return 0; } #if USEDCODE /* returns 0 if successful, 1 if insufficient memory. Does not change * the data in mtx. Orders the solvers_dvar list of the system to * match the permutation on the given mtx. It is assumed that * the org cols of mtx == dvar list position. Only the * dvars in range lo to hi of the dvar list are affected and * only these dvars should appear in the same org column range of * the mtx. This should not be called on blocks less than 3x3. */ static int reindex_dvars_from_mtx(slv_system_t sys, int32 lo, int32 hi, const mtx_matrix_t mtx) { struct dis_discrete **dvtmp, **dvp; int32 c,v,vlen; if (lo >= hi +1) return 0; /* job too small */ dvp = slv_get_solvers_dvar_list(sys); vlen = slv_get_num_solvers_dvars(sys); /* on dvtmp we DONT have the terminating null */ dvtmp = (struct dis_discrete **) ascmalloc(vlen*sizeof(struct dis_discrete *)); if (dvtmp == NULL) { return 1; } /* copy pointers to dvtmp in order desired and copy back rather than sort */ for (c=lo;c<=hi;c++) { v = mtx_col_to_org(mtx,c); dvtmp[c] = dvp[v]; } /* copying back and re-sindexing */ for (c=lo;c<=hi;c++) { dvp[c] = dvtmp[c]; dis_set_sindex(dvp[c],c); } ascfree(dvtmp); return 0; } /* returns 0 if successful, 1 if insufficient memory. Does not change * the data in mtx. Orders the solvers_logrel list of the system to * match the permutation on the given mtx. It is assumed that * the org rows of mtx == logrel list position. Only the *logrels in range lo to hi of the logrel list are affected and * only these logrels should appear in the same org row range of * the input mtx. This should not be called on blocks less than 3x3. */ static int reindex_logrels_from_mtx(slv_system_t sys, int32 lo, int32 hi, const mtx_matrix_t mtx) { struct logrel_relation **lrtmp, **lrp; int32 c,v,rlen; lrp = slv_get_solvers_logrel_list(sys); rlen = slv_get_num_solvers_logrels(sys); /* on lrtmp we DONT have the terminating null */ lrtmp = (struct logrel_relation **) ascmalloc(rlen*sizeof(struct logrel_relation *)); if (lrtmp == NULL) { return 1; } /* copy pointers to lrtmp in order desired and copy back rather than sort. * do this only in the row range of interest. */ for (c=lo;c<=hi;c++) { v = mtx_row_to_org(mtx,c); #if RIDEBUG FPRINTF(stderr,"Old logrel sindex(org) %d becoming sindex(cur) %d\n",v,c); #endif lrtmp[c] = lrp[v]; } /* copying back and re-sindexing */ for (c=lo;c<=hi;c++) { lrp[c] = lrtmp[c]; logrel_set_sindex(lrp[c],c); } ascfree(lrtmp); return 0; } #endif /* USEDCODE */ /* * returns 0 if ok, OTHERWISE if madness detected. */ #define SLBPDEBUG 0 int slv_log_block_partition(slv_system_t sys) { #if SLBPDEBUG FILE *fp; #endif struct logrel_relation **lrp; struct dis_discrete **dvp; mtx_region_t *newblocks; dof_t *d; int32 nrow,ncol; int32 ncolpfix, nrowpun, vlenmnv; mtx_matrix_t mtx; int32 order, rank; int32 c,len,dvlen,r,lrlen; dis_filter_t dvf; logrel_filter_t lrf; lrp = slv_get_solvers_logrel_list(sys); dvp = slv_get_solvers_dvar_list(sys); lrlen = slv_get_num_solvers_logrels(sys); dvlen = slv_get_num_solvers_dvars(sys); if (lrlen ==0 || dvlen == 0) return 1; order = MAX(lrlen,dvlen); lrf.matchbits = (LOGREL_ACTIVE); lrf.matchvalue = (LOGREL_ACTIVE); dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE); dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE); /* nrow plus unincluded logrels */ nrowpun = slv_count_solvers_logrels(sys,&lrf); /* ncol plus fixed discrete vars */ ncolpfix = slv_count_solvers_dvars(sys,&dvf); dvf.matchbits = (DIS_BVAR); dvf.matchvalue = 0; /* dvlen minus non boolean vars */ vlenmnv = dvlen - slv_count_solvers_dvars(sys,&dvf); lrf.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE); lrf.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE); dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_FIXED | DIS_ACTIVE); dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE); nrow = slv_count_solvers_logrels(sys,&lrf); ncol = slv_count_solvers_dvars(sys,&dvf); mtx = slv_get_sys_mtx(sys); mtx_set_order(mtx,order); if (slv_std_make_log_incidence_mtx(sys,mtx,&dvf,&lrf)) { FPRINTF(stderr, "slv_log_block_partition: failure in creating incidence matrix.\n"); mtx_destroy(mtx); return 1; } mtx_output_assign(mtx,lrlen,dvlen); rank = mtx_symbolic_rank(mtx); if (rank == 0 ) return 1; if (rank < nrow) { FPRINTF(stderr,"rank=1 as above */ len = mtx_number_of_blocks(mtx); newblocks = (mtx_region_t *)ascmalloc(len*sizeof(mtx_region_t)); if (newblocks == NULL) { mtx_destroy(mtx); return 2; } for (c = 0 ; c < len; c++) { mtx_block(mtx,c,&(newblocks[c])); } slv_set_solvers_log_blocks(sys,len,newblocks); /* give it to the system */ d = slv_get_log_dofdata(sys); d->structural_rank = rank; d->n_rows = nrow; d->n_cols = ncol; d->n_fixed = ncolpfix - ncol; d->n_unincluded = nrowpun - nrow; d->reorder.partition = 1; /* yes */ d->reorder.basis_selection = 0; /* none yet */ d->reorder.block_reordering = 0; /* none */ #if SLBPDEBUG fp = fopen("/tmp/sbp1.plot","w+"); if (fp !=NULL) { mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX); fclose(fp); } #endif len = vlenmnv - 1; /* row of last possible boolean variable */ for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */ r = mtx_col_to_org(mtx,c); if ( dis_fixed(dvp[r]) && dis_active(dvp[r])) { c++; } else { mtx_swap_cols(mtx,c,len); len--; } } /* Now, leftover columns are where we want them. * That is, unassigned incidence is next to assigned and fixed on right with * no columns which do not correspond to dis vars in the range 0..dvlen-1. */ /* there aren't any nonincident dvars in the solvers dvarlist, so * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */ /* rows 0->lrlen-1 should be actual logrels, though, and not phantoms. * At this point we should have * rows - assigned * - unassigned w/incidence * - unincluded and included w/out free incidence and inactive mixed. */ /* sort unincluded logrelations to bottom of good region */ /* this is needed because the idiot user may have fixed everything */ /* in a logrelation, but we want those rows not to be confused with */ /* unincluded rows.*/ len = lrlen - 1; /* row of last possible relation */ /* sort the inactive out of the unassigned/unincluded */ for (c=rank; len > c ; ) { r = mtx_row_to_org(mtx,c); if ( logrel_active(lrp[r])) { c++; } else { mtx_swap_rows(mtx,c,len); len--; } } /* sort the unincluded out of the unassigned */ len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */ for (c=rank; len > c ; ) { r = mtx_row_to_org(mtx,c); if ( logrel_included(lrp[r]) ) { c++; } else { mtx_swap_rows(mtx,c,len); len--; } } #if SBPDEBUG fp = fopen("/tmp/sbp2.plot","w+"); if (fp !=NULL) { mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX); fclose(fp); } #endif /* These functions are called for slv_block_partitioning, but I do not * know if I need them here: * * now, at last we have cols in the order we want the lists to * be handed to the solver. So, we need to reset the dvar sindex values * and reorder the lists. */ #if USEDCODE if (reindex_dvars_from_mtx(sys,0,dvlen-1,mtx)) { mtx_destroy(mtx); return 2; } #endif /* USEDCODE */ /* * now, at last we have rows jacobian in the order we want the lists to * be handed to the solvers. So, we need to reset the rel sindex values. */ #if USEDCODE if (reindex_logrels_from_mtx(sys,0,lrlen-1,mtx)) { mtx_destroy(mtx); return 2; } #endif /* USEDCODE */ return 0; }