DSDP
|
00001 #include "numchol.h" 00002 00003 int OdAlloc(int nnod, 00004 int nn0, 00005 char *info, 00006 order **rr) 00007 { 00008 order *r; 00009 int ierr=0; 00010 00011 r=(order*)calloc(1,sizeof(order)); 00012 if (!r) ExitProc(OutOfSpc,info); 00013 00014 r->nnod=nnod; 00015 r->nn0 =nn0; 00016 00017 ierr=iAlloc(nn0,info,&r->adjn); if(ierr) return 1; 00018 ierr=iAlloc(nnod,info,&r->rbeg); if(ierr) return 1; 00019 ierr=iAlloc(nnod,info,&r->rexs); if(ierr) return 1; 00020 ierr=iAlloc(nnod,info,&r->rlen); if(ierr) return 1; 00021 ierr=iAlloc(nnod,info,&r->rend); if(ierr) return 1; 00022 ierr=iAlloc(nnod,info,&r->pres); if(ierr) return 1; 00023 ierr=iAlloc(nnod,info,&r->succ); if(ierr) return 1; 00024 *rr=r; 00025 return (0); 00026 } /* OdAlloc */ 00027 00028 void OdFree(order **od) 00029 { 00030 order *r; 00031 00032 if (*od) { 00033 r=*od; 00034 iFree(&r->adjn); 00035 iFree(&r->rbeg); 00036 iFree(&r->rexs); 00037 iFree(&r->rlen); 00038 iFree(&r->rend); 00039 iFree(&r->pres); 00040 iFree(&r->succ); 00041 free(*od); 00042 *od=NULL; 00043 } 00044 } /* OdFree */ 00045 00046 void OdInit(order *od, 00047 int *nnzi) 00048 { 00049 int i,n=od->nnod; 00050 00051 if (n) { 00052 od->rexs[0]=nnzi[0]; 00053 od->rlen[0]=nnzi[0]; 00054 od->rbeg[0]=0; 00055 od->pres[0]=n; 00056 od->succ[0]=1; 00057 for(i=1; i<od->nnod; ++i) { 00058 od->pres[i]=i-1; 00059 od->succ[i]=i+1; 00060 od->rexs[i]=nnzi[i]; 00061 od->rlen[i]=nnzi[i]; 00062 od->rbeg[i]=od->rbeg[i-1]+od->rlen[i-1]; 00063 } 00064 00065 od->succ[n-1]=n; 00066 od->last =n-1; 00067 00068 od->raft=od->rbeg[n-1]+nnzi[n-1]; 00069 00070 if (od->raft>od->nn0) 00071 ExitProc(OutOfSpc,"InitMmd"); 00072 } 00073 } /* OdInit */ 00074 00075 static void OdSet(order *od, 00076 int allow_eli, 00077 xlist *elist, 00078 int *node_status, 00079 int *marker, 00080 int *isize, 00081 int *ilink, 00082 int *oinfo, 00083 int *osize, 00084 int *e, 00085 int *p) 00086 { 00087 int i,n,deg,*rbeg,*rexs,*rend; 00088 00089 n =od->nnod; 00090 rbeg=od->rbeg; 00091 rexs=od->rexs; 00092 rend=od->rend; 00093 *e =0; 00094 00095 for (i=0; i<n; i++) { 00096 isize[i] =0; 00097 ilink[i] =n; 00098 osize[i] =0; 00099 oinfo[i] =n; 00100 marker[i]=0; 00101 } 00102 00103 for(i=0; i<n; ++i) { 00104 rbeg[i] -= rexs[i]; 00105 rend[i] = 0; 00106 } 00107 00108 for(i=0; i<n; ++i) { 00109 deg = rexs[i]; 00110 if (!allow_eli||deg) { 00111 node_status[i]=1; 00112 XtPut(elist,i,deg); 00113 } 00114 else { 00115 node_status[i] = 0; 00116 marker[i] = TRUE; 00117 p[*e] = i; 00118 (*e)++; 00119 } 00120 } 00121 } /* OdSet */ 00122 00123 void OdIndex(order *od, 00124 int i, 00125 int j) 00126 { 00127 if (i!=j) { 00128 od->adjn[od->rbeg[i]++]=j; 00129 od->adjn[od->rbeg[j]++]=i; 00130 } 00131 } /* OdIndex */ 00132 00133 static void OdArriv(order *od, 00134 int *node_status, 00135 int *marker, 00136 int *isize, 00137 int x, 00138 int *xdeg, 00139 int *rsze, 00140 int *esze, 00141 int *rchset) 00142 { 00143 int *visited,i,n,y,z,l,s,t,f,stopt, 00144 stops,*adjn,*rbeg,*rexs,*rend; 00145 00146 n =od->nnod; 00147 adjn =od->adjn; 00148 rbeg =od->rbeg; 00149 rexs =od->rexs; 00150 rend =od->rend; 00151 *rsze=0; 00152 *esze=0; 00153 00154 if (rexs[x]) { 00155 l=n; 00156 00157 visited=marker; 00158 visited[x]=TRUE; 00159 00160 for(t=rbeg[x], stopt=rbeg[x]+rend[x]; t<stopt; ++t) { 00161 y=adjn[t]; 00162 00163 if (node_status[y]!=0) { 00164 l--; 00165 rchset[l]=y; 00166 visited[y]=TRUE; 00167 00168 for(s=rbeg[y], stops=rbeg[y]+rexs[y]; s<stops; ++s) { 00169 z=adjn[s]; 00170 00171 if (node_status[z]!=0) { 00172 if (!visited[z]) { 00173 visited[z]=TRUE; 00174 00175 rchset[*rsze]=z; 00176 (*rsze)++; 00177 } 00178 } 00179 } 00180 } 00181 } 00182 00183 f=rbeg[x]+rend[x]; 00184 for(t=f, stopt=rbeg[x]+rexs[x]; t<stopt; ++t) { 00185 y=adjn[t]; 00186 if (!visited[y]) { 00187 adjn[f++]=y; 00188 visited[y]=TRUE; 00189 00190 rchset[*rsze]=y; 00191 (*rsze)++; 00192 } 00193 } 00194 00195 rexs[x]=f-rbeg[x]; 00196 00197 *esze=n-l; 00198 visited[x] = FALSE; 00199 iZero(*rsze,visited,rchset); 00200 iZero(n-l,visited,rchset+l); 00201 } 00202 00203 if (xdeg) { 00204 *xdeg = *rsze+isize[x]; 00205 for(i=0; i<*rsze; ++i) 00206 *xdeg+=isize[rchset[i]]; 00207 } 00208 } /* OdArriv */ 00209 00210 static void OdRenew(order *od, 00211 int *ilink, 00212 int x, 00213 int xdeg, 00214 int *e, 00215 int *p) 00216 { 00217 int c,n; 00218 00219 n =od->nnod; 00220 od->ntot+=xdeg--; 00221 p[*e] =x; 00222 (*e)++; 00223 for(c=x; ilink[c]!=n; c=ilink[c]) { 00224 od->ntot+=xdeg--; 00225 p[*e] =ilink[c]; 00226 (*e)++; 00227 } 00228 } /* OdRenew */ 00229 00230 static void OdCheck(order *od, 00231 int *node_status) 00232 { 00233 int f,i,t,stopt,rnew,z,previous,n,*adjn, 00234 *rbeg,*rexs,*rlen,*rend,*pres,*succ; 00235 00236 n =od->nnod; 00237 adjn=od->adjn; 00238 rbeg=od->rbeg; 00239 rexs=od->rexs; 00240 rlen=od->rlen; 00241 rend=od->rend; 00242 pres=od->pres; 00243 succ=od->succ; 00244 00245 f=0; 00246 previous=n; 00247 for(i=od->head; i!=n; i=succ[i]) { 00248 if (node_status[i]!=0) { 00249 rnew=f; 00250 00251 for(t=rbeg[i], stopt=rbeg[i]+rend[i]; t<stopt; ++t) { 00252 z=adjn[t]; 00253 if (node_status[z]==3) 00254 adjn[f++]=z; 00255 } 00256 00257 rend[i]=f-rnew; 00258 00259 for(stopt=rbeg[i]+rexs[i]; t<stopt; ++t) { 00260 z=adjn[t]; 00261 if (node_status[z]!=0) 00262 adjn[f++]=z; 00263 } 00264 00265 rexs[i]=f-rnew; 00266 rlen[i]=rexs[i]; 00267 00268 rbeg[i]=rnew; 00269 00270 if (previous==n) { 00271 od->head=i; 00272 pres[i]=n; 00273 } 00274 00275 else { 00276 succ[previous]=i; 00277 pres[i]=previous; 00278 } 00279 previous=i; 00280 } 00281 } 00282 00283 if (previous!=n) { 00284 succ[previous]=n; 00285 od->raft=rbeg[previous]+rexs[previous]; 00286 } 00287 00288 od->last=previous; 00289 } /* OdCheck */ 00290 00291 static void OdAdd(order *od, 00292 int *node_status, 00293 int x, 00294 int newsze) 00295 { 00296 int n,*adjn,*rbeg,*rexs,*rlen,*pres,*succ; 00297 00298 n =od->nnod; 00299 adjn=od->adjn; 00300 rbeg=od->rbeg; 00301 rexs=od->rexs; 00302 rlen=od->rlen; 00303 pres=od->pres; 00304 succ=od->succ; 00305 00306 if (newsze<=rlen[x]) 00307 return; 00308 00309 if (od->raft+newsze>od->nn0) 00310 OdCheck(od,node_status); 00311 00312 if (od->raft+newsze>od->nn0) 00313 ExitProc(OutOfSpc,"OdAdd"); 00314 00315 if (pres[x]!=n) 00316 rlen[pres[x]]+=rlen[x]; 00317 00318 iCopy(rexs[x],adjn+rbeg[x],adjn+od->raft); 00319 rbeg[x]=od->raft; 00320 rlen[x]=newsze; 00321 od->raft+=newsze; 00322 00323 if (pres[x]==n) { 00324 if (succ[x]==n) 00325 od->head=x; 00326 else 00327 od->head=succ[x]; 00328 } 00329 00330 else { 00331 if (succ[x]==n) 00332 succ[pres[x]]=x; 00333 else 00334 succ[pres[x]]=succ[x]; 00335 } 00336 00337 if (succ[x]!=n) 00338 pres[succ[x]]=pres[x]; 00339 00340 if (od->last!=x) { 00341 succ[od->last]=x; 00342 pres[x]=od->last; 00343 } 00344 00345 succ[x] =n; 00346 od->last=x; 00347 } /* OdAdd */ 00348 00349 static int OdComb(order *od, 00350 int *node_status, 00351 int *marker, 00352 int *isize, 00353 int *ilink, 00354 int *osize, 00355 int xsize, 00356 int *xset) 00357 { 00358 int i,n,rnew,rlen,x,icur; 00359 00360 n =od->nnod; 00361 rlen=0; 00362 00363 if (xsize==0) 00364 rnew=n; 00365 else if (xsize==1) 00366 rnew=xset[0]; 00367 else { 00368 rnew=xset[0]; 00369 for(i=1; i<xsize; ++i) 00370 rlen+=1+isize[xset[i]]; 00371 00372 node_status[rnew]=1; 00373 osize[rnew]=0; 00374 00375 for(icur=rnew; ilink[icur]!=n; icur=ilink[icur]); 00376 isize[rnew]+=rlen; 00377 00378 for(i=1; i<xsize; ++i) { 00379 x=xset[i]; 00380 00381 node_status[x]=0; 00382 marker[x]=TRUE; 00383 00384 ilink[icur]=x; 00385 00386 for(icur=x; ilink[icur]!=n; icur=ilink[icur]); 00387 00388 isize[x]=0; 00389 } 00390 } 00391 00392 return (rnew); 00393 } /* OdComb */ 00394 00395 static int OdSelect(order *od, 00396 xlist *elist, 00397 int *node_status, 00398 int *marker, 00399 int *isize, 00400 int *ilink, 00401 int *oinfo, 00402 int *osize, 00403 int x, 00404 int *rsze, 00405 int *rchset, 00406 int *ibuf1, 00407 int *ibuf2, 00408 int *mask2, 00409 int *e, 00410 int *p) 00411 { 00412 int absorp,old,i,j,n,esze,y,z,l,f,t,stopt,s, 00413 o,stops,indsze,xdeg,e0,ssze,*slist,tsze, 00414 *tlist,sze,*adjn,*rbeg,*rexs,*rlen,*rend; 00415 00416 adjn =od->adjn; 00417 rbeg =od->rbeg; 00418 rexs =od->rexs; 00419 rlen =od->rlen; 00420 rend =od->rend; 00421 n =od->nnod; 00422 slist=ibuf1; 00423 00424 e0 = *e; 00425 OdArriv(od,node_status,marker,isize,x,&xdeg,rsze,&esze,rchset); 00426 00427 XtDel(elist,x); 00428 00429 OdRenew(od,ilink,x,xdeg,e,p); 00430 00431 for(i=n-esze; i<n; ++i) { 00432 node_status[rchset[i]]=0; 00433 marker[rchset[i]]=TRUE; 00434 } 00435 00436 marker[x]=TRUE; 00437 iSet(*rsze,TRUE,marker,rchset); 00438 00439 ssze=0; 00440 for(i=0; i<*rsze;) { 00441 y=rchset[i]; 00442 00443 if (node_status[y]==0||node_status[y]==3) 00444 ExitProc(SysError,NULL); 00445 00446 f=rbeg[y]; 00447 for(t=f, stopt=f+rend[y]; t<stopt; ++t) { 00448 z=adjn[t]; 00449 if (node_status[z]==3) { 00450 adjn[f++]=z; 00451 00452 if (!mask2[z]) { 00453 slist[ssze++]=z; 00454 mask2[z]=TRUE; 00455 } 00456 } 00457 } 00458 rend[y]=f-rbeg[y]; 00459 00460 for(stopt=rbeg[y]+rexs[y]; t<stopt; ++t) { 00461 z=adjn[t]; 00462 if (!marker[z]) 00463 adjn[f++]=z; 00464 } 00465 00466 rexs[y]=f-rbeg[y]; 00467 00468 if (rexs[y]==0) { 00469 OdRenew(od,ilink,y,xdeg-(*e-e0),e,p); 00470 node_status[y] = 0; 00471 marker[y] = TRUE; 00472 00473 (*rsze)--; 00474 iSwap(i,*rsze,rchset); 00475 } 00476 00477 else { 00478 if (rexs[y]>=rlen[y]) 00479 ExitProc(SysError,NULL); 00480 00481 if (rexs[y]>rend[y]) 00482 adjn[rbeg[y]+rexs[y]]=adjn[rbeg[y]+rend[y]]; 00483 00484 rexs[y]++; 00485 00486 adjn[rbeg[y]+rend[y]]=x; 00487 rend[y]++; 00488 00489 i++; 00490 } 00491 } 00492 00493 iSet(ssze,FALSE,mask2,slist); 00494 00495 if (*rsze==0) { 00496 node_status[x]=0; 00497 marker[x]=TRUE; 00498 } 00499 00500 else { 00501 node_status[x]=3; 00502 00503 rend[x]=0; 00504 rexs[x]=0; 00505 if (*rsze>rlen[x]) 00506 OdAdd(od,node_status,x,*rsze); 00507 00508 rexs[x]=*rsze; 00509 iCopy(*rsze,rchset,adjn+rbeg[x]); 00510 00511 tsze=0; 00512 tlist=ibuf2; 00513 for(i=0; i<ssze; ++i){ 00514 y=slist[i]; 00515 old=marker[y]; 00516 marker[y]=TRUE; 00517 00518 absorp=TRUE; 00519 00520 indsze=n; 00521 l=n; 00522 00523 f=rbeg[y]; 00524 for(t=f, stopt=f+rexs[y]; t<stopt; ++t) { 00525 z=adjn[t]; 00526 if (node_status[z]!=0) { 00527 adjn[f++]=z; 00528 00529 if (marker[z]) { 00530 l--; 00531 slist[l]=z; 00532 00533 if (!mask2[z]) { 00534 for(s=rbeg[z],stops=rbeg[z]+rexs[z]; 00535 s<stops &&marker[adjn[s]]; ++s); 00536 00537 if (s==stops) { 00538 indsze--; 00539 iSwap(l,indsze,slist); 00540 } 00541 00542 mask2[z]=TRUE; 00543 tlist[tsze++]=z; 00544 } 00545 } 00546 else 00547 absorp=FALSE; 00548 } 00549 } 00550 00551 marker[y]=old; 00552 rexs[y]=f-rbeg[y]; 00553 00554 if (indsze<n) { 00555 z=OdComb(od,node_status,marker, 00556 isize,ilink,osize, 00557 n-indsze,slist+indsze); 00558 00559 node_status[z]=1; 00560 00561 sze=0; 00562 for(j=l; j<indsze; ++j) { 00563 o=slist[j]; 00564 sze+=1+isize[o]; 00565 node_status[o]=2; 00566 oinfo[o]=z; 00567 } 00568 osize[z]=max(osize[z],sze); 00569 } 00570 00571 if (absorp) { 00572 node_status[y]=0; 00573 marker[y]=TRUE; 00574 } 00575 } 00576 00577 iSet(tsze,FALSE,mask2,tlist); 00578 } 00579 00580 marker[x]=(node_status[x]==0); 00581 00582 for(t=0; t<*rsze; ++t) { 00583 z=rchset[t]; 00584 marker[z]=(node_status[z]==0); 00585 } 00586 00587 return (FALSE); 00588 } /* OdSelect */ 00589 00590 static int OdOrder(order *od, 00591 int *node_status, 00592 int *marker, 00593 int *isize, 00594 int x, 00595 int *ibuf1) 00596 { 00597 int rsze,esze,deg; 00598 00599 OdArriv(od,node_status,marker,isize, 00600 x,°,&rsze,&esze,ibuf1); 00601 00602 return deg; 00603 } /* OdOrder */ 00604 00605 static void OdModf(order *od, 00606 xlist *elist, 00607 int *node_status, 00608 int *marker, 00609 int *isize, 00610 int *oinfo, 00611 int rsze, 00612 int *rchset, 00613 int *ibuf1) 00614 { 00615 00616 int i,x,deg; 00617 00618 for(i=0; i<rsze; ++i) { 00619 x=rchset[i]; 00620 if (node_status[x]==2) 00621 if (node_status[oinfo[x]]==0||node_status[oinfo[x]]==3) 00622 node_status[x]=1; 00623 00624 if (node_status[x]==1) { 00625 deg=OdOrder(od,node_status,marker,isize,x,ibuf1); 00626 XtPut(elist,x,deg-isize[x]); 00627 } 00628 00629 else 00630 XtDel(elist,x); 00631 } 00632 } /* mindeg_upddeg */ 00633 00634 void OdProc(order *od, 00635 xlist *xt, 00636 int *bbuf1, 00637 int *bbuf2, 00638 int *ibuf1, 00639 int *ibuf2, 00640 int *ibuf3, 00641 int *ibuf4, 00642 int *ibuf5, 00643 int *ibuf6, 00644 int *ibuf7, 00645 int *ibuf8, 00646 int *ibuf9, 00647 int *intbuf1, 00648 int *p) 00649 { 00650 int *mask2,total_fillin,use_mtmd=TRUE,*marker=bbuf1, 00651 *node_status,i,n,e,x,y,rsze,deg,mindeg,xsize,sze, 00652 j,*isize,*ilink,*oinfo,*osize,*rchset,*dependent, 00653 *slist; 00654 xlist *elist; 00655 00656 elist=xt; 00657 00658 isize =ibuf1; 00659 ilink =ibuf2; 00660 oinfo =ibuf3; 00661 osize =ibuf4; 00662 rchset =ibuf5; 00663 dependent=ibuf6; 00664 slist =ibuf7; 00665 00666 node_status=intbuf1; 00667 mask2=bbuf2; 00668 00669 OdSet(od,TRUE,elist,node_status,marker, 00670 isize,ilink,oinfo,osize,&e,p); 00671 00672 n=od->nnod; 00673 00674 iSet(n,0,dependent,NULL); 00675 00676 total_fillin=FALSE; 00677 for(; e<n && !total_fillin;) { 00678 00679 XtLeast(elist); 00680 00681 if (!XtGet(elist,&y,&mindeg)) { 00682 printf("\n No new nodes e=%d n=%d",e,n); 00683 printf(" Node status: "); 00684 00685 for(i=0; i<n; ++i) 00686 if (node_status[i]==1) 00687 printf("A\n"); 00688 else if (node_status[i]==2) 00689 printf("\n O%d: rlen=%d oinfo=%d\n", 00690 i,isize[i],oinfo[i]); 00691 00692 ExitProc(SysError,NULL); 00693 } 00694 00695 xsize=0; 00696 for(;use_mtmd;) { 00697 if (!XtGet(elist,&x,°)||deg>mindeg) 00698 break; 00699 00700 if (node_status[x]!=1) 00701 XtDel(elist,x); 00702 00703 else { 00704 XtSucc(elist); 00705 if (!dependent[x]) { 00706 00707 total_fillin = OdSelect(od,elist,node_status,marker, 00708 isize,ilink,oinfo,osize,x, 00709 &rsze,rchset,ibuf8,ibuf9,mask2, 00710 &e,p); 00711 00712 if (!total_fillin) { 00713 dependent[x]=2; 00714 slist[xsize++]=x; 00715 00716 for(i=0; i<rsze; ++i) { 00717 y=rchset[i]; 00718 if (!dependent[y]) { 00719 dependent[y]=1; 00720 slist[xsize++]=y; 00721 } 00722 } 00723 } 00724 00725 if (!use_mtmd) break; 00726 } 00727 } 00728 } 00729 00730 if (!total_fillin) { 00731 sze=0; 00732 for(j=0; j<xsize; ++j) { 00733 y=slist[j]; 00734 if (dependent[y]==1 && node_status[y]!=0) 00735 slist[sze++]=y; 00736 dependent[y]=0; 00737 } 00738 00739 OdModf(od,elist,node_status,marker, 00740 isize,oinfo,sze,slist,ibuf8); 00741 00742 } 00743 } 00744 00745 if (e<n) { 00746 sze=0; 00747 for(i=0; i<n; ++i) 00748 if (node_status[i]==2||node_status[i]==1) 00749 ibuf8[sze++]=i; 00750 00751 x = OdComb(od,node_status,marker,isize,ilink,osize, 00752 sze,ibuf8); 00753 00754 OdRenew(od,ilink,x,n-e-1,&e,p); 00755 } 00756 00757 } /* OdProc */ 00758 00759 int GetOrder(order *od, 00760 int *p) 00761 { 00762 int ierr=0,bbufs=2,*bbuf[2]={0}, 00763 ibufs=9,*ibuf[9]={0}, 00764 n,*iwmd; 00765 xlist *xt; 00766 00767 n=od->nnod; 00768 00769 ierr=XtAlloc(n,n+1,"xt, GetOrder",&xt); if(ierr) return FALSE; 00770 ierr=iAlloc(n,"ibuf21, GetOrder",&iwmd); if(ierr) return FALSE; 00771 00772 IptAlloc(ibufs,n,ibuf,"ibuf, GetOrder"); 00773 IptAlloc(bbufs,n,bbuf,"bbuf, GetOrder"); 00774 00775 OdProc(od,xt,ibuf[0],ibuf[1],ibuf[2],ibuf[3],ibuf[4],ibuf[5], 00776 ibuf[6],ibuf[7],ibuf[8],iwmd,bbuf[0],bbuf[1],p); 00777 00778 00779 /* 00780 XtFree(&xt); 00781 */ 00782 00783 free(xt->head); 00784 free(xt->port); 00785 free(xt->fwrd); 00786 free(xt->bwrd); 00787 free(xt); 00788 00789 iFree(&iwmd); 00790 IptFree(ibufs,ibuf); 00791 IptFree(bbufs,bbuf); 00792 return TRUE; 00793 } /* GetOrder */