DSDP
sdpsymb.c
1 #include "numchol.h"
2 
3 int IptAlloc(int m,
4  int n,
5  int *ipt[],
6  char *info)
7 {
8  int i;
9 
10  if (n) {
11  for (i=0; i<m; i++) {
12  ipt[i]=(int*)calloc(n,sizeof(int));
13  if (!ipt[i]){
14  ExitProc(OutOfSpc,info);
15  return 1;
16  }
17  }
18  }
19  return 0;
20 } /* AllocIptr */
21 
22 void IptFree(int m,
23  int *ipt[])
24 {
25  int i;
26 
27  for (i=0; i<m; i++)
28  iFree(&ipt[i]);
29 } /* IptFree */
30 
31 int LocIntPos(int n,
32  int i,
33  int *v)
34 {
35  int j;
36 
37  for(j=0; j<n && i!=v[j]; ++j);
38  return (j);
39 } /* LocIntPos */
40 
41 static void InsertDplRow(int i,
42  int ifir,
43  int *ifirst,
44  int *ilink)
45 {
46  int temp;
47 
48  temp=ifirst[ifir];
49  ifirst[ifir]=i;
50  ilink[i]=temp;
51 } /* InsertDplRow */
52 
53 void PermTransSym(int nrow,
54  int *fir,
55  int *nnz,
56  int *sub,
57  int *p,
58  int rowwise,
59  int *firt,
60  int *nnzt,
61  int *subt)
62 {
63  int i,j,s,t,stopt;
64 
65  iZero(nrow,nnzt,NULL);
66 
67  if (rowwise) {
68  if (p) {
69  for(s=0; s<nrow; ++s) {
70  j =p[s];
71  for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
72  i=p[sub[t]];
73  nnzt[max(i,j)]++;
74  }
75  }
76  }
77  else {
78  for(j=0; j<nrow; ++j) {
79  for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
80  i=sub[t];
81  nnzt[max(i,j)]++;
82  }
83  }
84  }
85  }
86 
87  else {
88  if (p) {
89  for(s=0; s<nrow; ++s) {
90  j =p[s];
91  for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
92  i=p[sub[t]];
93  nnzt[min(i,j)]++;
94  }
95  }
96  }
97  else {
98  for(j=0; j<nrow; ++j) {
99  for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
100  i=sub[t];
101  nnzt[min(i,j)]++;
102  }
103  }
104  }
105  }
106 
107  firt[0]=0;
108  for(i=1; i<nrow; ++i) {
109  firt[i] =firt[i-1] + nnzt[i-1];
110  nnzt[i-1]=0;
111  }
112  nnzt[nrow-1]=0;
113 
114  if (rowwise) {
115  if (p) {
116  for(s=0; s<nrow; ++s) {
117  j=p[s];
118  for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
119  i=p[sub[t]];
120  if (i>j) {
121  subt[firt[i]+nnzt[i]]=j;
122  nnzt[i]++;
123  }
124  else {
125  subt[firt[j]+nnzt[j]]=i;
126  nnzt[j]++;
127  }
128  }
129  }
130  }
131  else {
132  for(j=0; j<nrow; ++j) {
133  for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
134  i=sub[t];
135  if (i>j) {
136  subt[firt[i]+nnzt[i]]=j;
137  nnzt[i]++;
138  }
139  else {
140  subt[firt[j]+nnzt[j]]=i;
141  nnzt[j]++;
142  }
143  }
144  }
145  }
146  }
147 
148  else {
149  if (p) {
150  for(s=0; s<nrow; ++s) {
151  j=p[s];
152  for(t=fir[s], stopt=t+nnz[s]; t<stopt; ++t) {
153  i=p[sub[t]];
154  if (i<j) {
155  subt[firt[i]+nnzt[i]]=j;
156  nnzt[i]++;
157  }
158  else {
159  subt[firt[j]+nnzt[j]]=i;
160  nnzt[j]++;
161  }
162  }
163  }
164  }
165  else {
166  for(j=0; j<nrow; ++j) {
167  for(t=fir[j], stopt=t+nnz[j]; t<stopt; ++t) {
168  i=sub[t];
169  if (i<j) {
170  subt[firt[i]+nnzt[i]]=j;
171  nnzt[i]++;
172  }
173  else {
174  subt[firt[j]+nnzt[j]]=i;
175  nnzt[j]++;
176  }
177  }
178  }
179  }
180  }
181 } /* PermTransSym */
182 
183 static void LocDplRow(int nrow,
184  int ncol,
185  int fir[],
186  int sze[],
187  int *sub,
188  int map[],
189  int ifirst[],
190  int ilink[],
191  int ilist[],
192  int *icount,
193  int i1nrow[])
194 {
195  int i,rnew,n,oisze,isze,s,count,
196  temp,k,nexti,*cur=i1nrow;
197 
198  n =nrow;
199  *icount=0;
200 
201  for (i=0; i<nrow; i++) {
202  cur[i] =0;
203  ilink[i]=n;
204  ilist[i]=n;
205  }
206  iSet(ncol,n,ifirst,NULL);
207 
208  isze =0;
209  count=0;
210  oisze=isze;
211  rnew =n;
212  for(i=0; i<nrow; ++i) {
213  if (map)
214  for(; cur[i]<sze[i] && !map[sub[fir[i]+cur[i]]]; ++cur[i]);
215 
216  if ( cur[i]<sze[i] ) {
217  s=sub[fir[i]+cur[i]];
218  if ( ifirst[s]==n )
219  ilist[isze++]=s;
220 
221  InsertDplRow(i,s,ifirst,ilink);
222 
223  cur[i]++;
224  }
225 
226  else {
227  temp =rnew;
228  rnew =i;
229  ilink[i]=temp;
230  }
231  }
232 
233  for(k=oisze; k<isze; ++k) {
234  temp =ifirst[ilist[k]];
235  ifirst[ilist[k]]=n;
236  ilist[k] =temp;
237  }
238 
239  if (rnew!=n) {
240  count++;
241  ilist[nrow-count]=rnew;
242  }
243 
244  while(isze) {
245  isze--;
246  oisze =isze;
247 
248  i =ilist[isze];
249  ilist[isze]=n;
250 
251  if (i==n)
252  exit(0);
253 
254  rnew=n;
255  if (ilink[i]==n)
256  rnew=i;
257  else {
258  for(; i!=n; i=nexti) {
259  nexti =ilink[i];
260  ilink[i]=n;
261 
262  if ( map )
263  for(; cur[i]<sze[i] && !map[sub[fir[i]+cur[i]]]; ++cur[i]);
264 
265  if (cur[i]<sze[i]) {
266  s =sub[fir[i]+cur[i]];
267  cur[i]++;
268 
269  if (ifirst[s]==n)
270  ilist[isze++]=s;
271 
272  temp =ifirst[s];
273  ifirst[s]=i;
274  ilink[i] =temp;
275  }
276 
277  else {
278  temp =rnew;
279  rnew =i;
280  ilink[i]=temp;
281  }
282  }
283  }
284 
285  for(k=oisze; k<isze; ++k) {
286  temp =ifirst[ilist[k]];
287  ifirst[ilist[k]]=n;
288  ilist[k] =temp;
289  }
290 
291  if (rnew!=n) {
292  count++;
293  ilist[nrow-count]=rnew;
294  }
295  }
296 
297  *icount=count;
298  for(k=0; k<count; ++k)
299  ilist[k]=ilist[nrow-count+k];
300 } /* LocDplRow */
301 
302 static int CompIntElem(const void *e1,
303  const void *e2)
304 {
305  int *i1,*i2;
306 
307  i1=(int *) e1;
308  i2=(int *) e2;
309 
310  if (*i1<*i2)
311  return (-1);
312  else if(*i1>*i2)
313  return (1);
314  return (0);
315 } /* CompIntElem */
316 
317 static void iSort(int n,
318  int *x)
319 {
320  qsort((void *)x,n,sizeof(int),CompIntElem);
321 } /* iSort */
322 
323 static void DetectDenseNodes(chfac *sf,
324  int *i1nrow,
325  int *i2nrow,
326  int *i3nrow,
327  int *i4nrow,
328  int *i5nrow,
329  int *i6nrow)
330 {
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;
335 
336  if (!sf->nsnds||
337  !i1nrow || !i2nrow || !i3nrow ||
338  !i4nrow || !i5nrow || !i6nrow) {
339  sf->sdens=FALSE;
340  return;
341  }
342 
343  sf->sdens =TRUE;
344  fir =i1nrow;
345  sze =i2nrow;
346  ilist =i3nrow;
347  ilink =i4nrow;
348 
349  sf->nsndn=0;
350 
351  l=subg[sf->nsnds-1];
352  for(k=0; k+1<sf->nsnds; ++k) {
353  j=subg[k];
354  for(t=0; t<ujsze[j] && usub[ujbeg[j]+t]<l; ++t);
355 
356  fir[k] =ujbeg[j]+t;
357  sze[k] =ujsze[j]-t;
358  }
359 
360  LocDplRow(sf->nsnds-1,sf->nrow,fir,sze,usub,
361  NULL,
362  i6nrow,ilink,ilist,&ndens,i5nrow);
363 
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;
367 
368  nil =sf->nsnds-1;
369  sf->ndens =0;
370  sf->nsndn =0;
371  sf->dhead[0]=0;
372 
373  for(k=0; k<ndens; ++k) {
374  j=ilist[k];
375  if ( sze[j] ) {
376  sf->dhead[sf->ndens+1]=sf->dhead[sf->ndens];
377  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]++;
383  }
384  iSort(sf->dhead[sf->ndens]-sf->dhead[sf->ndens-1],
385  sf->dsub+sf->dhead[sf->ndens-1]);
386 
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]]];
389  }
390  }
391 } /* DetectDenseNodes */
392 
393 static int ChlSymb(chfac *sf,
394  int ulnnz)
395 {
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,
402  *subg=sf->subg;
403 
404  ierr=iAlloc(sf->snnz,"pssub, ChlSymb",&pssub);if(ierr) return FALSE;
405 
406  for(i=0; i<nrow; ++i)
407  invp[p[i]]=i;
408 
409  nnz=sf->uhead;
410  fir=sf->subg;
411 
412  PermTransSym(nrow,sf->shead,sf->ssize,sf->ssub,
413  invp,TRUE,fir,nnz,pssub);
414 
415  PermTransSym(nrow,fir,nnz,pssub,NULL,FALSE,
416  sf->shead,sf->ssize,sf->ssub);
417 
418  iFree(&pssub);
419 
420  k =ulnnz+nrow;
421  ierr =iAlloc(k,"usub, ChlSymb",&usub);if(ierr) return FALSE;
422  buf =usub+ulnnz;
423 
424  mask=sf->uhead;
425 
426  link=invp;
427 
428  iZero(nrow,mask,NULL);
429  iSet(nrow,nil,link,NULL);
430 
431  ffree =0;
432  sf->nsnds=0;
433  subg[0] =0;
434  for(i=0; i<nrow; ++i) {
435  sze =sf->ssize[i];
436  first=nil;
437  cur =link[i];
438  chksn=FALSE;
439 
440  if (cur==nil) {
441 
442  subg[sf->nsnds+1] =1 + subg[sf->nsnds];
443  ujsze[i] =sze;
444  ujbeg[i] =ffree;
445  ffree += sze;
446 
447  iCopy(sze,sf->ssub+sf->shead[i],usub+ujbeg[i]);
448  if (sze) {
449  first=usub[ujbeg[i]];
450  for(cur=first; link[cur]!=nil; cur=link[cur]);
451  link[cur] =sf->nsnds;
452  link[sf->nsnds]=nil;
453  }
454  sf->nsnds++;
455  }
456 
457  else {
458  mask[i]=1;
459 
460  iCopy(sze,sf->ssub+sf->shead[i],buf);
461  iSet(sze,1,mask,buf);
462 
463  for(; cur!=nil; cur=link[cur]) {
464  chksn |= (1+cur==sf->nsnds);
465  k =subg[cur];
466 
467  for(t=ujbeg[k], stopt=t+ujsze[k]; t<stopt; ++t) {
468  j=usub[t];
469  if ( j>i && !mask[j] ) {
470  buf[sze]=j;
471  mask[j] =1;
472  sze++;
473  }
474  }
475  }
476 
477  if (chksn) {
478  k =subg[sf->nsnds-1];
479  chksn=sze==( ujsze[k]-(subg[sf->nsnds]-subg[sf->nsnds-1]) );
480  }
481 
482  first =nrow;
483  mask[i]=0;
484  for(t=0; t<sze; ++t) {
485  j =buf[t];
486  mask[j]=0;
487  first =min(j,first);
488  }
489 
490  if (chksn) {
491  ipos=LocIntPos(ujsze[i-1],i,usub+ujbeg[i-1]);
492 
493  if (ipos==ujsze[i-1])
494  ExitProc(SysError,NULL);
495 
496  iSwap(ujbeg[i-1],ipos+ujbeg[i-1],usub);
497 
498  subg[sf->nsnds]++;
499  ujbeg[i] =ujbeg[i-1]+1;
500  ujsze[i] =ujsze[i-1]-1;
501 
502  if (usub[ujbeg[i]-1]!=i)
503  ExitProc(SysError,NULL);
504 
505  if ( first!=nil ) {
506  for(cur=first; link[cur]!=nil; cur=link[cur]);
507  link[cur] =sf->nsnds-1;
508  link[sf->nsnds-1]=nil;
509  }
510  }
511 
512  else {
513  subg[sf->nsnds+1] =1 + subg[sf->nsnds];
514  ujbeg[i] =ffree;
515  ujsze[i] =sze;
516  ffree += sze;
517 
518  if (ffree>ulnnz)
519  ExitProc(SysError,NULL);
520 
521  iCopy(sze,buf,usub+ujbeg[i]);
522 
523  if ( first!=nil ) {
524  for(cur=first; link[cur]!=nil; cur=link[cur]);
525  link[cur] =sf->nsnds;
526  link[sf->nsnds]=nil;
527  }
528  sf->nsnds++;
529  }
530  }
531 
532  if (ujsze[i]+1==nrow-i)
533  break;
534  }
535 
536  for(++i; i<nrow; ++i) {
537  ujsze[i] =ujsze[i-1]-1;
538  ujbeg[i]=ujbeg[i-1]+1;
539 
540  subg[sf->nsnds]++;
541  }
542 
543  ierr=iAlloc(ffree,"tusub, ChlSymb",&tusub);if(ierr) return FALSE;
544 
545  fir=buf;
546  nnz=sf->uhead;
547 
548  iZero(nrow,nnz,NULL);
549 
550  for(k=0; k<sf->nsnds; ++k) {
551  j=subg[k];
552  plusXs(ujsze[j],nnz,usub+ujbeg[j]);
553  }
554 
555  fir[0]=0;
556  for(k=1; k<nrow; ++k)
557  fir[k]=fir[k-1] + nnz[k-1];
558 
559  iZero(nrow,nnz,NULL);
560 
561  for(k=0; k<sf->nsnds; ++k) {
562  j=subg[k];
563  for(t=ujbeg[j], stopt=t+ujsze[j]; t<stopt; ++t) {
564  i =usub[t];
565  tusub[fir[i]+nnz[i]] =j;
566  nnz[i]++;
567  }
568  ujsze[j]=0;
569  }
570 
571  for(i=0; i<nrow; ++i) {
572  for(t=fir[i], stopt=t+nnz[i]; t<stopt; ++t) {
573  j =tusub[t];
574  usub[ujbeg[j]+ujsze[j]] =i;
575  ujsze[j]++;
576  }
577  }
578 
579  iFree(&tusub);
580 
581  if (ffree<=sf->ujnz) {
582  iCopy(ffree,usub,sf->usub);
583  iFree(&usub);
584  }
585 
586  else {
587  sf->ujnz=0;
588  iFree(&sf->usub);
589 
590  ierr=iAlloc(ffree,"sf->usub, ChlSymb",&sf->usub);if(ierr) return FALSE;
591  iCopy(ffree,usub,sf->usub);
592 
593  sf->ujnz=ffree;
594  iFree(&usub);
595  }
596 
597  ierr=iAlloc(4*nrow,"i1nrow, ChlSymb",&i1nrow);if(ierr) return FALSE;
598  i2nrow=i1nrow+nrow;
599  i3nrow=i2nrow+nrow;
600  i4nrow=i3nrow+nrow;
601 
602  DetectDenseNodes(sf,sf->uhead,sf->invp,
603  i1nrow,i2nrow,i3nrow,i4nrow);
604 
605  iFree(&i1nrow);
606 
607  sf->uhead[0]=0;
608  for(i=1; i<nrow; ++i)
609  sf->uhead[i]=sf->uhead[i-1]+sf->ujsze[i-1];
610 
611  for(i=0; i<nrow; ++i)
612  invp[p[i]]=i;
613 
614  for(k=0; k<sf->nsnds; ++k)
615  if ( subg[k]+1!=subg[k+1] )
616  break;
617 
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;
620  sf->factor=0;
621 
622  return TRUE;
623 } /* ChlSymb */
624 
625 int SymbProc(int* isze,
626  int* jsub,
627  int n,
628  chfac** sf)
629 {
630  int ierr,i,k,t,snnz,lnnz,*nnzi,nrow,resp;
631  chfac* cf;
632  order *od;
633 
634  /*
635  * set data for symmetric matrix to be factorized
636  */
637  ierr=CfcAlloc(n,"sdt->sf, SymbProc",&cf);if(ierr) return FALSE;
638 
639  nrow=cf->nrow;
640  for (snnz=0,i=0; i<nrow; i++)
641  snnz+=isze[i];
642  /*
643  if (!snnz)
644  return TRUE;
645  */
646  ierr=iAlloc(snnz,"cf, SymbProc",&cf->ssub); if(ierr) return FALSE;
647  cf->snnz=snnz;
648 
649  nnzi=cf->perm;
650  iZero(nrow,nnzi,NULL);
651 
652  k=0;
653  for (i=0; i<nrow; i++) {
654  snnz =isze[i];
655  cf->shead[i]=k;
656  cf->ssize[i]=snnz;
657  k+=snnz;
658  }
659  iCopy(k,jsub,cf->ssub);
660 
661  nnzi=cf->perm;
662  iZero(nrow,nnzi,NULL);
663 
664  for (i=0; i<nrow; i++) {
665  nnzi[i]+=cf->ssize[i];
666  plusXs(cf->ssize[i],nnzi,cf->ssub+cf->shead[i]);
667  }
668 
669  ierr=OdAlloc(nrow,2*cf->snnz,"od, PspSymbo",&od); if(ierr) return FALSE;
670  nnzi=cf->perm;
671 
672  OdInit(od,nnzi);
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]);
676 
677  GetOrder(od,cf->perm);
678 
679  lnnz=od->ntot;
680  OdFree(&od);
681 
682  resp=ChlSymb(cf,lnnz);
683  LvalAlloc(cf,"cf, PspSymb");
684  /* sdt->sf=cf; */
685 
686  *sf=cf;
687  return resp;
688 } /* SymbProc */
689 
690 int MchlSetup2(int m, chfac** A)
691 {
692  int ierr,i,j,k,lnnz,mnnz;
693  chfac *mf;
694 
695  ierr =CfcAlloc(m,NULL,&mf); if(ierr) return 1;
696  *A=mf;
697 
698  mnnz=m*(m-1)/2;
699  if (!mnnz && m>1)
700  return TRUE;
701 
702  lnnz=mnnz;
703  ierr=iAlloc(mnnz,NULL,&mf->ssub);if(ierr) return 1;
704 
705  mf->snnz=mnnz;
706  k=0;
707  for (i=0; i<m; i++)
708  {
709  mnnz =m-(i+1);
710  mf->shead[i] =k;
711  mf->ssize[i] =mnnz;
712  for (j=0; j<mnnz; j++)
713  mf->ssub[k+j]=(j+1+i);
714  k +=mnnz;
715  mf->perm[i]=i;
716  }
717 
718 
719  k=ChlSymb(mf,lnnz);
720 
721  iFree(&mf->ssub);
722  iFree(&mf->shead);
723  iFree(&mf->ssize);
724 
725  /* This part is different */
726  mf->alldense=1;
727  iFree(&mf->invp);
728  mf->invp=mf->perm;
729 
730  iFree(&mf->ujbeg);
731  mf->ujbeg=mf->perm;
732 
733  iFree(&mf->usub);
734  mf->usub=mf->perm+1;
735 
736  ierr=LvalAlloc(mf,"cf, PspSymb");if(ierr) return 1;
737 
738  return 0;
739 
740 } /* MchlSetup2 */