DSDP
sdpvec.c
Go to the documentation of this file.
1#include "dsdpsys.h"
2#include "dsdpvec.h"
3#include "dsdplapack.h"
8#if !defined (min)
9#define min(a,b) ((a <= b)? (a) : (b))
10#endif
11#if !defined (max)
12#define max(a,b) ((a >= b)? (a) : (b))
13#endif
14
15#define DSPPVecCheck(a,b) {if (a.dim != b.dim) return 1; if (a.dim>0 && (a.val==NULL || b.val==NULL) ) return 2;}
16
17static int nvecs=0;
18#undef __FUNCT__
19#define __FUNCT__ "DSDPVecCreateSeq"
20int DSDPVecCreate(DSDPVec *V){
21 int info;
22 info = DSDPVecCreateSeq(0,V);DSDPCHKERR(info);
23 return 0;
24}
25
26#undef __FUNCT__
27#define __FUNCT__ "DSDPVecCreateSeq"
28int DSDPVecCreateSeq(int n ,DSDPVec *V){
29 int info;
30 V->dim=n;
31 if (n>0){
32 nvecs++;
33 DSDPCALLOC2(&(V->val),double,n,&info);DSDPCHKERR(info);
34 if (V->val==NULL) return 1;
35 } else {
36 V->val=NULL;
37 }
38 return 0;
39}
40/*
41#undef __FUNCT__
42#define __FUNCT__ "DSDPVecCreateWArray"
43int DSDPVecCreateWArray(DSDPVec *V, double* vv, int n){
44 V->dim=n;
45 if (n>0){
46 V->val=vv;
47 } else {
48 V->val=NULL;
49 }
50 return 0;
51}
52*/
53#undef __FUNCT__
54#define __FUNCT__ "DSDPVecDestroy"
55int DSDPVecDestroy(DSDPVec *V){
56 int info;
57 if ((*V).val){
58 DSDPFREE(&(*V).val,&info);DSDPCHKERR(info);
59 nvecs--;
60 }
61
62 (*V).dim=0;
63 (*V).val=0;
64 return 0;
65}
66
67/*
68int DSDPVecGetSize(DSDPVec V, int *n){
69 *n=V.dim;
70 return 0;
71}
72
73int DSDPVecGetArray(DSDPVec V, double **dptr){
74 *dptr=V.val;
75 return 0;
76}
77
78int DSDPVecRestoreArray(DSDPVec V, double **dptr){
79 *dptr=0;
80 return 0;
81}
82*/
83
84
85#undef __FUNCT__
86#define __FUNCT__ "DSDPVecView"
87int DSDPVecView(DSDPVec vec){
88 int i;
89 for (i=0; i<vec.dim; i++){
90 printf("%3.3e ",vec.val[i]);
91 }
92 printf("\n");
93 return 0;
94}
95
96#undef __FUNCT__
97#define __FUNCT__ "DSDPVecISet"
98int DSDPVecISet(int* ival,DSDPVec V){
99 int i;
100 for (i=0;i<V.dim;i++){
101 V.val[i]=ival[i];
102 }
103 return 0;
104}
105
106#undef __FUNCT__
107#define __FUNCT__ "DSDPVecSetValue"
108int DSDPVecSetValue(DSDPVec V,int row,double value){
109 V.val[row]=value;
110 return 0;
111}
112
113#undef __FUNCT__
114#define __FUNCT__ "DSDPVecZero"
115int DSDPVecZero(DSDPVec V){
116 int n=V.dim;
117 double *v=V.val;
118 memset((void*)v,0,n*sizeof(double));
119 return 0;
120}
121
122
123#undef __FUNCT__
124#define __FUNCT__ "DSDPVecNormalize"
125int DSDPVecNormalize(DSDPVec V){
126 int info;
127 double vnorm;
128 info = DSDPVecNorm2(V,&vnorm);DSDPCHKERR(info);
129 if (vnorm==0){ return 1;}
130 vnorm=1.0/(vnorm);
131 info = DSDPVecScale(vnorm,V);DSDPCHKERR(info);
132 return 0;
133}
134
135#undef __FUNCT__
136#define __FUNCT__ "DSDPVecSetBasis"
137int DSDPVecSetBasis(DSDPVec V,int row){
138 int info;
139 info=DSDPVecZero(V);
140 V.val[row]=1.0;
141 return 0;
142}
143
144#undef __FUNCT__
145#define __FUNCT__ "DSDPVecCopy"
146int DSDPVecCopy( DSDPVec v1, DSDPVec v2){
147
148 int n=v1.dim;
149 double *val1=v1.val,*val2=v2.val;
150 DSPPVecCheck(v1,v2);
151 if (val1!=val2){
152 memcpy(val2,val1,n*sizeof(double));
153 }
154 return 0;
155}
156
157
158#undef __FUNCT__
159#define __FUNCT__ "DSDPVecSum"
160int DSDPVecSum( DSDPVec v, double *vnorm){
161 int i,n;
162 n = v.dim;
163 *vnorm = 0.0;
164 for (i=0; i<n; i++){
165 *vnorm += v.val[i];
166 }
167 if (*vnorm!=*vnorm) return 1;
168 return 0;
169}
170#undef __FUNCT__
171#define __FUNCT__ "DSDPVecNorm1"
172int DSDPVecNorm1( DSDPVec v, double *vnorm){
173 ffinteger N=v.dim,INCX=1;
174 *vnorm=0;
175 *vnorm=dasum(&N,v.val,&INCX);
176 if (*vnorm!=*vnorm) return 1;
177 return 0;
178}
179
180
181#undef __FUNCT__
182#define __FUNCT__ "DSDPVecDot"
183int DSDPVecDot(DSDPVec V1, DSDPVec V2, double *ans){
184 ffinteger ione=1, nn=V1.dim;
185 double *v1=V1.val,*v2=V2.val;
186 *ans=ddot(&nn,v1,&ione,v2,&ione);
187 if (*ans!=*ans) return 1;
188 return 0;
189}
190
191
192#undef __FUNCT__
193#define __FUNCT__ "DSDPVecNorm22"
194int DSDPVecNorm22( DSDPVec VV, double *vnorm){
195 ffinteger ione=1,nn=VV.dim;
196 double dd,*v=VV.val;
197 dd=dnrm2(&nn,v,&ione);
198 *vnorm = dd*dd;
199 if (*vnorm!=*vnorm) return 1;
200 return 0;
201}
202#undef __FUNCT__
203#define __FUNCT__ "DSDPVecNorm2"
204int DSDPVecNorm2( DSDPVec VV, double *vnorm){
205 ffinteger ione=1,nn=VV.dim;
206 double dd,*v=VV.val;
207 dd=dnrm2(&nn,v,&ione);
208 *vnorm = dd;
209 if (*vnorm!=*vnorm) return 1;
210 return 0;
211}
212
213#undef __FUNCT__
214#define __FUNCT__ "DSDPVecScale"
215int DSDPVecScale(double alpha, DSDPVec VV){
216 ffinteger ione=1,nn=VV.dim;
217 double *v=VV.val;
218 dscal(&nn,&alpha,v,&ione);
219 return 0;
220}
221
222#undef __FUNCT__
223#define __FUNCT__ "DSDPVecAXPY"
224int DSDPVecAXPY(double alpha, DSDPVec x, DSDPVec y){
225 ffinteger ione=1,nn=x.dim;
226 double *yy=y.val,*xx=x.val;
227 if (alpha==0) return 0;
228 daxpy(&nn,&alpha,xx,&ione,yy,&ione);
229 return 0;
230}
231
232/*
233#undef __FUNCT__
234#define __FUNCT__ "DSDPVecNorm22"
235int DSDPVecNorm22( DSDPVec v, double *vnorm){
236 int i,n=v.dim;
237 double *val=v.val;
238
239 *vnorm = 0.0;
240 for (i=0; i<n; i++){
241 *vnorm += val[i]*val[i];
242 }
243 return 0;
244}
245#undef __FUNCT__
246#define __FUNCT__ "DSDPVecNorm2"
247int DSDPVecNorm2( DSDPVec v, double *vnorm){
248 int info;
249 info=DSDPVecNorm22(v,vnorm); if (info) return 1;
250 if (*vnorm!=*vnorm) return 1;
251 *vnorm = sqrt(*vnorm);
252 return 0;
253}
254*/
255#undef __FUNCT__
256#define __FUNCT__ "DSDPVecNormInfinity"
257int DSDPVecNormInfinity( DSDPVec v, double *vnorm){
258
259 int i,n=v.dim;
260 double *val=v.val;
261
262 *vnorm = 0.0;
263
264 for (i=0; i<n; i++){
265 *vnorm = max(*vnorm,fabs(val[i]));
266 }
267 if (*vnorm!=*vnorm) return 1;
268
269 return 0;
270}
271/*
272#undef __FUNCT__
273#define __FUNCT__ "DSDPVecScale"
274int DSDPVecScale(double alpha, DSDPVec x){
275 int i,ii,n;
276 double *xx=x.val;
277 n=x.dim;
278
279 for (ii=0; ii<n/4; ++ii){
280 i=ii*4;
281 xx[i]*= alpha;
282 xx[i+1]*= alpha;
283 xx[i+2]*= alpha;
284 xx[i+3]*= alpha;
285 }
286 for (i=4*(n/4); i<n; ++i){
287 xx[i]*= alpha;
288 }
289 return 0;
290}
291
292#undef __FUNCT__
293#define __FUNCT__ "DSDPVecAXPY"
294int DSDPVecAXPY(double alpha, DSDPVec x, DSDPVec y){
295
296 int i,ii,n=x.dim;
297 double *yy=y.val,*xx=x.val;
298
299 DSPPVecCheck(x,y);
300
301 for (ii=0; ii<n/4; ++ii){
302 i=ii*4;
303 yy[i] += (alpha)*xx[i];
304 yy[i+1] += (alpha)*xx[i+1];
305 yy[i+2] += (alpha)*xx[i+2];
306 yy[i+3] += (alpha)*xx[i+3];
307 }
308 for (i=4*(n/4); i<n; ++i){
309 yy[i] += (alpha)*xx[i];
310 }
311
312 return 0;
313}
314*/
315#undef __FUNCT__
316#define __FUNCT__ "DSDPVecWAXPBY"
317int DSDPVecWAXPBY(DSDPVec w, double alpha, DSDPVec x, double beta, DSDPVec y){
318
319 int i,ii,n=x.dim;
320 double *yy=y.val,*xx=x.val,*ww=w.val;
321 DSPPVecCheck(x,y);
322 DSPPVecCheck(x,w);
323
324 for (ii=0; ii<n/4; ++ii){
325 i=ii*4;
326 ww[i] = (alpha)*xx[i] + (beta)*yy[i];
327 ww[i+1] = (alpha)*xx[i+1] + (beta)*yy[i+1];
328 ww[i+2] = (alpha)*xx[i+2] + (beta)*yy[i+2];
329 ww[i+3] = (alpha)*xx[i+3] + (beta)*yy[i+3];
330 }
331 for (i=4*(n/4); i<n; ++i){
332 ww[i] = (alpha)*xx[i] + (beta)*yy[i];
333 }
334
335 return 0;
336}
337
338#undef __FUNCT__
339#define __FUNCT__ "DSDPVecWAXPY"
340int DSDPVecWAXPY(DSDPVec w,double alpha, DSDPVec x, DSDPVec y){
341 int info;
342 info=DSDPVecCopy(y,w);
343 info=DSDPVecAXPY(alpha,x,w);
344 return 0;
345}
346
347
348#undef __FUNCT__
349#define __FUNCT__ "DSDPVecAYPX"
350int DSDPVecAYPX(double alpha, DSDPVec x, DSDPVec y){
351
352 int i,ii,n=x.dim;
353 double *yy=y.val,*xx=x.val;
354
355 DSPPVecCheck(x,y);
356 for (ii=0; ii<n/4; ++ii){
357 i=ii*4;
358 yy[i] = xx[i]+(alpha)*yy[i];
359 yy[i+1] = xx[i+1]+(alpha)*yy[i+1];
360 yy[i+2] = xx[i+2]+(alpha)*yy[i+2];
361 yy[i+3] = xx[i+3]+(alpha)*yy[i+3];
362 }
363 for (i=4*(n/4); i<n; ++i){
364 yy[i] = xx[i]+(alpha)*yy[i];
365 }
366
367 return 0;
368}
369
370#undef __FUNCT__
371#define __FUNCT__ "DSDPVecAYPX"
372int DSDPVecScaleCopy(DSDPVec x, double alpha, DSDPVec y){
373
374 int i,ii,n=x.dim;
375 double *yy=y.val,*xx=x.val;
376
377 DSPPVecCheck(x,y);
378 for (ii=0; ii<n/4; ++ii){
379 i=ii*4;
380 yy[i] = (alpha)*xx[i];
381 yy[i+1] = (alpha)*xx[i+1];
382 yy[i+2] = (alpha)*xx[i+2];
383 yy[i+3] = (alpha)*xx[i+3];
384 }
385 for (i=4*(n/4); i<n; ++i){
386 yy[i] = (alpha)*xx[i];
387 }
388
389 return 0;
390}
391
392#undef __FUNCT__
393#define __FUNCT__ "DSDPVecDuplicate"
394int DSDPVecDuplicate(DSDPVec V1,DSDPVec *V2){
395 int info,n=V1.dim;
396 info = DSDPVecCreateSeq(n ,V2);DSDPCHKERR(info);
397 return 0;
398}
399
400
401#undef __FUNCT__
402#define __FUNCT__ "DSDPVecSet"
403int DSDPVecSet(double alpha, DSDPVec V){
404
405 int i,ii,n=V.dim;
406 double *val=V.val;
407
408 if (alpha==0.0){
409 memset((void*)val,0,n*sizeof(double));
410 return 0;
411 }
412 for (ii=0; ii<n/4; ++ii){
413 i=ii*4;
414 val[i] = val[i+1] = val[i+2] = val[i+3] = alpha;
415 }
416 for (i=4*(n/4); i<n; ++i){
417 val[i]= alpha;
418 }
419 return 0;
420}
421
422/*
423#undef __FUNCT__
424#define __FUNCT__ "DSDPVecDot"
425int DSDPVecDot(DSDPVec V1, DSDPVec V2, double *ans){
426
427 int i,ii,m=V1.dim;
428 double *v1=V1.val,*v2=V2.val;
429
430 DSPPVecCheck(V1,V2);
431 *ans=0.0;
432 for (ii=0; ii<m/4; ++ii){
433 i=ii*4;
434 *ans += v1[i]*v2[i] + v1[i+1]*v2[i+1] + v1[i+2]*v2[i+2] + v1[i+3]*v2[i+3] ;
435 }
436 for (i=4*(m/4); i<m; ++i){
437 *ans += v1[i]*v2[i];
438 }
439 if (*ans!=*ans) return 1;
440
441 return 0;
442}
443*/
444
445#undef __FUNCT__
446#define __FUNCT__ "DSDPVecPointwiseMin"
447int DSDPVecPointwiseMin( DSDPVec V1, DSDPVec V2, DSDPVec V3){
448
449 int i,n=V1.dim;
450 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
451
452 DSPPVecCheck(V1,V3);
453 DSPPVecCheck(V2,V3);
454 for (i=0; i<n; ++i){
455 v3[i]=DSDPMin(v2[i],v1[i]);
456 }
457
458 return 0;
459}
460
461#undef __FUNCT__
462#define __FUNCT__ "DSDPVecPointwiseMax"
463int DSDPVecPointwiseMax( DSDPVec V1, DSDPVec V2, DSDPVec V3){
464
465 int i,n=V1.dim;
466 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
467
468 DSPPVecCheck(V1,V3);
469 DSPPVecCheck(V2,V3);
470 for (i=0; i<n; ++i){
471 v3[i]=DSDPMax(v2[i],v1[i]);
472 }
473 return 0;
474}
475
476#undef __FUNCT__
477#define __FUNCT__ "DSDPVecPointwiseMult"
478int DSDPVecPointwiseMult( DSDPVec V1, DSDPVec V2, DSDPVec V3){
479
480 int ii,i,n=V1.dim;
481 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
482
483 DSPPVecCheck(V1,V3);
484 DSPPVecCheck(V2,V3);
485 for (ii=0; ii<n/4; ++ii){
486 i=ii*4;
487 v3[i]=v1[i]*v2[i];
488 v3[i+1]=v1[i+1]*v2[i+1];
489 v3[i+2]=v1[i+2]*v2[i+2];
490 v3[i+3]=v1[i+3]*v2[i+3];
491 }
492 for (i=4*(n/4); i<n; i++){
493 v3[i]=v1[i]*v2[i];
494 }
495
496 return 0;
497}
498
499#undef __FUNCT__
500#define __FUNCT__ "DSDPVecPointwiseDivide"
501int DSDPVecPointwiseDivide( DSDPVec V1, DSDPVec V2, DSDPVec V3){
502
503 int ii,i,n=V1.dim;
504 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
505
506 DSPPVecCheck(V1,V3);
507 DSPPVecCheck(V2,V3);
508 for (ii=0; ii<n/4; ++ii){
509 i=ii*4;
510 v3[i]=v1[i]/v2[i]; v3[i+1]=v1[i+1]/v2[i+1]; v3[i+2]=v1[i+2]/v2[i+2]; v3[i+3]=v1[i+3]/v2[i+3];
511 }
512 for (i=4*(n/4); i<n; i++){
513 v3[i]=v1[i]/v2[i];
514 }
515
516 return 0;
517}
518
519#undef __FUNCT__
520#define __FUNCT__ "DSDPVecShift"
521int DSDPVecShift(double alpha, DSDPVec V){
522 int i,n=V.dim;
523 double *v=V.val;
524 for (i=0; i<n; i++){
525 v[i]+= alpha;
526 }
527
528 return 0;
529}
530
531#undef __FUNCT__
532#define __FUNCT__ "DSDPVecSemiNorm"
533int DSDPVecSemiNorm(DSDPVec V, double *ans){
534
535 int i;
536 double dtmp=0.0;
537
538 for (i=0; i<V.dim; i++){
539 dtmp=min(V.val[i],dtmp);
540 }
541 *ans = fabs(dtmp);
542 if (*ans!=*ans) return 1;
543
544 return 0;
545}
546
547#undef __FUNCT__
548#define __FUNCT__ "DSDPVecReciprocal"
549int DSDPVecReciprocal(DSDPVec V){
550
551 int i,n=V.dim;
552 double *val=V.val;
553
554 for (i=0; i<n; i++){
555 val[i]= 1.0/val[i];
556 }
557
558 return 0;
559}
560
561#undef __FUNCT__
562#define __FUNCT__ "DSDPVecReciprocalSqrt"
563int DSDPVecReciprocalSqrt(DSDPVec V){
564
565 int i,n=V.dim;
566 double *val=V.val;
567
568 for (i=0; i<n; i++){
569 val[i]= sqrt(1.0/val[i]);
570 }
571
572 return 0;
573}
574
575#undef __FUNCT__
576#define __FUNCT__ "DSDPVecAbsoluteValue"
577int DSDPVecAbsoluteValue(DSDPVec V){
578
579 int i,n=V.dim;
580 double *val=V.val;
581 for (i=0; i<n; i++){
582 val[i]=fabs(val[i]);
583 }
584 return 0;
585}
586
587
DSDP uses BLAS and LAPACK for many of its operations.
Error handling, printing, and profiling.
Vector operations used by the solver.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition dsdpvec.h:25