12 ipt[i]=(
int*)calloc(n,
sizeof(
int));
14 ExitProc(OutOfSpc,info);
37 for(j=0; j<n && i!=v[j]; ++j);
41static void InsertDplRow(
int i,
53void PermTransSym(
int nrow,
65 iZero(nrow,nnzt,NULL);
69 for(s=0; s<nrow; ++s) {
71 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
78 for(j=0; j<nrow; ++j) {
79 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
89 for(s=0; s<nrow; ++s) {
91 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
98 for(j=0; j<nrow; ++j) {
99 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
108 for(i=1; i<nrow; ++i) {
109 firt[i] =firt[i-1] + nnzt[i-1];
116 for(s=0; s<nrow; ++s) {
118 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
121 subt[firt[i]+nnzt[i]]=j;
125 subt[firt[j]+nnzt[j]]=i;
132 for(j=0; j<nrow; ++j) {
133 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
136 subt[firt[i]+nnzt[i]]=j;
140 subt[firt[j]+nnzt[j]]=i;
150 for(s=0; s<nrow; ++s) {
152 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
155 subt[firt[i]+nnzt[i]]=j;
159 subt[firt[j]+nnzt[j]]=i;
166 for(j=0; j<nrow; ++j) {
167 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
170 subt[firt[i]+nnzt[i]]=j;
174 subt[firt[j]+nnzt[j]]=i;
183static void LocDplRow(
int nrow,
195 int i,rnew,n,oisze,isze,s,count,
196 temp,k,nexti,*cur=i1nrow;
201 for (i=0; i<nrow; i++) {
206 iSet(ncol,n,ifirst,NULL);
212 for(i=0; i<nrow; ++i) {
214 for(; cur[i]<sze[i] && !map[sub[fir[i]+cur[i]]]; ++cur[i]);
216 if ( cur[i]<sze[i] ) {
217 s=sub[fir[i]+cur[i]];
221 InsertDplRow(i,s,ifirst,ilink);
233 for(k=oisze; k<isze; ++k) {
234 temp =ifirst[ilist[k]];
241 ilist[nrow-count]=rnew;
258 for(; i!=n; i=nexti) {
263 for(; cur[i]<sze[i] && !map[sub[fir[i]+cur[i]]]; ++cur[i]);
266 s =sub[fir[i]+cur[i]];
285 for(k=oisze; k<isze; ++k) {
286 temp =ifirst[ilist[k]];
293 ilist[nrow-count]=rnew;
298 for(k=0; k<count; ++k)
299 ilist[k]=ilist[nrow-count+k];
302static int CompIntElem(
const void *e1,
317static void iSort(
int n,
320 qsort((
void *)x,n,
sizeof(
int),CompIntElem);
323static void DetectDenseNodes(chfac *sf,
331 int ierr=0,j,k,l,t,ndens,nil=sf->nrow,
332 *subg=sf->subg,*ujbeg=sf->ujbeg,
333 *ujsze=sf->ujsze,*usub=sf->usub,
334 *fir,*sze,*ilist,*ilink;
337 !i1nrow || !i2nrow || !i3nrow ||
338 !i4nrow || !i5nrow || !i6nrow) {
352 for(k=0; k+1<sf->nsnds; ++k) {
354 for(t=0; t<ujsze[j] && usub[ujbeg[j]+t]<l; ++t);
360 LocDplRow(sf->nsnds-1,sf->nrow,fir,sze,usub,
362 i6nrow,ilink,ilist,&ndens,i5nrow);
364 ierr=iAlloc(ndens+1,
"sf->dhead, DetectDenseNodes",&sf->dhead);
if(ierr)
return;
365 ierr=iAlloc(sf->nsnds,
"sf->dsub, DetectDenseNodes",&sf->dsub);
if(ierr)
return;
366 ierr=iAlloc(sf->nsnds,
"sf->dbeg, DetectDenseNodes",&sf->dbeg);
if(ierr)
return;
373 for(k=0; k<ndens; ++k) {
376 sf->dhead[sf->ndens+1]=sf->dhead[sf->ndens];
378 for(; j!=nil; j=ilink[j]) {
379 sf->nsndn += sf->subg[j+1]-sf->subg[j];
380 sf->dsub[sf->dhead[sf->ndens]]=j;
381 sf->dbeg[sf->dhead[sf->ndens]]=fir[j]-ujbeg[subg[j]];
382 sf->dhead[sf->ndens]++;
384 iSort(sf->dhead[sf->ndens]-sf->dhead[sf->ndens-1],
385 sf->dsub+sf->dhead[sf->ndens-1]);
387 for(t=sf->dhead[sf->ndens-1]; t<sf->dhead[sf->ndens]; ++t)
388 sf->dbeg[t]=fir[sf->dsub[t]]-ujbeg[subg[sf->dsub[t]]];
393static int ChlSymb(chfac *sf,
396 int ierr=0,chksn,i,j,t,stopt,sze,first,cur,k,
397 ffree=0,ipos,nrow=sf->nrow,nil=nrow,
398 *nnz,*fir,*pssub,*link,*buf,*mask,
399 *usub,*tusub,*i1nrow,*i2nrow,*i3nrow,
400 *i4nrow,*p=sf->perm,*invp=sf->invp,
401 *ujbeg=sf->ujbeg,*ujsze=sf->ujsze,
404 ierr=iAlloc(sf->snnz,
"pssub, ChlSymb",&pssub);
if(ierr)
return FALSE;
406 for(i=0; i<nrow; ++i)
412 PermTransSym(nrow,sf->shead,sf->ssize,sf->ssub,
413 invp,TRUE,fir,nnz,pssub);
415 PermTransSym(nrow,fir,nnz,pssub,NULL,FALSE,
416 sf->shead,sf->ssize,sf->ssub);
421 ierr =iAlloc(k,
"usub, ChlSymb",&usub);
if(ierr)
return FALSE;
428 iZero(nrow,mask,NULL);
429 iSet(nrow,nil,link,NULL);
434 for(i=0; i<nrow; ++i) {
442 subg[sf->nsnds+1] =1 + subg[sf->nsnds];
447 iCopy(sze,sf->ssub+sf->shead[i],usub+ujbeg[i]);
449 first=usub[ujbeg[i]];
450 for(cur=first; link[cur]!=nil; cur=link[cur]);
451 link[cur] =sf->nsnds;
460 iCopy(sze,sf->ssub+sf->shead[i],buf);
461 iSet(sze,1,mask,buf);
463 for(; cur!=nil; cur=link[cur]) {
464 chksn |= (1+cur==sf->nsnds);
467 for(t=ujbeg[k], stopt=t+ujsze[k]; t<stopt; ++t) {
469 if ( j>i && !mask[j] ) {
478 k =subg[sf->nsnds-1];
479 chksn=sze==( ujsze[k]-(subg[sf->nsnds]-subg[sf->nsnds-1]) );
484 for(t=0; t<sze; ++t) {
491 ipos=LocIntPos(ujsze[i-1],i,usub+ujbeg[i-1]);
493 if (ipos==ujsze[i-1])
494 ExitProc(SysError,NULL);
496 iSwap(ujbeg[i-1],ipos+ujbeg[i-1],usub);
499 ujbeg[i] =ujbeg[i-1]+1;
500 ujsze[i] =ujsze[i-1]-1;
502 if (usub[ujbeg[i]-1]!=i)
503 ExitProc(SysError,NULL);
506 for(cur=first; link[cur]!=nil; cur=link[cur]);
507 link[cur] =sf->nsnds-1;
508 link[sf->nsnds-1]=nil;
513 subg[sf->nsnds+1] =1 + subg[sf->nsnds];
519 ExitProc(SysError,NULL);
521 iCopy(sze,buf,usub+ujbeg[i]);
524 for(cur=first; link[cur]!=nil; cur=link[cur]);
525 link[cur] =sf->nsnds;
532 if (ujsze[i]+1==nrow-i)
536 for(++i; i<nrow; ++i) {
537 ujsze[i] =ujsze[i-1]-1;
538 ujbeg[i]=ujbeg[i-1]+1;
543 ierr=iAlloc(ffree,
"tusub, ChlSymb",&tusub);
if(ierr)
return FALSE;
548 iZero(nrow,nnz,NULL);
550 for(k=0; k<sf->nsnds; ++k) {
552 plusXs(ujsze[j],nnz,usub+ujbeg[j]);
556 for(k=1; k<nrow; ++k)
557 fir[k]=fir[k-1] + nnz[k-1];
559 iZero(nrow,nnz,NULL);
561 for(k=0; k<sf->nsnds; ++k) {
563 for(t=ujbeg[j], stopt=t+ujsze[j]; t<stopt; ++t) {
565 tusub[fir[i]+nnz[i]] =j;
571 for(i=0; i<nrow; ++i) {
572 for(t=fir[i], stopt=t+nnz[i]; t<stopt; ++t) {
574 usub[ujbeg[j]+ujsze[j]] =i;
581 if (ffree<=sf->ujnz) {
582 iCopy(ffree,usub,sf->usub);
590 ierr=iAlloc(ffree,
"sf->usub, ChlSymb",&sf->usub);
if(ierr)
return FALSE;
591 iCopy(ffree,usub,sf->usub);
597 ierr=iAlloc(4*nrow,
"i1nrow, ChlSymb",&i1nrow);
if(ierr)
return FALSE;
602 DetectDenseNodes(sf,sf->uhead,sf->invp,
603 i1nrow,i2nrow,i3nrow,i4nrow);
608 for(i=1; i<nrow; ++i)
609 sf->uhead[i]=sf->uhead[i-1]+sf->ujsze[i-1];
611 for(i=0; i<nrow; ++i)
614 for(k=0; k<sf->nsnds; ++k)
615 if ( subg[k]+1!=subg[k+1] )
618 ierr=iAlloc(3*sf->n,NULL,&sf->iw);
if(ierr)
return FALSE;
619 ierr=dAlloc(2*sf->n,NULL,&sf->rw);
if(ierr)
return FALSE;
625int SymbProc(
int* isze,
630 int ierr,i,k,t,snnz,lnnz,*nnzi,nrow,resp;
637 ierr=CfcAlloc(n,
"sdt->sf, SymbProc",&cf);
if(ierr)
return FALSE;
640 for (snnz=0,i=0; i<nrow; i++)
646 ierr=iAlloc(snnz,
"cf, SymbProc",&cf->ssub);
if(ierr)
return FALSE;
650 iZero(nrow,nnzi,NULL);
653 for (i=0; i<nrow; i++) {
659 iCopy(k,jsub,cf->ssub);
662 iZero(nrow,nnzi,NULL);
664 for (i=0; i<nrow; i++) {
665 nnzi[i]+=cf->ssize[i];
666 plusXs(cf->ssize[i],nnzi,cf->ssub+cf->shead[i]);
669 ierr=OdAlloc(nrow,2*cf->snnz,
"od, PspSymbo",&od);
if(ierr)
return FALSE;
673 for (i=0; i<nrow; i++)
674 for (t=0; t<cf->ssize[i]; ++t)
675 OdIndex(od,i,cf->ssub[cf->shead[i]+t]);
677 GetOrder(od,cf->perm);
682 resp=ChlSymb(cf,lnnz);
683 LvalAlloc(cf,
"cf, PspSymb");
690int MchlSetup2(
int m, chfac** A)
692 int ierr,i,j,k,lnnz,mnnz;
695 ierr =CfcAlloc(m,NULL,&mf);
if(ierr)
return 1;
703 ierr=iAlloc(mnnz,NULL,&mf->ssub);
if(ierr)
return 1;
712 for (j=0; j<mnnz; j++)
713 mf->ssub[k+j]=(j+1+i);
736 ierr=LvalAlloc(mf,
"cf, PspSymb");
if(ierr)
return 1;