1 #include "numchol.h"
2
OdAlloc(int nnod,int nn0,char * info,order ** rr)3 int OdAlloc(int nnod,
4 int nn0,
5 char *info,
6 order **rr)
7 {
8 order *r;
9 int ierr=0;
10
11 r=(order*)calloc(1,sizeof(order));
12 if (!r) ExitProc(OutOfSpc,info);
13
14 r->nnod=nnod;
15 r->nn0 =nn0;
16
17 ierr=iAlloc(nn0,info,&r->adjn); if(ierr) return 1;
18 ierr=iAlloc(nnod,info,&r->rbeg); if(ierr) return 1;
19 ierr=iAlloc(nnod,info,&r->rexs); if(ierr) return 1;
20 ierr=iAlloc(nnod,info,&r->rlen); if(ierr) return 1;
21 ierr=iAlloc(nnod,info,&r->rend); if(ierr) return 1;
22 ierr=iAlloc(nnod,info,&r->pres); if(ierr) return 1;
23 ierr=iAlloc(nnod,info,&r->succ); if(ierr) return 1;
24 *rr=r;
25 return (0);
26 } /* OdAlloc */
27
OdFree(order ** od)28 void OdFree(order **od)
29 {
30 order *r;
31
32 if (*od) {
33 r=*od;
34 iFree(&r->adjn);
35 iFree(&r->rbeg);
36 iFree(&r->rexs);
37 iFree(&r->rlen);
38 iFree(&r->rend);
39 iFree(&r->pres);
40 iFree(&r->succ);
41 free(*od);
42 *od=NULL;
43 }
44 } /* OdFree */
45
OdInit(order * od,int * nnzi)46 void OdInit(order *od,
47 int *nnzi)
48 {
49 int i,n=od->nnod;
50
51 if (n) {
52 od->rexs[0]=nnzi[0];
53 od->rlen[0]=nnzi[0];
54 od->rbeg[0]=0;
55 od->pres[0]=n;
56 od->succ[0]=1;
57 for(i=1; i<od->nnod; ++i) {
58 od->pres[i]=i-1;
59 od->succ[i]=i+1;
60 od->rexs[i]=nnzi[i];
61 od->rlen[i]=nnzi[i];
62 od->rbeg[i]=od->rbeg[i-1]+od->rlen[i-1];
63 }
64
65 od->succ[n-1]=n;
66 od->last =n-1;
67
68 od->raft=od->rbeg[n-1]+nnzi[n-1];
69
70 if (od->raft>od->nn0)
71 ExitProc(OutOfSpc,"InitMmd");
72 }
73 } /* OdInit */
74
OdSet(order * od,int allow_eli,xlist * elist,int * node_status,int * marker,int * isize,int * ilink,int * oinfo,int * osize,int * e,int * p)75 static void OdSet(order *od,
76 int allow_eli,
77 xlist *elist,
78 int *node_status,
79 int *marker,
80 int *isize,
81 int *ilink,
82 int *oinfo,
83 int *osize,
84 int *e,
85 int *p)
86 {
87 int i,n,deg,*rbeg,*rexs,*rend;
88
89 n =od->nnod;
90 rbeg=od->rbeg;
91 rexs=od->rexs;
92 rend=od->rend;
93 *e =0;
94
95 for (i=0; i<n; i++) {
96 isize[i] =0;
97 ilink[i] =n;
98 osize[i] =0;
99 oinfo[i] =n;
100 marker[i]=0;
101 }
102
103 for(i=0; i<n; ++i) {
104 rbeg[i] -= rexs[i];
105 rend[i] = 0;
106 }
107
108 for(i=0; i<n; ++i) {
109 deg = rexs[i];
110 if (!allow_eli||deg) {
111 node_status[i]=1;
112 XtPut(elist,i,deg);
113 }
114 else {
115 node_status[i] = 0;
116 marker[i] = TRUE;
117 p[*e] = i;
118 (*e)++;
119 }
120 }
121 } /* OdSet */
122
OdIndex(order * od,int i,int j)123 void OdIndex(order *od,
124 int i,
125 int j)
126 {
127 if (i!=j) {
128 od->adjn[od->rbeg[i]++]=j;
129 od->adjn[od->rbeg[j]++]=i;
130 }
131 } /* OdIndex */
132
OdArriv(order * od,int * node_status,int * marker,int * isize,int x,int * xdeg,int * rsze,int * esze,int * rchset)133 static void OdArriv(order *od,
134 int *node_status,
135 int *marker,
136 int *isize,
137 int x,
138 int *xdeg,
139 int *rsze,
140 int *esze,
141 int *rchset)
142 {
143 int *visited,i,n,y,z,l,s,t,f,stopt,
144 stops,*adjn,*rbeg,*rexs,*rend;
145
146 n =od->nnod;
147 adjn =od->adjn;
148 rbeg =od->rbeg;
149 rexs =od->rexs;
150 rend =od->rend;
151 *rsze=0;
152 *esze=0;
153
154 if (rexs[x]) {
155 l=n;
156
157 visited=marker;
158 visited[x]=TRUE;
159
160 for(t=rbeg[x], stopt=rbeg[x]+rend[x]; t<stopt; ++t) {
161 y=adjn[t];
162
163 if (node_status[y]!=0) {
164 l--;
165 rchset[l]=y;
166 visited[y]=TRUE;
167
168 for(s=rbeg[y], stops=rbeg[y]+rexs[y]; s<stops; ++s) {
169 z=adjn[s];
170
171 if (node_status[z]!=0) {
172 if (!visited[z]) {
173 visited[z]=TRUE;
174
175 rchset[*rsze]=z;
176 (*rsze)++;
177 }
178 }
179 }
180 }
181 }
182
183 f=rbeg[x]+rend[x];
184 for(t=f, stopt=rbeg[x]+rexs[x]; t<stopt; ++t) {
185 y=adjn[t];
186 if (!visited[y]) {
187 adjn[f++]=y;
188 visited[y]=TRUE;
189
190 rchset[*rsze]=y;
191 (*rsze)++;
192 }
193 }
194
195 rexs[x]=f-rbeg[x];
196
197 *esze=n-l;
198 visited[x] = FALSE;
199 iZero(*rsze,visited,rchset);
200 iZero(n-l,visited,rchset+l);
201 }
202
203 if (xdeg) {
204 *xdeg = *rsze+isize[x];
205 for(i=0; i<*rsze; ++i)
206 *xdeg+=isize[rchset[i]];
207 }
208 } /* OdArriv */
209
OdRenew(order * od,int * ilink,int x,int xdeg,int * e,int * p)210 static void OdRenew(order *od,
211 int *ilink,
212 int x,
213 int xdeg,
214 int *e,
215 int *p)
216 {
217 int c,n;
218
219 n =od->nnod;
220 od->ntot+=xdeg--;
221 p[*e] =x;
222 (*e)++;
223 for(c=x; ilink[c]!=n; c=ilink[c]) {
224 od->ntot+=xdeg--;
225 p[*e] =ilink[c];
226 (*e)++;
227 }
228 } /* OdRenew */
229
OdCheck(order * od,int * node_status)230 static void OdCheck(order *od,
231 int *node_status)
232 {
233 int f,i,t,stopt,rnew,z,previous,n,*adjn,
234 *rbeg,*rexs,*rlen,*rend,*pres,*succ;
235
236 n =od->nnod;
237 adjn=od->adjn;
238 rbeg=od->rbeg;
239 rexs=od->rexs;
240 rlen=od->rlen;
241 rend=od->rend;
242 pres=od->pres;
243 succ=od->succ;
244
245 f=0;
246 previous=n;
247 for(i=od->head; i!=n; i=succ[i]) {
248 if (node_status[i]!=0) {
249 rnew=f;
250
251 for(t=rbeg[i], stopt=rbeg[i]+rend[i]; t<stopt; ++t) {
252 z=adjn[t];
253 if (node_status[z]==3)
254 adjn[f++]=z;
255 }
256
257 rend[i]=f-rnew;
258
259 for(stopt=rbeg[i]+rexs[i]; t<stopt; ++t) {
260 z=adjn[t];
261 if (node_status[z]!=0)
262 adjn[f++]=z;
263 }
264
265 rexs[i]=f-rnew;
266 rlen[i]=rexs[i];
267
268 rbeg[i]=rnew;
269
270 if (previous==n) {
271 od->head=i;
272 pres[i]=n;
273 }
274
275 else {
276 succ[previous]=i;
277 pres[i]=previous;
278 }
279 previous=i;
280 }
281 }
282
283 if (previous!=n) {
284 succ[previous]=n;
285 od->raft=rbeg[previous]+rexs[previous];
286 }
287
288 od->last=previous;
289 } /* OdCheck */
290
OdAdd(order * od,int * node_status,int x,int newsze)291 static void OdAdd(order *od,
292 int *node_status,
293 int x,
294 int newsze)
295 {
296 int n,*adjn,*rbeg,*rexs,*rlen,*pres,*succ;
297
298 n =od->nnod;
299 adjn=od->adjn;
300 rbeg=od->rbeg;
301 rexs=od->rexs;
302 rlen=od->rlen;
303 pres=od->pres;
304 succ=od->succ;
305
306 if (newsze<=rlen[x])
307 return;
308
309 if (od->raft+newsze>od->nn0)
310 OdCheck(od,node_status);
311
312 if (od->raft+newsze>od->nn0)
313 ExitProc(OutOfSpc,"OdAdd");
314
315 if (pres[x]!=n)
316 rlen[pres[x]]+=rlen[x];
317
318 iCopy(rexs[x],adjn+rbeg[x],adjn+od->raft);
319 rbeg[x]=od->raft;
320 rlen[x]=newsze;
321 od->raft+=newsze;
322
323 if (pres[x]==n) {
324 if (succ[x]==n)
325 od->head=x;
326 else
327 od->head=succ[x];
328 }
329
330 else {
331 if (succ[x]==n)
332 succ[pres[x]]=x;
333 else
334 succ[pres[x]]=succ[x];
335 }
336
337 if (succ[x]!=n)
338 pres[succ[x]]=pres[x];
339
340 if (od->last!=x) {
341 succ[od->last]=x;
342 pres[x]=od->last;
343 }
344
345 succ[x] =n;
346 od->last=x;
347 } /* OdAdd */
348
OdComb(order * od,int * node_status,int * marker,int * isize,int * ilink,int * osize,int xsize,int * xset)349 static int OdComb(order *od,
350 int *node_status,
351 int *marker,
352 int *isize,
353 int *ilink,
354 int *osize,
355 int xsize,
356 int *xset)
357 {
358 int i,n,rnew,rlen,x,icur;
359
360 n =od->nnod;
361 rlen=0;
362
363 if (xsize==0)
364 rnew=n;
365 else if (xsize==1)
366 rnew=xset[0];
367 else {
368 rnew=xset[0];
369 for(i=1; i<xsize; ++i)
370 rlen+=1+isize[xset[i]];
371
372 node_status[rnew]=1;
373 osize[rnew]=0;
374
375 for(icur=rnew; ilink[icur]!=n; icur=ilink[icur]);
376 isize[rnew]+=rlen;
377
378 for(i=1; i<xsize; ++i) {
379 x=xset[i];
380
381 node_status[x]=0;
382 marker[x]=TRUE;
383
384 ilink[icur]=x;
385
386 for(icur=x; ilink[icur]!=n; icur=ilink[icur]);
387
388 isize[x]=0;
389 }
390 }
391
392 return (rnew);
393 } /* OdComb */
394
OdSelect(order * od,xlist * elist,int * node_status,int * marker,int * isize,int * ilink,int * oinfo,int * osize,int x,int * rsze,int * rchset,int * ibuf1,int * ibuf2,int * mask2,int * e,int * p)395 static int OdSelect(order *od,
396 xlist *elist,
397 int *node_status,
398 int *marker,
399 int *isize,
400 int *ilink,
401 int *oinfo,
402 int *osize,
403 int x,
404 int *rsze,
405 int *rchset,
406 int *ibuf1,
407 int *ibuf2,
408 int *mask2,
409 int *e,
410 int *p)
411 {
412 int absorp,old,i,j,n,esze,y,z,l,f,t,stopt,s,
413 o,stops,indsze,xdeg,e0,ssze,*slist,tsze,
414 *tlist,sze,*adjn,*rbeg,*rexs,*rlen,*rend;
415
416 adjn =od->adjn;
417 rbeg =od->rbeg;
418 rexs =od->rexs;
419 rlen =od->rlen;
420 rend =od->rend;
421 n =od->nnod;
422 slist=ibuf1;
423
424 e0 = *e;
425 OdArriv(od,node_status,marker,isize,x,&xdeg,rsze,&esze,rchset);
426
427 XtDel(elist,x);
428
429 OdRenew(od,ilink,x,xdeg,e,p);
430
431 for(i=n-esze; i<n; ++i) {
432 node_status[rchset[i]]=0;
433 marker[rchset[i]]=TRUE;
434 }
435
436 marker[x]=TRUE;
437 iSet(*rsze,TRUE,marker,rchset);
438
439 ssze=0;
440 for(i=0; i<*rsze;) {
441 y=rchset[i];
442
443 if (node_status[y]==0||node_status[y]==3)
444 ExitProc(SysError,NULL);
445
446 f=rbeg[y];
447 for(t=f, stopt=f+rend[y]; t<stopt; ++t) {
448 z=adjn[t];
449 if (node_status[z]==3) {
450 adjn[f++]=z;
451
452 if (!mask2[z]) {
453 slist[ssze++]=z;
454 mask2[z]=TRUE;
455 }
456 }
457 }
458 rend[y]=f-rbeg[y];
459
460 for(stopt=rbeg[y]+rexs[y]; t<stopt; ++t) {
461 z=adjn[t];
462 if (!marker[z])
463 adjn[f++]=z;
464 }
465
466 rexs[y]=f-rbeg[y];
467
468 if (rexs[y]==0) {
469 OdRenew(od,ilink,y,xdeg-(*e-e0),e,p);
470 node_status[y] = 0;
471 marker[y] = TRUE;
472
473 (*rsze)--;
474 iSwap(i,*rsze,rchset);
475 }
476
477 else {
478 if (rexs[y]>=rlen[y])
479 ExitProc(SysError,NULL);
480
481 if (rexs[y]>rend[y])
482 adjn[rbeg[y]+rexs[y]]=adjn[rbeg[y]+rend[y]];
483
484 rexs[y]++;
485
486 adjn[rbeg[y]+rend[y]]=x;
487 rend[y]++;
488
489 i++;
490 }
491 }
492
493 iSet(ssze,FALSE,mask2,slist);
494
495 if (*rsze==0) {
496 node_status[x]=0;
497 marker[x]=TRUE;
498 }
499
500 else {
501 node_status[x]=3;
502
503 rend[x]=0;
504 rexs[x]=0;
505 if (*rsze>rlen[x])
506 OdAdd(od,node_status,x,*rsze);
507
508 rexs[x]=*rsze;
509 iCopy(*rsze,rchset,adjn+rbeg[x]);
510
511 tsze=0;
512 tlist=ibuf2;
513 for(i=0; i<ssze; ++i){
514 y=slist[i];
515 old=marker[y];
516 marker[y]=TRUE;
517
518 absorp=TRUE;
519
520 indsze=n;
521 l=n;
522
523 f=rbeg[y];
524 for(t=f, stopt=f+rexs[y]; t<stopt; ++t) {
525 z=adjn[t];
526 if (node_status[z]!=0) {
527 adjn[f++]=z;
528
529 if (marker[z]) {
530 l--;
531 slist[l]=z;
532
533 if (!mask2[z]) {
534 for(s=rbeg[z],stops=rbeg[z]+rexs[z];
535 s<stops &&marker[adjn[s]]; ++s);
536
537 if (s==stops) {
538 indsze--;
539 iSwap(l,indsze,slist);
540 }
541
542 mask2[z]=TRUE;
543 tlist[tsze++]=z;
544 }
545 }
546 else
547 absorp=FALSE;
548 }
549 }
550
551 marker[y]=old;
552 rexs[y]=f-rbeg[y];
553
554 if (indsze<n) {
555 z=OdComb(od,node_status,marker,
556 isize,ilink,osize,
557 n-indsze,slist+indsze);
558
559 node_status[z]=1;
560
561 sze=0;
562 for(j=l; j<indsze; ++j) {
563 o=slist[j];
564 sze+=1+isize[o];
565 node_status[o]=2;
566 oinfo[o]=z;
567 }
568 osize[z]=max(osize[z],sze);
569 }
570
571 if (absorp) {
572 node_status[y]=0;
573 marker[y]=TRUE;
574 }
575 }
576
577 iSet(tsze,FALSE,mask2,tlist);
578 }
579
580 marker[x]=(node_status[x]==0);
581
582 for(t=0; t<*rsze; ++t) {
583 z=rchset[t];
584 marker[z]=(node_status[z]==0);
585 }
586
587 return (FALSE);
588 } /* OdSelect */
589
OdOrder(order * od,int * node_status,int * marker,int * isize,int x,int * ibuf1)590 static int OdOrder(order *od,
591 int *node_status,
592 int *marker,
593 int *isize,
594 int x,
595 int *ibuf1)
596 {
597 int rsze,esze,deg;
598
599 OdArriv(od,node_status,marker,isize,
600 x,°,&rsze,&esze,ibuf1);
601
602 return deg;
603 } /* OdOrder */
604
OdModf(order * od,xlist * elist,int * node_status,int * marker,int * isize,int * oinfo,int rsze,int * rchset,int * ibuf1)605 static void OdModf(order *od,
606 xlist *elist,
607 int *node_status,
608 int *marker,
609 int *isize,
610 int *oinfo,
611 int rsze,
612 int *rchset,
613 int *ibuf1)
614 {
615
616 int i,x,deg;
617
618 for(i=0; i<rsze; ++i) {
619 x=rchset[i];
620 if (node_status[x]==2)
621 if (node_status[oinfo[x]]==0||node_status[oinfo[x]]==3)
622 node_status[x]=1;
623
624 if (node_status[x]==1) {
625 deg=OdOrder(od,node_status,marker,isize,x,ibuf1);
626 XtPut(elist,x,deg-isize[x]);
627 }
628
629 else
630 XtDel(elist,x);
631 }
632 } /* mindeg_upddeg */
633
OdProc(order * od,xlist * xt,int * bbuf1,int * bbuf2,int * ibuf1,int * ibuf2,int * ibuf3,int * ibuf4,int * ibuf5,int * ibuf6,int * ibuf7,int * ibuf8,int * ibuf9,int * intbuf1,int * p)634 void OdProc(order *od,
635 xlist *xt,
636 int *bbuf1,
637 int *bbuf2,
638 int *ibuf1,
639 int *ibuf2,
640 int *ibuf3,
641 int *ibuf4,
642 int *ibuf5,
643 int *ibuf6,
644 int *ibuf7,
645 int *ibuf8,
646 int *ibuf9,
647 int *intbuf1,
648 int *p)
649 {
650 int *mask2,total_fillin,use_mtmd=TRUE,*marker=bbuf1,
651 *node_status,i,n,e,x,y,rsze,deg,mindeg,xsize,sze,
652 j,*isize,*ilink,*oinfo,*osize,*rchset,*dependent,
653 *slist;
654 xlist *elist;
655
656 elist=xt;
657
658 isize =ibuf1;
659 ilink =ibuf2;
660 oinfo =ibuf3;
661 osize =ibuf4;
662 rchset =ibuf5;
663 dependent=ibuf6;
664 slist =ibuf7;
665
666 node_status=intbuf1;
667 mask2=bbuf2;
668
669 OdSet(od,TRUE,elist,node_status,marker,
670 isize,ilink,oinfo,osize,&e,p);
671
672 n=od->nnod;
673
674 iSet(n,0,dependent,NULL);
675
676 total_fillin=FALSE;
677 for(; e<n && !total_fillin;) {
678
679 XtLeast(elist);
680
681 if (!XtGet(elist,&y,&mindeg)) {
682 printf("\n No new nodes e=%d n=%d",e,n);
683 printf(" Node status: ");
684
685 for(i=0; i<n; ++i)
686 if (node_status[i]==1)
687 printf("A\n");
688 else if (node_status[i]==2)
689 printf("\n O%d: rlen=%d oinfo=%d\n",
690 i,isize[i],oinfo[i]);
691
692 ExitProc(SysError,NULL);
693 }
694
695 xsize=0;
696 for(;use_mtmd;) {
697 if (!XtGet(elist,&x,°)||deg>mindeg)
698 break;
699
700 if (node_status[x]!=1)
701 XtDel(elist,x);
702
703 else {
704 XtSucc(elist);
705 if (!dependent[x]) {
706
707 total_fillin = OdSelect(od,elist,node_status,marker,
708 isize,ilink,oinfo,osize,x,
709 &rsze,rchset,ibuf8,ibuf9,mask2,
710 &e,p);
711
712 if (!total_fillin) {
713 dependent[x]=2;
714 slist[xsize++]=x;
715
716 for(i=0; i<rsze; ++i) {
717 y=rchset[i];
718 if (!dependent[y]) {
719 dependent[y]=1;
720 slist[xsize++]=y;
721 }
722 }
723 }
724
725 if (!use_mtmd) break;
726 }
727 }
728 }
729
730 if (!total_fillin) {
731 sze=0;
732 for(j=0; j<xsize; ++j) {
733 y=slist[j];
734 if (dependent[y]==1 && node_status[y]!=0)
735 slist[sze++]=y;
736 dependent[y]=0;
737 }
738
739 OdModf(od,elist,node_status,marker,
740 isize,oinfo,sze,slist,ibuf8);
741
742 }
743 }
744
745 if (e<n) {
746 sze=0;
747 for(i=0; i<n; ++i)
748 if (node_status[i]==2||node_status[i]==1)
749 ibuf8[sze++]=i;
750
751 x = OdComb(od,node_status,marker,isize,ilink,osize,
752 sze,ibuf8);
753
754 OdRenew(od,ilink,x,n-e-1,&e,p);
755 }
756
757 } /* OdProc */
758
GetOrder(order * od,int * p)759 int GetOrder(order *od,
760 int *p)
761 {
762 int ierr=0,bbufs=2,*bbuf[2]={0},
763 ibufs=9,*ibuf[9]={0},
764 n,*iwmd;
765 xlist *xt;
766
767 n=od->nnod;
768
769 ierr=XtAlloc(n,n+1,"xt, GetOrder",&xt); if(ierr) return FALSE;
770 ierr=iAlloc(n,"ibuf21, GetOrder",&iwmd); if(ierr) return FALSE;
771
772 IptAlloc(ibufs,n,ibuf,"ibuf, GetOrder");
773 IptAlloc(bbufs,n,bbuf,"bbuf, GetOrder");
774
775 OdProc(od,xt,ibuf[0],ibuf[1],ibuf[2],ibuf[3],ibuf[4],ibuf[5],
776 ibuf[6],ibuf[7],ibuf[8],iwmd,bbuf[0],bbuf[1],p);
777
778
779 /*
780 XtFree(&xt);
781 */
782
783 free(xt->head);
784 free(xt->port);
785 free(xt->fwrd);
786 free(xt->bwrd);
787 free(xt);
788
789 iFree(&iwmd);
790 IptFree(ibufs,ibuf);
791 IptFree(bbufs,bbuf);
792 return TRUE;
793 } /* GetOrder */
794