DSDP
sdpsymb.c
1#include "numchol.h"
2
3int 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
22void 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
31int 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
41static 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
53void 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
183static 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
302static 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
317static void iSort(int n,
318 int *x)
319{
320 qsort((void *)x,n,sizeof(int),CompIntElem);
321} /* iSort */
322
323static 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
393static 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
625int 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
690int 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 */