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 
17 static int nvecs=0;
18 #undef __FUNCT__
19 #define __FUNCT__ "DSDPVecCreateSeq"
20 int 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"
28 int 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"
43 int 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"
55 int 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 /*
68 int DSDPVecGetSize(DSDPVec V, int *n){
69  *n=V.dim;
70  return 0;
71 }
72 
73 int DSDPVecGetArray(DSDPVec V, double **dptr){
74  *dptr=V.val;
75  return 0;
76 }
77 
78 int DSDPVecRestoreArray(DSDPVec V, double **dptr){
79  *dptr=0;
80  return 0;
81 }
82 */
83 
84 
85 #undef __FUNCT__
86 #define __FUNCT__ "DSDPVecView"
87 int 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"
98 int 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"
108 int DSDPVecSetValue(DSDPVec V,int row,double value){
109  V.val[row]=value;
110  return 0;
111 }
112 
113 #undef __FUNCT__
114 #define __FUNCT__ "DSDPVecZero"
115 int 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"
125 int 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"
137 int 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"
146 int 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"
160 int 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"
172 int 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"
183 int 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"
194 int 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"
204 int 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"
215 int 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"
224 int 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"
235 int 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"
247 int 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"
257 int 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"
274 int 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"
294 int 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"
317 int 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"
340 int 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"
350 int 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"
372 int 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"
394 int 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"
403 int 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"
425 int 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"
447 int 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"
463 int 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"
478 int 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"
501 int 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"
521 int 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"
533 int 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"
549 int 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"
563 int 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"
577 int 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 
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
Error handling, printing, and profiling.
Vector operations used by the solver.
DSDP uses BLAS and LAPACK for many of its operations.