DSDP
|
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 } /* AllocIptr */ 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 } /* IptFree */ 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 } /* LocIntPos */ 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 } /* InsertDplRow */ 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 } /* PermTransSym */ 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 } /* LocDplRow */ 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 } /* CompIntElem */ 00316 00317 static void iSort(int n, 00318 int *x) 00319 { 00320 qsort((void *)x,n,sizeof(int),CompIntElem); 00321 } /* iSort */ 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 } /* DetectDenseNodes */ 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 } /* ChlSymb */ 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 * set data for symmetric matrix to be factorized 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 if (!snnz) 00644 return TRUE; 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 /* sdt->sf=cf; */ 00685 00686 *sf=cf; 00687 return resp; 00688 } /* SymbProc */ 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 /* This part is different */ 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 } /* MchlSetup2 */