00001 #include "dsdpsys.h"
00002 #include "dsdpvec.h"
00003 #include "dsdplapack.h"
00008 #if !defined (min)
00009 #define min(a,b) ((a <= b)? (a) : (b))
00010 #endif
00011 #if !defined (max)
00012 #define max(a,b) ((a >= b)? (a) : (b))
00013 #endif
00014
00015 #define DSPPVecCheck(a,b) {if (a.dim != b.dim) return 1; if (a.dim>0 && (a.val==NULL || b.val==NULL) ) return 2;}
00016
00017 static int nvecs=0;
00018 #undef __FUNCT__
00019 #define __FUNCT__ "DSDPVecCreateSeq"
00020 int DSDPVecCreate(DSDPVec *V){
00021 int info;
00022 info = DSDPVecCreateSeq(0,V);DSDPCHKERR(info);
00023 return 0;
00024 }
00025
00026 #undef __FUNCT__
00027 #define __FUNCT__ "DSDPVecCreateSeq"
00028 int DSDPVecCreateSeq(int n ,DSDPVec *V){
00029 int info;
00030 V->dim=n;
00031 if (n>0){
00032 nvecs++;
00033 DSDPCALLOC2(&(V->val),double,n,&info);DSDPCHKERR(info);
00034 if (V->val==NULL) return 1;
00035 } else {
00036 V->val=NULL;
00037 }
00038 return 0;
00039 }
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 #undef __FUNCT__
00054 #define __FUNCT__ "DSDPVecDestroy"
00055 int DSDPVecDestroy(DSDPVec *V){
00056 int info;
00057 if ((*V).val){
00058 DSDPFREE(&(*V).val,&info);DSDPCHKERR(info);
00059 nvecs--;
00060 }
00061
00062 (*V).dim=0;
00063 (*V).val=0;
00064 return 0;
00065 }
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085 #undef __FUNCT__
00086 #define __FUNCT__ "DSDPVecView"
00087 int DSDPVecView(DSDPVec vec){
00088 int i;
00089 for (i=0; i<vec.dim; i++){
00090 printf("%3.3e ",vec.val[i]);
00091 }
00092 printf("\n");
00093 return 0;
00094 }
00095
00096 #undef __FUNCT__
00097 #define __FUNCT__ "DSDPVecISet"
00098 int DSDPVecISet(int* ival,DSDPVec V){
00099 int i;
00100 for (i=0;i<V.dim;i++){
00101 V.val[i]=ival[i];
00102 }
00103 return 0;
00104 }
00105
00106 #undef __FUNCT__
00107 #define __FUNCT__ "DSDPVecSetValue"
00108 int DSDPVecSetValue(DSDPVec V,int row,double value){
00109 V.val[row]=value;
00110 return 0;
00111 }
00112
00113 #undef __FUNCT__
00114 #define __FUNCT__ "DSDPVecZero"
00115 int DSDPVecZero(DSDPVec V){
00116 int n=V.dim;
00117 double *v=V.val;
00118 memset((void*)v,0,n*sizeof(double));
00119 return 0;
00120 }
00121
00122
00123 #undef __FUNCT__
00124 #define __FUNCT__ "DSDPVecNormalize"
00125 int DSDPVecNormalize(DSDPVec V){
00126 int info;
00127 double vnorm;
00128 info = DSDPVecNorm2(V,&vnorm);DSDPCHKERR(info);
00129 if (vnorm==0){ return 1;}
00130 vnorm=1.0/(vnorm);
00131 info = DSDPVecScale(vnorm,V);DSDPCHKERR(info);
00132 return 0;
00133 }
00134
00135 #undef __FUNCT__
00136 #define __FUNCT__ "DSDPVecSetBasis"
00137 int DSDPVecSetBasis(DSDPVec V,int row){
00138 int info;
00139 info=DSDPVecZero(V);
00140 V.val[row]=1.0;
00141 return 0;
00142 }
00143
00144 #undef __FUNCT__
00145 #define __FUNCT__ "DSDPVecCopy"
00146 int DSDPVecCopy( DSDPVec v1, DSDPVec v2){
00147
00148 int n=v1.dim;
00149 double *val1=v1.val,*val2=v2.val;
00150 DSPPVecCheck(v1,v2);
00151 if (val1!=val2){
00152 memcpy(val2,val1,n*sizeof(double));
00153 }
00154 return 0;
00155 }
00156
00157
00158 #undef __FUNCT__
00159 #define __FUNCT__ "DSDPVecSum"
00160 int DSDPVecSum( DSDPVec v, double *vnorm){
00161 int i,n;
00162 n = v.dim;
00163 *vnorm = 0.0;
00164 for (i=0; i<n; i++){
00165 *vnorm += v.val[i];
00166 }
00167 if (*vnorm!=*vnorm) return 1;
00168 return 0;
00169 }
00170 #undef __FUNCT__
00171 #define __FUNCT__ "DSDPVecNorm1"
00172 int DSDPVecNorm1( DSDPVec v, double *vnorm){
00173 ffinteger N=v.dim,INCX=1;
00174 *vnorm=0;
00175 *vnorm=dasum(&N,v.val,&INCX);
00176 if (*vnorm!=*vnorm) return 1;
00177 return 0;
00178 }
00179
00180
00181 #undef __FUNCT__
00182 #define __FUNCT__ "DSDPVecDot"
00183 int DSDPVecDot(DSDPVec V1, DSDPVec V2, double *ans){
00184 ffinteger ione=1, nn=V1.dim;
00185 double *v1=V1.val,*v2=V2.val;
00186 *ans=ddot(&nn,v1,&ione,v2,&ione);
00187 if (*ans!=*ans) return 1;
00188 return 0;
00189 }
00190
00191
00192 #undef __FUNCT__
00193 #define __FUNCT__ "DSDPVecNorm22"
00194 int DSDPVecNorm22( DSDPVec VV, double *vnorm){
00195 ffinteger ione=1,nn=VV.dim;
00196 double dd,*v=VV.val;
00197 dd=dnrm2(&nn,v,&ione);
00198 *vnorm = dd*dd;
00199 if (*vnorm!=*vnorm) return 1;
00200 return 0;
00201 }
00202 #undef __FUNCT__
00203 #define __FUNCT__ "DSDPVecNorm2"
00204 int DSDPVecNorm2( DSDPVec VV, double *vnorm){
00205 ffinteger ione=1,nn=VV.dim;
00206 double dd,*v=VV.val;
00207 dd=dnrm2(&nn,v,&ione);
00208 *vnorm = dd;
00209 if (*vnorm!=*vnorm) return 1;
00210 return 0;
00211 }
00212
00213 #undef __FUNCT__
00214 #define __FUNCT__ "DSDPVecScale"
00215 int DSDPVecScale(double alpha, DSDPVec VV){
00216 ffinteger ione=1,nn=VV.dim;
00217 double *v=VV.val;
00218 dscal(&nn,&alpha,v,&ione);
00219 return 0;
00220 }
00221
00222 #undef __FUNCT__
00223 #define __FUNCT__ "DSDPVecAXPY"
00224 int DSDPVecAXPY(double alpha, DSDPVec x, DSDPVec y){
00225 ffinteger ione=1,nn=x.dim;
00226 double *yy=y.val,*xx=x.val;
00227 if (alpha==0) return 0;
00228 daxpy(&nn,&alpha,xx,&ione,yy,&ione);
00229 return 0;
00230 }
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 #undef __FUNCT__
00256 #define __FUNCT__ "DSDPVecNormInfinity"
00257 int DSDPVecNormInfinity( DSDPVec v, double *vnorm){
00258
00259 int i,n=v.dim;
00260 double *val=v.val;
00261
00262 *vnorm = 0.0;
00263
00264 for (i=0; i<n; i++){
00265 *vnorm = max(*vnorm,fabs(val[i]));
00266 }
00267 if (*vnorm!=*vnorm) return 1;
00268
00269 return 0;
00270 }
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 #undef __FUNCT__
00316 #define __FUNCT__ "DSDPVecWAXPBY"
00317 int DSDPVecWAXPBY(DSDPVec w, double alpha, DSDPVec x, double beta, DSDPVec y){
00318
00319 int i,ii,n=x.dim;
00320 double *yy=y.val,*xx=x.val,*ww=w.val;
00321 DSPPVecCheck(x,y);
00322 DSPPVecCheck(x,w);
00323
00324 for (ii=0; ii<n/4; ++ii){
00325 i=ii*4;
00326 ww[i] = (alpha)*xx[i] + (beta)*yy[i];
00327 ww[i+1] = (alpha)*xx[i+1] + (beta)*yy[i+1];
00328 ww[i+2] = (alpha)*xx[i+2] + (beta)*yy[i+2];
00329 ww[i+3] = (alpha)*xx[i+3] + (beta)*yy[i+3];
00330 }
00331 for (i=4*(n/4); i<n; ++i){
00332 ww[i] = (alpha)*xx[i] + (beta)*yy[i];
00333 }
00334
00335 return 0;
00336 }
00337
00338 #undef __FUNCT__
00339 #define __FUNCT__ "DSDPVecWAXPY"
00340 int DSDPVecWAXPY(DSDPVec w,double alpha, DSDPVec x, DSDPVec y){
00341 int info;
00342 info=DSDPVecCopy(y,w);
00343 info=DSDPVecAXPY(alpha,x,w);
00344 return 0;
00345 }
00346
00347
00348 #undef __FUNCT__
00349 #define __FUNCT__ "DSDPVecAYPX"
00350 int DSDPVecAYPX(double alpha, DSDPVec x, DSDPVec y){
00351
00352 int i,ii,n=x.dim;
00353 double *yy=y.val,*xx=x.val;
00354
00355 DSPPVecCheck(x,y);
00356 for (ii=0; ii<n/4; ++ii){
00357 i=ii*4;
00358 yy[i] = xx[i]+(alpha)*yy[i];
00359 yy[i+1] = xx[i+1]+(alpha)*yy[i+1];
00360 yy[i+2] = xx[i+2]+(alpha)*yy[i+2];
00361 yy[i+3] = xx[i+3]+(alpha)*yy[i+3];
00362 }
00363 for (i=4*(n/4); i<n; ++i){
00364 yy[i] = xx[i]+(alpha)*yy[i];
00365 }
00366
00367 return 0;
00368 }
00369
00370 #undef __FUNCT__
00371 #define __FUNCT__ "DSDPVecAYPX"
00372 int DSDPVecScaleCopy(DSDPVec x, double alpha, DSDPVec y){
00373
00374 int i,ii,n=x.dim;
00375 double *yy=y.val,*xx=x.val;
00376
00377 DSPPVecCheck(x,y);
00378 for (ii=0; ii<n/4; ++ii){
00379 i=ii*4;
00380 yy[i] = (alpha)*xx[i];
00381 yy[i+1] = (alpha)*xx[i+1];
00382 yy[i+2] = (alpha)*xx[i+2];
00383 yy[i+3] = (alpha)*xx[i+3];
00384 }
00385 for (i=4*(n/4); i<n; ++i){
00386 yy[i] = (alpha)*xx[i];
00387 }
00388
00389 return 0;
00390 }
00391
00392 #undef __FUNCT__
00393 #define __FUNCT__ "DSDPVecDuplicate"
00394 int DSDPVecDuplicate(DSDPVec V1,DSDPVec *V2){
00395 int info,n=V1.dim;
00396 info = DSDPVecCreateSeq(n ,V2);DSDPCHKERR(info);
00397 return 0;
00398 }
00399
00400
00401 #undef __FUNCT__
00402 #define __FUNCT__ "DSDPVecSet"
00403 int DSDPVecSet(double alpha, DSDPVec V){
00404
00405 int i,ii,n=V.dim;
00406 double *val=V.val;
00407
00408 if (alpha==0.0){
00409 memset((void*)val,0,n*sizeof(double));
00410 return 0;
00411 }
00412 for (ii=0; ii<n/4; ++ii){
00413 i=ii*4;
00414 val[i] = val[i+1] = val[i+2] = val[i+3] = alpha;
00415 }
00416 for (i=4*(n/4); i<n; ++i){
00417 val[i]= alpha;
00418 }
00419 return 0;
00420 }
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 #undef __FUNCT__
00446 #define __FUNCT__ "DSDPVecPointwiseMin"
00447 int DSDPVecPointwiseMin( DSDPVec V1, DSDPVec V2, DSDPVec V3){
00448
00449 int i,n=V1.dim;
00450 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
00451
00452 DSPPVecCheck(V1,V3);
00453 DSPPVecCheck(V2,V3);
00454 for (i=0; i<n; ++i){
00455 v3[i]=DSDPMin(v2[i],v1[i]);
00456 }
00457
00458 return 0;
00459 }
00460
00461 #undef __FUNCT__
00462 #define __FUNCT__ "DSDPVecPointwiseMax"
00463 int DSDPVecPointwiseMax( DSDPVec V1, DSDPVec V2, DSDPVec V3){
00464
00465 int i,n=V1.dim;
00466 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
00467
00468 DSPPVecCheck(V1,V3);
00469 DSPPVecCheck(V2,V3);
00470 for (i=0; i<n; ++i){
00471 v3[i]=DSDPMax(v2[i],v1[i]);
00472 }
00473 return 0;
00474 }
00475
00476 #undef __FUNCT__
00477 #define __FUNCT__ "DSDPVecPointwiseMult"
00478 int DSDPVecPointwiseMult( DSDPVec V1, DSDPVec V2, DSDPVec V3){
00479
00480 int ii,i,n=V1.dim;
00481 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
00482
00483 DSPPVecCheck(V1,V3);
00484 DSPPVecCheck(V2,V3);
00485 for (ii=0; ii<n/4; ++ii){
00486 i=ii*4;
00487 v3[i]=v1[i]*v2[i];
00488 v3[i+1]=v1[i+1]*v2[i+1];
00489 v3[i+2]=v1[i+2]*v2[i+2];
00490 v3[i+3]=v1[i+3]*v2[i+3];
00491 }
00492 for (i=4*(n/4); i<n; i++){
00493 v3[i]=v1[i]*v2[i];
00494 }
00495
00496 return 0;
00497 }
00498
00499 #undef __FUNCT__
00500 #define __FUNCT__ "DSDPVecPointwiseDivide"
00501 int DSDPVecPointwiseDivide( DSDPVec V1, DSDPVec V2, DSDPVec V3){
00502
00503 int ii,i,n=V1.dim;
00504 double *v1=V1.val,*v2=V2.val,*v3=V3.val;
00505
00506 DSPPVecCheck(V1,V3);
00507 DSPPVecCheck(V2,V3);
00508 for (ii=0; ii<n/4; ++ii){
00509 i=ii*4;
00510 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];
00511 }
00512 for (i=4*(n/4); i<n; i++){
00513 v3[i]=v1[i]/v2[i];
00514 }
00515
00516 return 0;
00517 }
00518
00519 #undef __FUNCT__
00520 #define __FUNCT__ "DSDPVecShift"
00521 int DSDPVecShift(double alpha, DSDPVec V){
00522 int i,n=V.dim;
00523 double *v=V.val;
00524 for (i=0; i<n; i++){
00525 v[i]+= alpha;
00526 }
00527
00528 return 0;
00529 }
00530
00531 #undef __FUNCT__
00532 #define __FUNCT__ "DSDPVecSemiNorm"
00533 int DSDPVecSemiNorm(DSDPVec V, double *ans){
00534
00535 int i;
00536 double dtmp=0.0;
00537
00538 for (i=0; i<V.dim; i++){
00539 dtmp=min(V.val[i],dtmp);
00540 }
00541 *ans = fabs(dtmp);
00542 if (*ans!=*ans) return 1;
00543
00544 return 0;
00545 }
00546
00547 #undef __FUNCT__
00548 #define __FUNCT__ "DSDPVecReciprocal"
00549 int DSDPVecReciprocal(DSDPVec V){
00550
00551 int i,n=V.dim;
00552 double *val=V.val;
00553
00554 for (i=0; i<n; i++){
00555 val[i]= 1.0/val[i];
00556 }
00557
00558 return 0;
00559 }
00560
00561 #undef __FUNCT__
00562 #define __FUNCT__ "DSDPVecReciprocalSqrt"
00563 int DSDPVecReciprocalSqrt(DSDPVec V){
00564
00565 int i,n=V.dim;
00566 double *val=V.val;
00567
00568 for (i=0; i<n; i++){
00569 val[i]= sqrt(1.0/val[i]);
00570 }
00571
00572 return 0;
00573 }
00574
00575 #undef __FUNCT__
00576 #define __FUNCT__ "DSDPVecAbsoluteValue"
00577 int DSDPVecAbsoluteValue(DSDPVec V){
00578
00579 int i,n=V.dim;
00580 double *val=V.val;
00581 for (i=0; i<n; i++){
00582 val[i]=fabs(val[i]);
00583 }
00584 return 0;
00585 }
00586
00587