00001 #include "numchol.h"
00002
00003 int IptAlloc(int m,
00004 int n,
00005 int *ipt[],
00006 char *info)
00007 {
00008 int i;
00009
00010 if (n) {
00011 for (i=0; i<m; i++) {
00012 ipt[i]=(int*)calloc(n,sizeof(int));
00013 if (!ipt[i]){
00014 ExitProc(OutOfSpc,info);
00015 return 1;
00016 }
00017 }
00018 }
00019 return 0;
00020 }
00021
00022 void IptFree(int m,
00023 int *ipt[])
00024 {
00025 int i;
00026
00027 for (i=0; i<m; i++)
00028 iFree(&ipt[i]);
00029 }
00030
00031 int LocIntPos(int n,
00032 int i,
00033 int *v)
00034 {
00035 int j;
00036
00037 for(j=0; j<n && i!=v[j]; ++j);
00038 return (j);
00039 }
00040
00041 static void InsertDplRow(int i,
00042 int ifir,
00043 int *ifirst,
00044 int *ilink)
00045 {
00046 int temp;
00047
00048 temp=ifirst[ifir];
00049 ifirst[ifir]=i;
00050 ilink[i]=temp;
00051 }
00052
00053 void PermTransSym(int nrow,
00054 int *fir,
00055 int *nnz,
00056 int *sub,
00057 int *p,
00058 int rowwise,
00059 int *firt,
00060 int *nnzt,
00061 int *subt)
00062 {
00063 int i,j,s,t,stopt;
00064
00065 iZero(nrow,nnzt,NULL);
00066
00067 if (rowwise) {
00068 if (p) {
00069 for(s=0; s<nrow; ++s) {
00070 j =p[s];
00071 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
00072 i=p[sub[t]];
00073 nnzt[max(i,j)]++;
00074 }
00075 }
00076 }
00077 else {
00078 for(j=0; j<nrow; ++j) {
00079 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
00080 i=sub[t];
00081 nnzt[max(i,j)]++;
00082 }
00083 }
00084 }
00085 }
00086
00087 else {
00088 if (p) {
00089 for(s=0; s<nrow; ++s) {
00090 j =p[s];
00091 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
00092 i=p[sub[t]];
00093 nnzt[min(i,j)]++;
00094 }
00095 }
00096 }
00097 else {
00098 for(j=0; j<nrow; ++j) {
00099 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
00100 i=sub[t];
00101 nnzt[min(i,j)]++;
00102 }
00103 }
00104 }
00105 }
00106
00107 firt[0]=0;
00108 for(i=1; i<nrow; ++i) {
00109 firt[i] =firt[i-1] + nnzt[i-1];
00110 nnzt[i-1]=0;
00111 }
00112 nnzt[nrow-1]=0;
00113
00114 if (rowwise) {
00115 if (p) {
00116 for(s=0; s<nrow; ++s) {
00117 j=p[s];
00118 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
00119 i=p[sub[t]];
00120 if (i>j) {
00121 subt[firt[i]+nnzt[i]]=j;
00122 nnzt[i]++;
00123 }
00124 else {
00125 subt[firt[j]+nnzt[j]]=i;
00126 nnzt[j]++;
00127 }
00128 }
00129 }
00130 }
00131 else {
00132 for(j=0; j<nrow; ++j) {
00133 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
00134 i=sub[t];
00135 if (i>j) {
00136 subt[firt[i]+nnzt[i]]=j;
00137 nnzt[i]++;
00138 }
00139 else {
00140 subt[firt[j]+nnzt[j]]=i;
00141 nnzt[j]++;
00142 }
00143 }
00144 }
00145 }
00146 }
00147
00148 else {
00149 if (p) {
00150 for(s=0; s<nrow; ++s) {
00151 j=p[s];
00152 for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
00153 i=p[sub[t]];
00154 if (i<j) {
00155 subt[firt[i]+nnzt[i]]=j;
00156 nnzt[i]++;
00157 }
00158 else {
00159 subt[firt[j]+nnzt[j]]=i;
00160 nnzt[j]++;
00161 }
00162 }
00163 }
00164 }
00165 else {
00166 for(j=0; j<nrow; ++j) {
00167 for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
00168 i=sub[t];
00169 if (i<j) {
00170 subt[firt[i]+nnzt[i]]=j;
00171 nnzt[i]++;
00172 }
00173 else {
00174 subt[firt[j]+nnzt[j]]=i;
00175 nnzt[j]++;
00176 }
00177 }
00178 }
00179 }
00180 }
00181 }
00182
00183 static void LocDplRow(int nrow,
00184 int ncol,
00185 int fir[],
00186 int sze[],
00187 int *sub,
00188 int map[],
00189 int ifirst[],
00190 int ilink[],
00191 int ilist[],
00192 int *icount,
00193 int i1nrow[])
00194 {
00195 int i,rnew,n,oisze,isze,s,count,
00196 temp,k,nexti,*cur=i1nrow;
00197
00198 n =nrow;
00199 *icount=0;
00200
00201 for (i=0; i<nrow; i++) {
00202 cur[i] =0;
00203 ilink[i]=n;
00204 ilist[i]=n;
00205 }
00206 iSet(ncol,n,ifirst,NULL);
00207
00208 isze =0;
00209 count=0;
00210 oisze=isze;
00211 rnew =n;
00212 for(i=0; i<nrow; ++i) {
00213 if (map)
00214 for(; cur[i]<sze[i] && !map[sub[fir[i]+cur[i]]]; ++cur[i]);
00215
00216 if ( cur[i]<sze[i] ) {
00217 s=sub[fir[i]+cur[i]];
00218 if ( ifirst[s]==n )
00219 ilist[isze++]=s;
00220
00221 InsertDplRow(i,s,ifirst,ilink);
00222
00223 cur[i]++;
00224 }
00225
00226 else {
00227 temp =rnew;
00228 rnew =i;
00229 ilink[i]=temp;
00230 }
00231 }
00232
00233 for(k=oisze; k<isze; ++k) {
00234 temp =ifirst[ilist[k]];
00235 ifirst[ilist[k]]=n;
00236 ilist[k] =temp;
00237 }
00238
00239 if (rnew!=n) {
00240 count++;
00241 ilist[nrow-count]=rnew;
00242 }
00243
00244 while(isze) {
00245 isze--;
00246 oisze =isze;
00247
00248 i =ilist[isze];
00249 ilist[isze]=n;
00250
00251 if (i==n)
00252 exit(0);
00253
00254 rnew=n;
00255 if (ilink[i]==n)
00256 rnew=i;
00257 else {
00258 for(; i!=n; i=nexti) {
00259 nexti =ilink[i];
00260 ilink[i]=n;
00261
00262 if ( map )
00263 for(; cur[i]<sze[i] && !map[sub[fir[i]+cur[i]]]; ++cur[i]);
00264
00265 if (cur[i]<sze[i]) {
00266 s =sub[fir[i]+cur[i]];
00267 cur[i]++;
00268
00269 if (ifirst[s]==n)
00270 ilist[isze++]=s;
00271
00272 temp =ifirst[s];
00273 ifirst[s]=i;
00274 ilink[i] =temp;
00275 }
00276
00277 else {
00278 temp =rnew;
00279 rnew =i;
00280 ilink[i]=temp;
00281 }
00282 }
00283 }
00284
00285 for(k=oisze; k<isze; ++k) {
00286 temp =ifirst[ilist[k]];
00287 ifirst[ilist[k]]=n;
00288 ilist[k] =temp;
00289 }
00290
00291 if (rnew!=n) {
00292 count++;
00293 ilist[nrow-count]=rnew;
00294 }
00295 }
00296
00297 *icount=count;
00298 for(k=0; k<count; ++k)
00299 ilist[k]=ilist[nrow-count+k];
00300 }
00301
00302 static int CompIntElem(const void *e1,
00303 const void *e2)
00304 {
00305 int *i1,*i2;
00306
00307 i1=(int *) e1;
00308 i2=(int *) e2;
00309
00310 if (*i1<*i2)
00311 return (-1);
00312 else if(*i1>*i2)
00313 return (1);
00314 return (0);
00315 }
00316
00317 static void iSort(int n,
00318 int *x)
00319 {
00320 qsort((void *)x,n,sizeof(int),CompIntElem);
00321 }
00322
00323 static void DetectDenseNodes(chfac *sf,
00324 int *i1nrow,
00325 int *i2nrow,
00326 int *i3nrow,
00327 int *i4nrow,
00328 int *i5nrow,
00329 int *i6nrow)
00330 {
00331 int ierr=0,j,k,l,t,ndens,nil=sf->nrow,
00332 *subg=sf->subg,*ujbeg=sf->ujbeg,
00333 *ujsze=sf->ujsze,*usub=sf->usub,
00334 *fir,*sze,*ilist,*ilink;
00335
00336 if (!sf->nsnds||
00337 !i1nrow || !i2nrow || !i3nrow ||
00338 !i4nrow || !i5nrow || !i6nrow) {
00339 sf->sdens=FALSE;
00340 return;
00341 }
00342
00343 sf->sdens =TRUE;
00344 fir =i1nrow;
00345 sze =i2nrow;
00346 ilist =i3nrow;
00347 ilink =i4nrow;
00348
00349 sf->nsndn=0;
00350
00351 l=subg[sf->nsnds-1];
00352 for(k=0; k+1<sf->nsnds; ++k) {
00353 j=subg[k];
00354 for(t=0; t<ujsze[j] && usub[ujbeg[j]+t]<l; ++t);
00355
00356 fir[k] =ujbeg[j]+t;
00357 sze[k] =ujsze[j]-t;
00358 }
00359
00360 LocDplRow(sf->nsnds-1,sf->nrow,fir,sze,usub,
00361 NULL,
00362 i6nrow,ilink,ilist,&ndens,i5nrow);
00363
00364 ierr=iAlloc(ndens+1, "sf->dhead, DetectDenseNodes",&sf->dhead);if(ierr) return;
00365 ierr=iAlloc(sf->nsnds,"sf->dsub, DetectDenseNodes",&sf->dsub);if(ierr) return;
00366 ierr=iAlloc(sf->nsnds,"sf->dbeg, DetectDenseNodes",&sf->dbeg);if(ierr) return;
00367
00368 nil =sf->nsnds-1;
00369 sf->ndens =0;
00370 sf->nsndn =0;
00371 sf->dhead[0]=0;
00372
00373 for(k=0; k<ndens; ++k) {
00374 j=ilist[k];
00375 if ( sze[j] ) {
00376 sf->dhead[sf->ndens+1]=sf->dhead[sf->ndens];
00377 sf->ndens++;
00378 for(; j!=nil; j=ilink[j]) {
00379 sf->nsndn += sf->subg[j+1]-sf->subg[j];
00380 sf->dsub[sf->dhead[sf->ndens]]=j;
00381 sf->dbeg[sf->dhead[sf->ndens]]=fir[j]-ujbeg[subg[j]];
00382 sf->dhead[sf->ndens]++;
00383 }
00384 iSort(sf->dhead[sf->ndens]-sf->dhead[sf->ndens-1],
00385 sf->dsub+sf->dhead[sf->ndens-1]);
00386
00387 for(t=sf->dhead[sf->ndens-1]; t<sf->dhead[sf->ndens]; ++t)
00388 sf->dbeg[t]=fir[sf->dsub[t]]-ujbeg[subg[sf->dsub[t]]];
00389 }
00390 }
00391 }
00392
00393 static int ChlSymb(chfac *sf,
00394 int ulnnz)
00395 {
00396 int ierr=0,chksn,i,j,t,stopt,sze,first,cur,k,
00397 ffree=0,ipos,nrow=sf->nrow,nil=nrow,
00398 *nnz,*fir,*pssub,*link,*buf,*mask,
00399 *usub,*tusub,*i1nrow,*i2nrow,*i3nrow,
00400 *i4nrow,*p=sf->perm,*invp=sf->invp,
00401 *ujbeg=sf->ujbeg,*ujsze=sf->ujsze,
00402 *subg=sf->subg;
00403
00404 ierr=iAlloc(sf->snnz,"pssub, ChlSymb",&pssub);if(ierr) return FALSE;
00405
00406 for(i=0; i<nrow; ++i)
00407 invp[p[i]]=i;
00408
00409 nnz=sf->uhead;
00410 fir=sf->subg;
00411
00412 PermTransSym(nrow,sf->shead,sf->ssize,sf->ssub,
00413 invp,TRUE,fir,nnz,pssub);
00414
00415 PermTransSym(nrow,fir,nnz,pssub,NULL,FALSE,
00416 sf->shead,sf->ssize,sf->ssub);
00417
00418 iFree(&pssub);
00419
00420 k =ulnnz+nrow;
00421 ierr =iAlloc(k,"usub, ChlSymb",&usub);if(ierr) return FALSE;
00422 buf =usub+ulnnz;
00423
00424 mask=sf->uhead;
00425
00426 link=invp;
00427
00428 iZero(nrow,mask,NULL);
00429 iSet(nrow,nil,link,NULL);
00430
00431 ffree =0;
00432 sf->nsnds=0;
00433 subg[0] =0;
00434 for(i=0; i<nrow; ++i) {
00435 sze =sf->ssize[i];
00436 first=nil;
00437 cur =link[i];
00438 chksn=FALSE;
00439
00440 if (cur==nil) {
00441
00442 subg[sf->nsnds+1] =1 + subg[sf->nsnds];
00443 ujsze[i] =sze;
00444 ujbeg[i] =ffree;
00445 ffree += sze;
00446
00447 iCopy(sze,sf->ssub+sf->shead[i],usub+ujbeg[i]);
00448 if (sze) {
00449 first=usub[ujbeg[i]];
00450 for(cur=first; link[cur]!=nil; cur=link[cur]);
00451 link[cur] =sf->nsnds;
00452 link[sf->nsnds]=nil;
00453 }
00454 sf->nsnds++;
00455 }
00456
00457 else {
00458 mask[i]=1;
00459
00460 iCopy(sze,sf->ssub+sf->shead[i],buf);
00461 iSet(sze,1,mask,buf);
00462
00463 for(; cur!=nil; cur=link[cur]) {
00464 chksn |= (1+cur==sf->nsnds);
00465 k =subg[cur];
00466
00467 for(t=ujbeg[k], stopt=t+ujsze[k]; t<stopt; ++t) {
00468 j=usub[t];
00469 if ( j>i && !mask[j] ) {
00470 buf[sze]=j;
00471 mask[j] =1;
00472 sze++;
00473 }
00474 }
00475 }
00476
00477 if (chksn) {
00478 k =subg[sf->nsnds-1];
00479 chksn=sze==( ujsze[k]-(subg[sf->nsnds]-subg[sf->nsnds-1]) );
00480 }
00481
00482 first =nrow;
00483 mask[i]=0;
00484 for(t=0; t<sze; ++t) {
00485 j =buf[t];
00486 mask[j]=0;
00487 first =min(j,first);
00488 }
00489
00490 if (chksn) {
00491 ipos=LocIntPos(ujsze[i-1],i,usub+ujbeg[i-1]);
00492
00493 if (ipos==ujsze[i-1])
00494 ExitProc(SysError,NULL);
00495
00496 iSwap(ujbeg[i-1],ipos+ujbeg[i-1],usub);
00497
00498 subg[sf->nsnds]++;
00499 ujbeg[i] =ujbeg[i-1]+1;
00500 ujsze[i] =ujsze[i-1]-1;
00501
00502 if (usub[ujbeg[i]-1]!=i)
00503 ExitProc(SysError,NULL);
00504
00505 if ( first!=nil ) {
00506 for(cur=first; link[cur]!=nil; cur=link[cur]);
00507 link[cur] =sf->nsnds-1;
00508 link[sf->nsnds-1]=nil;
00509 }
00510 }
00511
00512 else {
00513 subg[sf->nsnds+1] =1 + subg[sf->nsnds];
00514 ujbeg[i] =ffree;
00515 ujsze[i] =sze;
00516 ffree += sze;
00517
00518 if (ffree>ulnnz)
00519 ExitProc(SysError,NULL);
00520
00521 iCopy(sze,buf,usub+ujbeg[i]);
00522
00523 if ( first!=nil ) {
00524 for(cur=first; link[cur]!=nil; cur=link[cur]);
00525 link[cur] =sf->nsnds;
00526 link[sf->nsnds]=nil;
00527 }
00528 sf->nsnds++;
00529 }
00530 }
00531
00532 if (ujsze[i]+1==nrow-i)
00533 break;
00534 }
00535
00536 for(++i; i<nrow; ++i) {
00537 ujsze[i] =ujsze[i-1]-1;
00538 ujbeg[i]=ujbeg[i-1]+1;
00539
00540 subg[sf->nsnds]++;
00541 }
00542
00543 ierr=iAlloc(ffree,"tusub, ChlSymb",&tusub);if(ierr) return FALSE;
00544
00545 fir=buf;
00546 nnz=sf->uhead;
00547
00548 iZero(nrow,nnz,NULL);
00549
00550 for(k=0; k<sf->nsnds; ++k) {
00551 j=subg[k];
00552 plusXs(ujsze[j],nnz,usub+ujbeg[j]);
00553 }
00554
00555 fir[0]=0;
00556 for(k=1; k<nrow; ++k)
00557 fir[k]=fir[k-1] + nnz[k-1];
00558
00559 iZero(nrow,nnz,NULL);
00560
00561 for(k=0; k<sf->nsnds; ++k) {
00562 j=subg[k];
00563 for(t=ujbeg[j], stopt=t+ujsze[j]; t<stopt; ++t) {
00564 i =usub[t];
00565 tusub[fir[i]+nnz[i]] =j;
00566 nnz[i]++;
00567 }
00568 ujsze[j]=0;
00569 }
00570
00571 for(i=0; i<nrow; ++i) {
00572 for(t=fir[i], stopt=t+nnz[i]; t<stopt; ++t) {
00573 j =tusub[t];
00574 usub[ujbeg[j]+ujsze[j]] =i;
00575 ujsze[j]++;
00576 }
00577 }
00578
00579 iFree(&tusub);
00580
00581 if (ffree<=sf->ujnz) {
00582 iCopy(ffree,usub,sf->usub);
00583 iFree(&usub);
00584 }
00585
00586 else {
00587 sf->ujnz=0;
00588 iFree(&sf->usub);
00589
00590 ierr=iAlloc(ffree,"sf->usub, ChlSymb",&sf->usub);if(ierr) return FALSE;
00591 iCopy(ffree,usub,sf->usub);
00592
00593 sf->ujnz=ffree;
00594 iFree(&usub);
00595 }
00596
00597 ierr=iAlloc(4*nrow,"i1nrow, ChlSymb",&i1nrow);if(ierr) return FALSE;
00598 i2nrow=i1nrow+nrow;
00599 i3nrow=i2nrow+nrow;
00600 i4nrow=i3nrow+nrow;
00601
00602 DetectDenseNodes(sf,sf->uhead,sf->invp,
00603 i1nrow,i2nrow,i3nrow,i4nrow);
00604
00605 iFree(&i1nrow);
00606
00607 sf->uhead[0]=0;
00608 for(i=1; i<nrow; ++i)
00609 sf->uhead[i]=sf->uhead[i-1]+sf->ujsze[i-1];
00610
00611 for(i=0; i<nrow; ++i)
00612 invp[p[i]]=i;
00613
00614 for(k=0; k<sf->nsnds; ++k)
00615 if ( subg[k]+1!=subg[k+1] )
00616 break;
00617
00618 ierr=iAlloc(3*sf->n,NULL,&sf->iw);if(ierr) return FALSE;
00619 ierr=dAlloc(2*sf->n,NULL,&sf->rw);if(ierr) return FALSE;
00620 sf->factor=0;
00621
00622 return TRUE;
00623 }
00624
00625 int SymbProc(int* isze,
00626 int* jsub,
00627 int n,
00628 chfac** sf)
00629 {
00630 int ierr,i,k,t,snnz,lnnz,*nnzi,nrow,resp;
00631 chfac* cf;
00632 order *od;
00633
00634
00635
00636
00637 ierr=CfcAlloc(n,"sdt->sf, SymbProc",&cf);if(ierr) return FALSE;
00638
00639 nrow=cf->nrow;
00640 for (snnz=0,i=0; i<nrow; i++)
00641 snnz+=isze[i];
00642
00643
00644
00645
00646 ierr=iAlloc(snnz,"cf, SymbProc",&cf->ssub); if(ierr) return FALSE;
00647 cf->snnz=snnz;
00648
00649 nnzi=cf->perm;
00650 iZero(nrow,nnzi,NULL);
00651
00652 k=0;
00653 for (i=0; i<nrow; i++) {
00654 snnz =isze[i];
00655 cf->shead[i]=k;
00656 cf->ssize[i]=snnz;
00657 k+=snnz;
00658 }
00659 iCopy(k,jsub,cf->ssub);
00660
00661 nnzi=cf->perm;
00662 iZero(nrow,nnzi,NULL);
00663
00664 for (i=0; i<nrow; i++) {
00665 nnzi[i]+=cf->ssize[i];
00666 plusXs(cf->ssize[i],nnzi,cf->ssub+cf->shead[i]);
00667 }
00668
00669 ierr=OdAlloc(nrow,2*cf->snnz,"od, PspSymbo",&od); if(ierr) return FALSE;
00670 nnzi=cf->perm;
00671
00672 OdInit(od,nnzi);
00673 for (i=0; i<nrow; i++)
00674 for (t=0; t<cf->ssize[i]; ++t)
00675 OdIndex(od,i,cf->ssub[cf->shead[i]+t]);
00676
00677 GetOrder(od,cf->perm);
00678
00679 lnnz=od->ntot;
00680 OdFree(&od);
00681
00682 resp=ChlSymb(cf,lnnz);
00683 LvalAlloc(cf,"cf, PspSymb");
00684
00685
00686 *sf=cf;
00687 return resp;
00688 }
00689
00690 int MchlSetup2(int m, chfac** A)
00691 {
00692 int ierr,i,j,k,lnnz,mnnz;
00693 chfac *mf;
00694
00695 ierr =CfcAlloc(m,NULL,&mf); if(ierr) return 1;
00696 *A=mf;
00697
00698 mnnz=m*(m-1)/2;
00699 if (!mnnz && m>1)
00700 return TRUE;
00701
00702 lnnz=mnnz;
00703 ierr=iAlloc(mnnz,NULL,&mf->ssub);if(ierr) return 1;
00704
00705 mf->snnz=mnnz;
00706 k=0;
00707 for (i=0; i<m; i++)
00708 {
00709 mnnz =m-(i+1);
00710 mf->shead[i] =k;
00711 mf->ssize[i] =mnnz;
00712 for (j=0; j<mnnz; j++)
00713 mf->ssub[k+j]=(j+1+i);
00714 k +=mnnz;
00715 mf->perm[i]=i;
00716 }
00717
00718
00719 k=ChlSymb(mf,lnnz);
00720
00721 iFree(&mf->ssub);
00722 iFree(&mf->shead);
00723 iFree(&mf->ssize);
00724
00725
00726 mf->alldense=1;
00727 iFree(&mf->invp);
00728 mf->invp=mf->perm;
00729
00730 iFree(&mf->ujbeg);
00731 mf->ujbeg=mf->perm;
00732
00733 iFree(&mf->usub);
00734 mf->usub=mf->perm+1;
00735
00736 ierr=LvalAlloc(mf,"cf, PspSymb");if(ierr) return 1;
00737
00738 return 0;
00739
00740 }