DSDP
sdpmatx.c
1 #include "numchol.h"
2 void dCopy(int,double*,double*);
3 
4 static void SolFwdSnode(chfac *sf,
5  int snde,
6  int f,
7  int l,
8  double x[])
9 {
10  int i,t,sze,*ls,*subg=sf->subg,
11  *ujbeg=sf->ujbeg,*uhead=sf->uhead,
12  *usub=sf->usub;
13  double xi,*l1,*diag=sf->diag,*uval=sf->uval;
14 
15  f += subg[snde];
16  l += subg[snde];
17 
18  for(i=f; i<l; ++i)
19  {
20  x[i] /= diag[i];
21  xi = x[i];
22 
23  ls = usub+ujbeg[i];
24  l1 = uval+uhead[i];
25  sze = l-i-1;
26 
27  for(t=0; t<sze; ++t)
28  x[ls[t]] -= l1[t]*xi;
29  }
30 } /* SolFwdSnode */
31 
32 static void SolBward(int nrow,
33  double diag[],
34  double uval[],
35  int fir[],
36  double x[])
37 {
38  int i,t,sze;
39  double x1,x2,rtemp,
40  *x0,*l1,*l2;
41 
42  for(i=nrow; i;) {
43  for(; i>1; --i) {
44  -- i;
45  l1 = uval+fir[i-1]+1;
46  l2 = uval+fir[i ]+0;
47  sze = nrow-i-1;
48  x1 = 0.0;
49  x2 = 0.0;
50  x0 = x+1+i;
51 
52  for(t=0; t<sze; ++t)
53  {
54  rtemp = x0[t];
55 
56  x1 += l1[t]*rtemp;
57  x2 += l2[t]*rtemp;
58  }
59 
60  x[i] -= x2/diag[i];
61  x[i-1] -= (uval[fir[i-1]]*x[i]+x1)/diag[i-1];
62  }
63 
64  for(; i;) {
65  -- i;
66  l1 = uval+fir[i];
67  sze = nrow-i-1;
68  x1 = 0.0;
69  x0 = x+1+i;
70 
71  for(t=0; t<sze; ++t)
72  x1 += l1[t]*x0[t];
73 
74  x[i] -= x1/diag[i];
75  }
76  }
77 } /* SolBward */
78 
79 void ChlSolveForwardPrivate(chfac *sf,
80  double x[])
81 {
82  int k,s,t,sze,f,l,itemp,*ls,
83  *subg=sf->subg,*ujsze=sf->ujsze,*usub=sf->usub,
84  *ujbeg=sf->ujbeg,*uhead=sf->uhead;
85  double rtemp1,rtemp2,rtemp3,rtemp4,
86  rtemp5,rtemp6,rtemp7,rtemp8,
87  *l1,*l3,*l2,*l4,*l5,*l6,*l7,*l8,
88  *uval=sf->uval;
89 
90  for(s=0; s<sf->nsnds; ++s) {
91  f = subg[s];
92  l = subg[s+1];
93 
94  SolFwdSnode(sf,s,0,l-f,x);
95 
96  itemp = l-f-1;
97  ls = usub+ujbeg[f]+itemp;
98  sze = ujsze[f]-itemp;
99  k = f;
100 
101  itemp = l-1;
102  for(; k+7<l; k+=8) {
103  l1 = uval+uhead[k+0]+itemp-(k+0);
104  l2 = uval+uhead[k+1]+itemp-(k+1);
105  l3 = uval+uhead[k+2]+itemp-(k+2);
106  l4 = uval+uhead[k+3]+itemp-(k+3);
107  l5 = uval+uhead[k+4]+itemp-(k+4);
108  l6 = uval+uhead[k+5]+itemp-(k+5);
109  l7 = uval+uhead[k+6]+itemp-(k+6);
110  l8 = uval+uhead[k+7]+itemp-(k+7);
111 
112  rtemp1 = x[k+0];
113  rtemp2 = x[k+1];
114  rtemp3 = x[k+2];
115  rtemp4 = x[k+3];
116  rtemp5 = x[k+4];
117  rtemp6 = x[k+5];
118  rtemp7 = x[k+6];
119  rtemp8 = x[k+7];
120 
121  for(t=0; t<sze; ++t)
122  x[ls[t]] -= rtemp1*l1[t]
123  + rtemp2*l2[t]
124  + rtemp3*l3[t]
125  + rtemp4*l4[t]
126  + rtemp5*l5[t]
127  + rtemp6*l6[t]
128  + rtemp7*l7[t]
129  + rtemp8*l8[t];
130  }
131 
132  for(; k+3<l; k+=4) {
133  l1 = uval+uhead[k+0]+itemp-(k+0);
134  l2 = uval+uhead[k+1]+itemp-(k+1);
135  l3 = uval+uhead[k+2]+itemp-(k+2);
136  l4 = uval+uhead[k+3]+itemp-(k+3);
137 
138  rtemp1 = x[k+0];
139  rtemp2 = x[k+1];
140  rtemp3 = x[k+2];
141  rtemp4 = x[k+3];
142 
143  for(t=0; t<sze; ++t)
144  x[ls[t]] -= rtemp1*l1[t]
145  + rtemp2*l2[t]
146  + rtemp3*l3[t]
147  + rtemp4*l4[t];
148  }
149 
150  for(; k+1<l; k+=2) {
151  l1 = uval+uhead[k+0]+itemp-(k+0);
152  l2 = uval+uhead[k+1]+itemp-(k+1);
153 
154  rtemp1 = x[k+0];
155  rtemp2 = x[k+1];
156 
157  for(t=0; t<sze; ++t)
158  x[ls[t]] -= rtemp1*l1[t]
159  + rtemp2*l2[t];
160  }
161 
162 
163  for(; k<l; ++k) {
164  l1 = uval+uhead[k+0]+itemp-(k+0);
165 
166  rtemp1 = x[k+0];
167 
168  for(t=0; t<sze; ++t)
169  x[ls[t]] -= rtemp1*l1[t];
170  }
171  }
172 
173 }
174 
175 void ChlSolveBackwardPrivate(chfac *sf,
176  double x[],
177  double b[])
178 {
179  /* Note: x: input, or left hand side b: output, or solution */
180 
181  int i,s,t,sze,f,l,*ls,
182  *subg=sf->subg,*ujsze=sf->ujsze,*usub=sf->usub,
183  *ujbeg=sf->ujbeg,*uhead=sf->uhead;
184  double x1,x2,*l1,*l2,rtemp1,
185  *diag=sf->diag,*uval=sf->uval;
186 
187 
188  if (sf->nsnds) {
189  s = sf->nsnds - 1;
190  f = subg[s];
191  l = subg[s+1];
192 
193  dCopy(l-f,x+f,b+f);
194 
195  SolBward(l-f,diag+f,uval,uhead+f,b+f);
196 
197  s = sf->nsnds-1;
198 
199  for(; s>=1; --s) {
200  f = subg[s-1];
201  l = subg[s];
202  i = l;
203 
204  for(; i>1+f; --i) {
205  -- i;
206  ls = usub+ujbeg[i];
207  l1 = uval+uhead[i-1]+1;
208  l2 = uval+uhead[i ]+0;
209  sze = ujsze[i];
210  x1 = 0.0;
211  x2 = 0.0;
212 
213  for(t=0; t<sze; ++t) {
214  rtemp1 = b[ls[t]];
215 
216  x1 += l1[t]*rtemp1;
217  x2 += l2[t]*rtemp1;
218  }
219 
220  b[i] = x[i ] - x2 / diag[i];
221  b[i-1] = x[i-1] - (x1 + uval[uhead[i-1]]*b[i]) / diag[i-1];
222  }
223 
224  for(; i>f;) {
225  -- i;
226  l1 = uval+uhead[i];
227  ls = usub+ujbeg[i];
228  sze = ujsze[i];
229  x1 = 0.0;
230 
231  for(t=0; t<sze; ++t)
232  x1+= l1[t]*b[ls[t]];
233 
234  b[i] = x[i] - x1/diag[i];
235  }
236  }
237  }
238 
239 }
240 
241 /* Everything right for permuted system */
242 void ChlSolveForward2(chfac *sf,
243  double b[],
244  double x[]){
245  int i,nrow=sf->nrow;
246  double *sqrtdiag=sf->sqrtdiag;
247  ChlSolveForwardPrivate(sf,b);
248  for(i=0; i<nrow; ++i){
249  x[i] = b[i]*sqrtdiag[i]; /* x[i] = b[i]*sqrt(sf->diag[i]); */
250  }
251 }
252 
253 void ChlSolveBackward2(chfac *sf,
254  double b[],
255  double x[]){
256  int i,nrow=sf->nrow;
257  double *sqrtdiag=sf->sqrtdiag;
258  for(i=0; i<nrow; ++i){
259  x[i] = b[i]/(sqrtdiag[i]); /* x[i] = b[i]/sqrt(sf->diag[i]) ; */
260  }
261  ChlSolveBackwardPrivate(sf,x,b);
262  memcpy(x,b,nrow*sizeof(double));
263 }
264 
265 /* These routines together will solve an equation correctly */
266 void ChlSolveForward(chfac *sf,
267  double b[],
268  double x[]){
269  int i,nrow=sf->nrow,*perm=sf->perm;
270  double *w=sf->rw,*sqrtdiag=sf->sqrtdiag;
271  for(i=0; i<nrow; ++i)
272  w[i] = b[perm[i]];
273  ChlSolveForwardPrivate(sf,w);
274  for(i=0; i<nrow; ++i){
275  x[i] = w[i]*sqrtdiag[i]; /* x[i] = w[i]*sqrt(sf->diag[i]); */
276  }
277 
278 }
279 
280 void ChlSolveBackward(chfac *sf,
281  double b[],
282  double x[]){
283  int i,nrow=sf->nrow,*invp=sf->invp;
284  double *w=sf->rw,*sqrtdiag=sf->sqrtdiag;
285  for(i=0; i<nrow; ++i){
286  x[i] = b[i]/sqrtdiag[i];
287  }
288  ChlSolveBackwardPrivate(sf,x,w);
289  for(i=0; i<nrow; ++i)
290  x[i] = w[invp[i]];
291 
292 }
293 
294 void ChlSolve(chfac *sf,
295  double b[],
296  double x[]){
297  int i,nrow=sf->nrow,*perm=sf->perm,*invp=sf->invp;
298  double *rw=sf->rw;
299  /*
300  ChlSolveForward(sf,b,w,x);
301  ChlSolveBackward(sf,w,x,b);
302  */
303  for(i=0; i<nrow; ++i)
304  x[i] = b[perm[i]];
305 
306  ChlSolveForwardPrivate(sf,x);
307  ChlSolveBackwardPrivate(sf,x,rw);
308 
309  for(i=0; i<nrow; ++i)
310  x[i] = rw[invp[i]]; /* x[i] = b[i]; */
311  return;
312 }
313 
314 
315 void ForwSubst(chfac *sf,
316  double b[],
317  double x[])
318 {
319  int i,k,s,t,sze,f,l,itemp,*ls,
320  *subg=sf->subg,*ujsze=sf->ujsze,*usub=sf->usub,
321  *ujbeg=sf->ujbeg,*uhead=sf->uhead;
322  double rtemp1,rtemp2,rtemp3,rtemp4,
323  rtemp5,rtemp6,rtemp7,rtemp8,
324  *l1,*l3,*l2,*l4,*l5,*l6,*l7,*l8,
325  *diag=sf->diag,*uval=sf->uval;
326 
327  for(i=0; i<sf->nrow; ++i)
328  x[i] = b[sf->perm[i]];
329 
330  for(s=0; s<sf->nsnds; ++s) {
331  f = subg[s];
332  l = subg[s+1];
333 
334  SolFwdSnode(sf,s,0,l-f,x);
335 
336  itemp = l-f-1;
337  ls = usub+ujbeg[f]+itemp;
338  sze = ujsze[f]-itemp;
339  k = f;
340 
341  itemp = l-1;
342  for(; k+7<l; k+=8) {
343  l1 = uval+uhead[k+0]+itemp-(k+0);
344  l2 = uval+uhead[k+1]+itemp-(k+1);
345  l3 = uval+uhead[k+2]+itemp-(k+2);
346  l4 = uval+uhead[k+3]+itemp-(k+3);
347  l5 = uval+uhead[k+4]+itemp-(k+4);
348  l6 = uval+uhead[k+5]+itemp-(k+5);
349  l7 = uval+uhead[k+6]+itemp-(k+6);
350  l8 = uval+uhead[k+7]+itemp-(k+7);
351 
352  rtemp1 = x[k+0];
353  rtemp2 = x[k+1];
354  rtemp3 = x[k+2];
355  rtemp4 = x[k+3];
356  rtemp5 = x[k+4];
357  rtemp6 = x[k+5];
358  rtemp7 = x[k+6];
359  rtemp8 = x[k+7];
360 
361  for(t=0; t<sze; ++t)
362  x[ls[t]] -= rtemp1*l1[t]
363  + rtemp2*l2[t]
364  + rtemp3*l3[t]
365  + rtemp4*l4[t]
366  + rtemp5*l5[t]
367  + rtemp6*l6[t]
368  + rtemp7*l7[t]
369  + rtemp8*l8[t];
370  }
371 
372  for(; k+3<l; k+=4) {
373  l1 = uval+uhead[k+0]+itemp-(k+0);
374  l2 = uval+uhead[k+1]+itemp-(k+1);
375  l3 = uval+uhead[k+2]+itemp-(k+2);
376  l4 = uval+uhead[k+3]+itemp-(k+3);
377 
378  rtemp1 = x[k+0];
379  rtemp2 = x[k+1];
380  rtemp3 = x[k+2];
381  rtemp4 = x[k+3];
382 
383  for(t=0; t<sze; ++t)
384  x[ls[t]] -= rtemp1*l1[t]
385  + rtemp2*l2[t]
386  + rtemp3*l3[t]
387  + rtemp4*l4[t];
388  }
389 
390  for(; k+1<l; k+=2) {
391  l1 = uval+uhead[k+0]+itemp-(k+0);
392  l2 = uval+uhead[k+1]+itemp-(k+1);
393 
394  rtemp1 = x[k+0];
395  rtemp2 = x[k+1];
396 
397  for(t=0; t<sze; ++t)
398  x[ls[t]] -= rtemp1*l1[t]
399  + rtemp2*l2[t];
400  }
401 
402 
403  for(; k<l; ++k) {
404  l1 = uval+uhead[k+0]+itemp-(k+0);
405 
406  rtemp1 = x[k+0];
407 
408  for(t=0; t<sze; ++t)
409  x[ls[t]] -= rtemp1*l1[t];
410  }
411  }
412 
413  for (i=0; i<sf->nrow; i++){
414  x[i] = x[i] * sqrt( fabs(diag[i]) );
415  }
416 
417 } /* ForwSubst */
418 
419 
420 
421 static void mulSnod(chfac *sf,
422  int snde,
423  int f,
424  int l,
425  double *b,
426  double *x)
427 {
428  int i,t,sze,*ls,*subg,*ujbeg,*uhead,*usub;
429  double xi,*l1,*diag,*uval;
430 
431  subg =sf->subg;
432  ujbeg=sf->ujbeg;
433  uhead=sf->uhead;
434  usub =sf->usub;
435  diag =sf->diag;
436  uval =sf->uval;
437 
438  f += subg[snde];
439  l += subg[snde];
440 
441  for(i=f; i<l; ++i) {
442  xi =b[i];
443  ls =usub+ujbeg[i];
444  l1 =uval+uhead[i];
445  sze =l-i-1;
446  x[i]+=xi*diag[i];
447  for(t=0; t<sze; ++t)
448  x[ls[t]]+=l1[t]*xi;
449  }
450 } /* mulSnod */
451 
452 void GetUhat(chfac *sf,
453  double *b,
454  double *x)
455  /* If S = L L^T, then b = L x */
456 {
457  int i,k,n,s,t,sze,f,l,itemp,*ls,
458  *subg,*ujsze,*usub,*ujbeg,*uhead;
459  double rtemp1,rtemp2,rtemp3,rtemp4,
460  rtemp5,rtemp6,rtemp7,rtemp8,
461  *l1,*l3,*l2,*l4,*l5,*l6,*l7,*l8,
462  *diag,*uval;
463 
464  n =sf->nrow;
465  subg =sf->subg;
466  ujsze=sf->ujsze;
467  usub =sf->usub;
468  ujbeg=sf->ujbeg;
469  uhead=sf->uhead;
470  diag =sf->diag;
471  uval =sf->uval;
472 
473  for (i=0; i<n; i++) {
474  if (diag[i]>0)
475  x[i]=b[i]/sqrt(diag[i]);
476  else x[i]=b[i]/sqrt(-diag[i]);
477  b[i]=0.0;
478  }
479 
480  for (s=0; s<sf->nsnds; s++) {
481  f=subg[s];
482  l=subg[s+1];
483 
484  mulSnod(sf,s,0,l-f,x,b);
485 
486  itemp=l-f-1;
487  ls =usub+ujbeg[f]+itemp;
488  sze =ujsze[f]-itemp;
489  k =f;
490 
491  itemp=l-1;
492  for(; k+7<l; k+=8) {
493  l1 =uval+uhead[k+0]+itemp-(k+0);
494  l2 =uval+uhead[k+1]+itemp-(k+1);
495  l3 =uval+uhead[k+2]+itemp-(k+2);
496  l4 =uval+uhead[k+3]+itemp-(k+3);
497  l5 =uval+uhead[k+4]+itemp-(k+4);
498  l6 =uval+uhead[k+5]+itemp-(k+5);
499  l7 =uval+uhead[k+6]+itemp-(k+6);
500  l8 =uval+uhead[k+7]+itemp-(k+7);
501 
502  rtemp1=x[k+0];
503  rtemp2=x[k+1];
504  rtemp3=x[k+2];
505  rtemp4=x[k+3];
506  rtemp5=x[k+4];
507  rtemp6=x[k+5];
508  rtemp7=x[k+6];
509  rtemp8=x[k+7];
510 
511  for(t=0; t<sze; ++t)
512  b[ls[t]]+= rtemp1*l1[t]
513  +rtemp2*l2[t]
514  +rtemp3*l3[t]
515  +rtemp4*l4[t]
516  +rtemp5*l5[t]
517  +rtemp6*l6[t]
518  +rtemp7*l7[t]
519  +rtemp8*l8[t];
520  }
521 
522 
523  for(; k+3<l; k+=4) {
524  l1 =uval+uhead[k+0]+itemp-(k+0);
525  l2 =uval+uhead[k+1]+itemp-(k+1);
526  l3 =uval+uhead[k+2]+itemp-(k+2);
527  l4 =uval+uhead[k+3]+itemp-(k+3);
528 
529  rtemp1=x[k+0];
530  rtemp2=x[k+1];
531  rtemp3=x[k+2];
532  rtemp4=x[k+3];
533 
534  for(t=0; t<sze; ++t)
535  b[ls[t]]+= rtemp1*l1[t]
536  +rtemp2*l2[t]
537  +rtemp3*l3[t]
538  +rtemp4*l4[t];
539  }
540 
541  for(; k+1<l; k+=2) {
542  l1 =uval+uhead[k+0]+itemp-(k+0);
543  l2 =uval+uhead[k+1]+itemp-(k+1);
544 
545  rtemp1=x[k+0];
546  rtemp2=x[k+1];
547 
548  for(t=0; t<sze; ++t)
549  b[ls[t]]+= rtemp1*l1[t]
550  +rtemp2*l2[t];
551  }
552 
553 
554  for(; k<l; ++k) {
555  l1 =uval+uhead[k+0]+itemp-(k+0);
556 
557  rtemp1=x[k+0];
558 
559  for(t=0; t<sze; ++t)
560  b[ls[t]]+=rtemp1*l1[t];
561  }
562  }
563 
564  for (i=0; i<n; i++)
565  x[ sf->invp[i] ]=b[i];
566 } /* GetUhat */
567 
568 
569 
570 
571 int MatSolve4(chfac*sf, double *b, double *x,int n){
572 
573  memcpy(x,b,n*sizeof(double));
574  ChlSolve(sf, b, x);
575 
576  return 0;
577 }
578 
579 int Mat4GetDiagonal(chfac*sf, double *b,int n){
580 
581  int i,*invp=sf->invp;
582  double *diag=sf->diag;
583 
584  for (i=0; i<n; i++, invp++){
585  b[i]=diag[*invp];
586  }
587  return 0;
588 }
589 
590 int Mat4SetDiagonal(chfac*sf, double *b,int n){
591 
592  int i,*invp=sf->invp;
593  double *diag=sf->diag;
594  for (i=0; i<n; i++){
595  diag[invp[i]]=b[i];
596  }
597  return 0;
598 }
599 
600 int Mat4AddDiagonal(chfac*sf, double *b,int n){
601 
602  int i,*invp=sf->invp;
603  double *diag=sf->diag;
604  for (i=0; i<n; i++){
605  diag[invp[i]]+=b[i];
606  }
607  return 0;
608 }
609 
610 int MatAddDiagonalElement(chfac*sf, int row, double dd){
611 
612  int *invp=sf->invp;
613  double *diag=sf->diag;
614  diag[invp[row]]+=dd;
615  return 0;
616 }
617 
618 
619 int MatMult4(chfac *sf, double *x, double *y, int n){
620 
621  int i,j,*invp=sf->invp,*perm=sf->perm;
622  int *usub=sf->usub,*ujbeg=sf->ujbeg,*uhead=sf->uhead, *ujsze=sf->ujsze;
623  int *iptr,k1,k2;
624  double dd,*sval,*diag=sf->diag,*uval=sf->uval;
625 
626  for (i=0; i<n; i++){
627  y[i] = diag[invp[i]] * x[i];
628  }
629  for (i=0; i<n; ++i){
630 
631  iptr=usub + ujbeg[i];
632  sval=uval + uhead[i];
633  k1=perm[i];
634  for (j=0; j<ujsze[i]; j++){
635  dd=sval[j];
636  if (fabs(dd)> 1e-15){
637  k2=perm[iptr[j]];
638  y[k1] += dd * x[k2];
639  y[k2] += dd * x[k1];
640  }
641  }
642  }
643  return 0;
644 }
645 
646 
647 static void setXYind2(int nnz,
648  double* y,
649  double* x,
650  int* s,
651  int* invp)
652 {
653  int i;
654 
655  for(i=0; i<nnz; ++i) {
656  x[i]=y[ invp[ s[i] ] ];
657  y[ invp[s[i]] ]=0.0;
658  }
659 } /* setXYind */
660 
661 static void setColi(chfac* cl,
662  int i,
663  double* ai)
664 {
665  setXYind2(cl->ujsze[i],ai,
666  cl->uval+cl->uhead[i],
667  cl->usub+cl->ujbeg[i],
668  cl->perm);
669 } /* setColi */
670 
671 static void setXYind2add(int nnz,
672  double dd,
673  double* y,
674  double* x,
675  int* s,
676  int* invp)
677 {
678  int i;
679 
680  for(i=0; i<nnz; ++i) {
681  x[i]+=dd*y[ invp[ s[i] ] ];
682  y[ invp[s[i]] ]=0.0;
683  }
684 } /* setXYind */
685 
686 static void setColi2(chfac* cl,
687  int i,double dd,
688  double* ai)
689 {
690  setXYind2add(cl->ujsze[i],dd,ai,
691  cl->uval+cl->uhead[i],
692  cl->usub+cl->ujbeg[i],
693  cl->perm);
694 } /* setColi */
695 
696 
697 static void getXYind2(int nnz,
698  double* y,
699  double* x,
700  int* s,
701  int* invp)
702 {
703  int i;
704 
705  for(i=0; i<nnz; ++i) {
706  y[ invp[s[i]] ]=x[i];
707  }
708 } /* setXYind */
709 
710 
711 int Mat4View(chfac *sf){
712 
713  int i,j,n=sf->nrow;
714  double *v=sf->rw;
715  for (i=0; i<n; i++){
716  for (j=0;j<n;++j) v[j]=0.0;
717  getXYind2(sf->ujsze[i],v,
718  sf->uval+sf->uhead[i],
719  sf->usub+sf->ujbeg[i],
720  sf->perm);
721  v[i]=sf->diag[sf->invp[i]];
722  printf("Row %d, ",i);
723  for (j=0;j<n;j++){
724  if (v[j]!=0) printf(" %d: %4.4e ",j,v[j]);
725  }
726  printf("\n");
727  }
728 
729  return 0;
730 
731 }
732 
733 int MatZeroEntries4(chfac *sf){
734 
735  int i,n=sf->n;
736  double *rw=sf->rw;
737  memset((void*)(sf->diag),0,n*sizeof(double));
738  memset((void*)(rw),0,n*sizeof(double));
739 
740  for (i=0; i<n; i++){
741  setColi(sf,i,rw);
742  }
743 
744  return 0;
745 }
746 
747 
748 
749 int MatSetColumn4(chfac *cl, double *val, int col){
750 
751  int pcol=cl->invp[col];
752 
753  cl->diag[pcol]=val[col];
754  val[col]=0.0;
755  setColi(cl,pcol,val);
756 
757  return 0;
758 } /* SetColumn */
759 
760 int MatAddColumn4(chfac *cl, double dd, double *val, int col){
761 
762  int pcol=cl->invp[col];
763 
764  cl->diag[pcol]+=dd*val[col];
765  val[col]=0.0;
766  setColi2(cl,pcol,dd,val);
767 
768  return 0;
769 } /* SetColumn */
770 
771 
772 int MatSetValue4(chfac *cl, int row,int col,double val, int setmode){
773 
774  int i;
775  double* x=cl->uval+cl->uhead[col];
776  int* s=cl->usub+cl->ujbeg[col];
777  int nnz=cl->ujsze[col];
778  int insertmode=1,addmode=2;
779 
780  if (row<0 || col<0 || row>=cl->n || col>=cl->n){
781  printf("CHol set Value error: Row: %d, COl: %d \n",row,col);
782  return 1;
783  }
784 
785  if (setmode==insertmode&&row==col){
786  cl->diag[cl->invp[col]]=val;
787  } else if (setmode==addmode&&row==col) {
788  cl->diag[cl->invp[col]]+=val;
789  } else if (setmode==insertmode){
790  for(i=0; i<nnz; ++i) {
791  if (s[i]==row){
792  x[i]=val;
793  }
794  }
795  } else if (setmode==addmode){
796  for(i=0; i<nnz; ++i) {
797  if (s[i]==row){
798  x[i]+=val;
799  }
800  }
801  } else {
802  return 1;
803  }
804 
805  return 0;
806 } /* SetValue */
807 
808 
809 int Mat4DiagonalShift(chfac*sf, double shift){
810  int i,n=sf->nrow;
811  double *diag=sf->diag;
812  for (i=0; i<n; i++){
813  diag[i] += shift;
814  }
815  return 0;
816 }
817 
818 int Mat4LogDet(chfac*sf, double *dd){
819  int i,n=sf->nrow;
820  double *diag=sf->diag,ddd=0;
821  for (i=0; i<n; i++){
822  if (diag[i]<=0) return 1;
823  ddd+=log(diag[i]);
824  }
825  *dd=ddd;
826  return 0;
827 }
828