DSDP
sdpmatx.c
1#include "numchol.h"
2void dCopy(int,double*,double*);
3
4static 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
32static 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
79void 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
175void 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 */
242void 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
253void 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 */
266void 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
280void 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
294void 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
315void 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
421static 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
452void 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
571int 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
579int 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
590int 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
600int 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
610int 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
619int 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
647static 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
661static 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
671static 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
686static 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
697static 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
711int 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
733int 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
749int 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
760int 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
772int 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
809int 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
818int 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