DSDP
|
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 #undef __FUNCT__ 00042 #define __FUNCT__ "DSDPVecCreateWArray" 00043 int DSDPVecCreateWArray(DSDPVec *V, double* vv, int n){ 00044 V->dim=n; 00045 if (n>0){ 00046 V->val=vv; 00047 } else { 00048 V->val=NULL; 00049 } 00050 return 0; 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 int DSDPVecGetSize(DSDPVec V, int *n){ 00069 *n=V.dim; 00070 return 0; 00071 } 00072 00073 int DSDPVecGetArray(DSDPVec V, double **dptr){ 00074 *dptr=V.val; 00075 return 0; 00076 } 00077 00078 int DSDPVecRestoreArray(DSDPVec V, double **dptr){ 00079 *dptr=0; 00080 return 0; 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 #undef __FUNCT__ 00234 #define __FUNCT__ "DSDPVecNorm22" 00235 int DSDPVecNorm22( DSDPVec v, double *vnorm){ 00236 int i,n=v.dim; 00237 double *val=v.val; 00238 00239 *vnorm = 0.0; 00240 for (i=0; i<n; i++){ 00241 *vnorm += val[i]*val[i]; 00242 } 00243 return 0; 00244 } 00245 #undef __FUNCT__ 00246 #define __FUNCT__ "DSDPVecNorm2" 00247 int DSDPVecNorm2( DSDPVec v, double *vnorm){ 00248 int info; 00249 info=DSDPVecNorm22(v,vnorm); if (info) return 1; 00250 if (*vnorm!=*vnorm) return 1; 00251 *vnorm = sqrt(*vnorm); 00252 return 0; 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 #undef __FUNCT__ 00273 #define __FUNCT__ "DSDPVecScale" 00274 int DSDPVecScale(double alpha, DSDPVec x){ 00275 int i,ii,n; 00276 double *xx=x.val; 00277 n=x.dim; 00278 00279 for (ii=0; ii<n/4; ++ii){ 00280 i=ii*4; 00281 xx[i]*= alpha; 00282 xx[i+1]*= alpha; 00283 xx[i+2]*= alpha; 00284 xx[i+3]*= alpha; 00285 } 00286 for (i=4*(n/4); i<n; ++i){ 00287 xx[i]*= alpha; 00288 } 00289 return 0; 00290 } 00291 00292 #undef __FUNCT__ 00293 #define __FUNCT__ "DSDPVecAXPY" 00294 int DSDPVecAXPY(double alpha, DSDPVec x, DSDPVec y){ 00295 00296 int i,ii,n=x.dim; 00297 double *yy=y.val,*xx=x.val; 00298 00299 DSPPVecCheck(x,y); 00300 00301 for (ii=0; ii<n/4; ++ii){ 00302 i=ii*4; 00303 yy[i] += (alpha)*xx[i]; 00304 yy[i+1] += (alpha)*xx[i+1]; 00305 yy[i+2] += (alpha)*xx[i+2]; 00306 yy[i+3] += (alpha)*xx[i+3]; 00307 } 00308 for (i=4*(n/4); i<n; ++i){ 00309 yy[i] += (alpha)*xx[i]; 00310 } 00311 00312 return 0; 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 #undef __FUNCT__ 00424 #define __FUNCT__ "DSDPVecDot" 00425 int DSDPVecDot(DSDPVec V1, DSDPVec V2, double *ans){ 00426 00427 int i,ii,m=V1.dim; 00428 double *v1=V1.val,*v2=V2.val; 00429 00430 DSPPVecCheck(V1,V2); 00431 *ans=0.0; 00432 for (ii=0; ii<m/4; ++ii){ 00433 i=ii*4; 00434 *ans += v1[i]*v2[i] + v1[i+1]*v2[i+1] + v1[i+2]*v2[i+2] + v1[i+3]*v2[i+3] ; 00435 } 00436 for (i=4*(m/4); i<m; ++i){ 00437 *ans += v1[i]*v2[i]; 00438 } 00439 if (*ans!=*ans) return 1; 00440 00441 return 0; 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