DSDP
|
00001 00002 #include "dsdpcone_impl.h" 00003 #include "dsdpsys.h" 00004 #include "dsdp5.h" 00005 00006 #include <stdio.h> 00007 #include <stdlib.h> 00008 #include <math.h> 00009 #include <string.h> 00014 typedef struct { 00015 int nrow; 00016 int ncol; 00017 int owndata; 00018 00019 const double *an; 00020 const int *col; 00021 const int *nnz; 00022 int *nzrows; 00023 int nnzrows; 00024 00025 } smatx; 00026 00031 struct LPCone_C{ 00032 smatx *A,*AT; 00033 DSDPVec C; 00034 DSDPVec PS,DS,X; 00035 double sscale; 00036 double r; 00037 double muscale; 00038 DSDPVec Y,WY,WY2,WX,WX2; 00039 double *xout; 00040 int n,m; 00041 }; 00042 00043 00044 static int CreateSpRowMatWdata(int,int,const double[], const int[], const int[],smatx **); 00045 static int SpRowMatMult(smatx*,const double[], int , double[], int); 00046 static int SpRowMatMultTrans(smatx *, const double[],int, double[],int); 00047 static int SpRowMatGetRowVector(smatx*, int, double*,int); 00048 static int SpRowMatGetScaledRowVector(smatx*, int, const double[], double*, int); 00049 static int SpRowMatDestroy(smatx*); 00050 static int SpRowMatView(smatx*); 00051 /* 00052 static int SpRowMatGetSize(smatx *, int *, int *); 00053 static int SpRowMatZero(smatx*); 00054 static int SpRowMatAddRowMultiple(smatx*, int, double, double[], int); 00055 */ 00056 static int SpRowIMultAdd(smatx*,int*,int,int *,int); 00057 static int SpRowMatRowNnz(smatx*, int, int*,int); 00058 static int SpRowMatNorm2(smatx*, int, double*); 00059 00060 00061 #undef __FUNCT__ 00062 #define __FUNCT__ "LPConeSetUp" 00063 static int LPConeSetup(void *dcone,DSDPVec y){ 00064 int m,info; 00065 LPCone lpcone=(LPCone)dcone; 00066 DSDPFunctionBegin; 00067 if (lpcone->n<1) return 0; 00068 m=lpcone->m; 00069 info=DSDPVecCreateSeq(m+2,&lpcone->WY);DSDPCHKERR(info); 00070 info=DSDPVecDuplicate(lpcone->WY,&lpcone->WY2);DSDPCHKERR(info); 00071 info=DSDPVecDuplicate(lpcone->WY,&lpcone->Y);DSDPCHKERR(info); 00072 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX);DSDPCHKERR(info); 00073 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX2);DSDPCHKERR(info); 00074 info=DSDPVecDuplicate(lpcone->C,&lpcone->PS);DSDPCHKERR(info); 00075 info=DSDPVecDuplicate(lpcone->C,&lpcone->DS);DSDPCHKERR(info); 00076 info=DSDPVecDuplicate(lpcone->C,&lpcone->X);DSDPCHKERR(info); 00077 DSDPFunctionReturn(0); 00078 } 00079 00080 #undef __FUNCT__ 00081 #define __FUNCT__ "LPConeSetUp2" 00082 static int LPConeSetup2(void *dcone, DSDPVec Y, DSDPSchurMat M){ 00083 LPCone lpcone=(LPCone)dcone; 00084 DSDPFunctionBegin; 00085 DSDPLogInfo(0,19,"Setup LP Cone of dimension: %d\n",lpcone->n); 00086 DSDPFunctionReturn(0); 00087 } 00088 00089 00090 #undef __FUNCT__ 00091 #define __FUNCT__ "LPConeDestroy" 00092 static int LPConeDestroy(void *dcone){ 00093 int info; 00094 LPCone lpcone=(LPCone)dcone; 00095 DSDPFunctionBegin; 00096 if (lpcone->n<1) return 0; 00097 info=DSDPVecDestroy(&lpcone->DS);DSDPCHKERR(info); 00098 info=DSDPVecDestroy(&lpcone->PS);DSDPCHKERR(info); 00099 info=DSDPVecDestroy(&lpcone->C);DSDPCHKERR(info); 00100 info=DSDPVecDestroy(&lpcone->X);DSDPCHKERR(info); 00101 info=SpRowMatDestroy(lpcone->A);DSDPCHKERR(info); 00102 info=DSDPVecDestroy(&lpcone->WX2);DSDPCHKERR(info); 00103 info=DSDPVecDestroy(&lpcone->WY2);DSDPCHKERR(info); 00104 info=DSDPVecDestroy(&lpcone->WY);DSDPCHKERR(info); 00105 info=DSDPVecDestroy(&lpcone->Y);DSDPCHKERR(info); 00106 info=DSDPVecDestroy(&lpcone->WX);DSDPCHKERR(info); 00107 DSDPFREE(&lpcone,&info);DSDPCHKERR(info); 00108 DSDPFunctionReturn(0); 00109 } 00110 00111 #undef __FUNCT__ 00112 #define __FUNCT__ "LPConeSize" 00113 static int LPConeSize(void *dcone, double *n){ 00114 LPCone lpcone=(LPCone)dcone; 00115 DSDPFunctionBegin; 00116 *n=lpcone->muscale*lpcone->n; 00117 DSDPFunctionReturn(0); 00118 } 00119 00120 00121 00122 #undef __FUNCT__ 00123 #define __FUNCT__ "LPComputeAX" 00124 static int LPComputeAX( LPCone lpcone,DSDPVec X, DSDPVec Y){ 00125 int info,m=lpcone->m,n=lpcone->n; 00126 double *x,*y,ppobj; 00127 smatx *A=lpcone->A; 00128 DSDPFunctionBegin; 00129 if (lpcone->n<1) return 0; 00130 info=DSDPVecGetSize(X,&n);DSDPCHKERR(info); 00131 info=DSDPVecDot(lpcone->C,X,&ppobj);DSDPCHKERR(info); 00132 info=DSDPVecSetC(Y,ppobj); 00133 info=DSDPVecSum(X,&ppobj);DSDPCHKERR(info); 00134 info=DSDPVecSetR(Y,ppobj*lpcone->r);DSDPCHKERR(info); 00135 info=DSDPVecGetArray(Y,&y);DSDPCHKERR(info); 00136 info=DSDPVecGetArray(X,&x);DSDPCHKERR(info); 00137 info=SpRowMatMult(A,x,n,y+1,m); 00138 info=DSDPVecRestoreArray(X,&x);DSDPCHKERR(info); 00139 info=DSDPVecRestoreArray(Y,&y);DSDPCHKERR(info); 00140 DSDPFunctionReturn(0); 00141 } 00142 00143 #undef __FUNCT__ 00144 #define __FUNCT__ "LPComputeATY" 00145 static int LPComputeATY(LPCone lpcone,DSDPVec Y, DSDPVec S){ 00146 int info,m=lpcone->m,n=lpcone->n; 00147 double cc,r,*s,*y; 00148 DSDPVec C=lpcone->C; 00149 smatx *A=lpcone->A; 00150 DSDPFunctionBegin; 00151 if (lpcone->n<1) return 0; 00152 info=DSDPVecGetC(Y,&cc);DSDPCHKERR(info); 00153 info=DSDPVecGetR(Y,&r);DSDPCHKERR(info); 00154 info=DSDPVecGetSize(S,&n);DSDPCHKERR(info); 00155 info=DSDPVecGetArray(S,&s);DSDPCHKERR(info); 00156 info=DSDPVecGetArray(Y,&y);DSDPCHKERR(info); 00157 info=SpRowMatMultTrans(A,y+1,m,s,n); DSDPCHKERR(info); 00158 info=DSDPVecRestoreArray(S,&s);DSDPCHKERR(info); 00159 info=DSDPVecRestoreArray(Y,&y);DSDPCHKERR(info); 00160 info=DSDPVecAXPY(cc,C,S);DSDPCHKERR(info); 00161 info=DSDPVecShift(r*lpcone->r,S);DSDPCHKERR(info); 00162 info=DSDPVecScale(-1.0,S);DSDPCHKERR(info); 00163 DSDPFunctionReturn(0); 00164 } 00165 #undef __FUNCT__ 00166 #define __FUNCT__ "LPConeHessian" 00167 static int LPConeHessian(void* dcone, double mu, DSDPSchurMat M, 00168 DSDPVec vrhs1, DSDPVec vrhs2){ 00169 int info,i,m,n,ncols; 00170 double r=1.0,*wx,*wx2; 00171 LPCone lpcone=(LPCone)dcone; 00172 DSDPVec WX=lpcone->WX,WX2=lpcone->WX2,WY=lpcone->WY,WY2=lpcone->WY2,S=lpcone->DS; 00173 smatx *A=lpcone->A; 00174 00175 DSDPFunctionBegin; 00176 if (lpcone->n<1) return 0; 00177 mu*=lpcone->muscale; 00178 info=DSDPVecGetSize(vrhs1,&m);DSDPCHKERR(info); 00179 info=DSDPVecGetSize(WX,&n);DSDPCHKERR(info); 00180 info=DSDPVecSet(mu,WX2);DSDPCHKERR(info); 00181 info=DSDPVecPointwiseDivide(WX2,S,WX2);DSDPCHKERR(info); 00182 info=DSDPVecPointwiseDivide(WX2,S,WX2);DSDPCHKERR(info); 00183 for (i=0;i<m;i++){ 00184 00185 info=DSDPSchurMatRowColumnScaling(M,i,WY2,&ncols); DSDPCHKERR(info); 00186 if (ncols==0) continue; 00187 00188 if (i==0){ 00189 info=DSDPVecPointwiseMult(lpcone->C,WX2,WX); DSDPCHKERR(info); 00190 } else if (i==m-1){ 00191 info=DSDPVecScaleCopy(WX2,r,WX); DSDPCHKERR(info); 00192 } else { 00193 info=DSDPVecGetArray(WX,&wx); 00194 info=DSDPVecGetArray(WX2,&wx2);DSDPCHKERR(info); 00195 info=SpRowMatGetScaledRowVector(A,i-1,wx2,wx,n); 00196 info=DSDPVecRestoreArray(WX,&wx); 00197 info=DSDPVecRestoreArray(WX2,&wx2); 00198 } 00199 00200 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info); 00201 00202 info=DSDPVecPointwiseMult(WY2,WY,WY);DSDPCHKERR(info); 00203 00204 info=DSDPSchurMatAddRow(M,i,1.0,WY);DSDPCHKERR(info); 00205 } 00206 00207 /* Compute AS^{-1} */ 00208 info=DSDPVecSet(mu,WX);DSDPCHKERR(info); 00209 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info); 00210 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info); 00211 00212 info=DSDPSchurMatDiagonalScaling(M, WY2);DSDPCHKERR(info); 00213 info=DSDPVecPointwiseMult(WY2,WY,WY);DSDPCHKERR(info); 00214 info=DSDPVecAXPY(1.0,WY,vrhs2);DSDPCHKERR(info); 00215 00216 DSDPFunctionReturn(0); 00217 } 00218 00219 #undef __FUNCT__ 00220 #define __FUNCT__ "LPConeHessian" 00221 static int LPConeRHS(void* dcone, double mu, DSDPVec vrow, 00222 DSDPVec vrhs1, DSDPVec vrhs2){ 00223 int info; 00224 LPCone lpcone=(LPCone)dcone; 00225 DSDPVec WX=lpcone->WX,WY=lpcone->WY,S=lpcone->DS; 00226 00227 DSDPFunctionBegin; 00228 if (lpcone->n<1) return 0; 00229 mu*=lpcone->muscale; 00230 00231 /* Compute AS^{-1} */ 00232 info=DSDPVecSet(mu,WX);DSDPCHKERR(info); 00233 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info); 00234 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info); 00235 00236 info=DSDPVecPointwiseMult(vrow,WY,WY);DSDPCHKERR(info); 00237 info=DSDPVecAXPY(1.0,WY,vrhs2);DSDPCHKERR(info); 00238 00239 DSDPFunctionReturn(0); 00240 } 00241 00242 #undef __FUNCT__ 00243 #define __FUNCT__ "LPConeMultiply" 00244 static int LPConeMultiply(void* dcone, double mu, DSDPVec vrow, DSDPVec vin, DSDPVec vout){ 00245 int info; 00246 LPCone lpcone=(LPCone)dcone; 00247 DSDPVec WX=lpcone->WX,S=lpcone->DS,WY=lpcone->WY; 00248 DSDPFunctionBegin; 00249 if (lpcone->n<1) return 0; 00250 mu*=lpcone->muscale; 00251 info=LPComputeATY(lpcone,vin,WX);DSDPCHKERR(info); 00252 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info); 00253 info=DSDPVecScale(mu,WX);DSDPCHKERR(info); 00254 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info); 00255 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info); 00256 info=DSDPVecPointwiseMult(WY,vrow,WY);DSDPCHKERR(info); 00257 info=DSDPVecAXPY(1.0,WY,vout);DSDPCHKERR(info); 00258 DSDPFunctionReturn(0); 00259 } 00260 00261 #undef __FUNCT__ 00262 #define __FUNCT__ "LPConeSetX" 00263 static int LPConeSetX(void* dcone,double mu, DSDPVec Y,DSDPVec DY){ 00264 DSDPFunctionBegin; 00265 DSDPFunctionReturn(0); 00266 } 00267 00268 #undef __FUNCT__ 00269 #define __FUNCT__ "LPConeX" 00270 static int LPConeX(void* dcone,double mu, DSDPVec Y,DSDPVec DY, 00271 DSDPVec AX , double *tracexs){ 00272 int info; 00273 double dtracexs; 00274 LPCone lpcone=(LPCone)dcone; 00275 DSDPVec S=lpcone->PS,WX=lpcone->WX,X=lpcone->X,DS=lpcone->DS,WY=lpcone->WY; 00276 double *xx,*xout=lpcone->xout; 00277 DSDPFunctionBegin; 00278 00279 if (lpcone->n<1) return 0; 00280 mu*=lpcone->muscale; 00281 info=LPComputeATY(lpcone,Y,S);DSDPCHKERR(info); 00282 00283 info=DSDPVecSet(1.0,WX); 00284 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info); 00285 00286 info=LPComputeATY(lpcone,DY,DS);DSDPCHKERR(info); 00287 info=DSDPVecPointwiseMult(WX,DS,X);DSDPCHKERR(info); 00288 00289 info=DSDPVecScale(-mu,WX);DSDPCHKERR(info); 00290 info=DSDPVecPointwiseMult(WX,X,X);DSDPCHKERR(info); 00291 info=DSDPVecAXPY(-1.0,WX,X);DSDPCHKERR(info); 00292 info=DSDPVecGetArray(X,&xx);DSDPCHKERR(info); 00293 for (info=0;info<lpcone->n;info++){ 00294 if (xx[info]<0) xx[info]=0; 00295 } 00296 info=DSDPVecRestoreArray(X,&xx);DSDPCHKERR(info); 00297 info=LPComputeAX(lpcone,X,WY);DSDPCHKERR(info); 00298 info=DSDPVecAXPY(1.0,WY,AX);DSDPCHKERR(info); 00299 info=DSDPVecDot(S,X,&dtracexs);DSDPCHKERR(info); 00300 *tracexs+=dtracexs; 00301 info=DSDPVecGetArray(X,&xx);DSDPCHKERR(info); 00302 if (xout){ 00303 for (info=0;info<lpcone->n;info++){ 00304 if (xout){ xout[info]=xx[info]; } 00305 } 00306 } 00307 info=DSDPVecRestoreArray(X,&xx);DSDPCHKERR(info); 00308 DSDPFunctionReturn(0); 00309 } 00310 00311 00312 #undef __FUNCT__ 00313 #define __FUNCT__ "LPConeS" 00314 static int LPConeS(void* dcone, DSDPVec Y, DSDPDualFactorMatrix flag, 00315 DSDPTruth *psdefinite){ 00316 int i,n,info; 00317 double *s; 00318 LPCone lpcone=(LPCone)dcone; 00319 DSDPVec S; 00320 00321 DSDPFunctionBegin; 00322 00323 if (lpcone->n<1) return 0; 00324 if (flag==DUAL_FACTOR){ 00325 S=lpcone->DS; 00326 } else { 00327 S=lpcone->PS; 00328 } 00329 00330 info=DSDPVecCopy(Y,lpcone->Y);DSDPCHKERR(info); 00331 info=LPComputeATY(lpcone,Y,S);DSDPCHKERR(info); 00332 info=DSDPVecGetC(Y,&lpcone->sscale); 00333 info=DSDPVecGetSize(S,&n);DSDPCHKERR(info); 00334 info=DSDPVecGetArray(S,&s);DSDPCHKERR(info); 00335 *psdefinite=DSDP_TRUE; 00336 for (i=0;i<n;i++){ if (s[i]<=0) *psdefinite=DSDP_FALSE;} 00337 info=DSDPVecRestoreArray(S,&s);DSDPCHKERR(info); 00338 00339 DSDPFunctionReturn(0); 00340 } 00341 #undef __FUNCT__ 00342 #define __FUNCT__ "LPConeInvertS" 00343 static int LPConeInvertS(void* dcone){ 00344 DSDPFunctionBegin; 00345 DSDPFunctionReturn(0); 00346 } 00347 00348 #undef __FUNCT__ 00349 #define __FUNCT__ "LPConeComputeMaxStepLength" 00350 static int LPConeComputeMaxStepLength(void* dcone, DSDPVec DY, DSDPDualFactorMatrix flag, double *maxsteplength){ 00351 int i,n,info; 00352 double *s,*ds,mstep=1.0e200; 00353 LPCone lpcone=(LPCone)dcone; 00354 DSDPVec S,DS=lpcone->WX; 00355 DSDPFunctionBegin; 00356 if (lpcone->n<1) return 0; 00357 if (flag==DUAL_FACTOR){ 00358 S=lpcone->DS; 00359 } else { 00360 S=lpcone->PS; 00361 } 00362 00363 info=LPComputeATY(lpcone,DY,DS);DSDPCHKERR(info); 00364 00365 info=DSDPVecGetSize(DS,&n);DSDPCHKERR(info); 00366 info=DSDPVecGetArray(S,&s);DSDPCHKERR(info); 00367 info=DSDPVecGetArray(DS,&ds);DSDPCHKERR(info); 00368 for (i=0;i<n;i++) if (ds[i]<0){mstep=DSDPMin(mstep,-s[i]/ds[i]);} 00369 info=DSDPVecRestoreArray(S,&s);DSDPCHKERR(info); 00370 info=DSDPVecRestoreArray(DS,&ds);DSDPCHKERR(info); 00371 00372 *maxsteplength=mstep; 00373 00374 DSDPFunctionReturn(0); 00375 } 00376 00377 00378 #undef __FUNCT__ 00379 #define __FUNCT__ "LPConePotential" 00380 static int LPConePotential(void* dcone, double *logobj, double *logdet){ 00381 int i,n,info; 00382 double *s,mu,sumlog=0; 00383 LPCone lpcone=(LPCone)dcone; 00384 DSDPFunctionBegin; 00385 if (lpcone->n<1) return 0; 00386 mu=lpcone->muscale; 00387 info=DSDPVecGetArray(lpcone->DS,&s);DSDPCHKERR(info); 00388 info=DSDPVecGetSize(lpcone->DS,&n);DSDPCHKERR(info); 00389 for (i=0;i<n;i++){ 00390 sumlog+= mu*log(s[i]); 00391 } 00392 info=DSDPVecRestoreArray(lpcone->DS,&s);DSDPCHKERR(info); 00393 *logdet=sumlog; 00394 *logobj=0; 00395 DSDPFunctionReturn(0); 00396 } 00397 00398 #undef __FUNCT__ 00399 #define __FUNCT__ "LPConeSparsity" 00400 static int LPConeSparsity(void *dcone,int row, int *tnnz, int rnnz[], int m){ 00401 int info,*wi,n; 00402 double *wd; 00403 LPCone lpcone=(LPCone)dcone; 00404 DSDPVec W=lpcone->WX; 00405 DSDPFunctionBegin; 00406 if (lpcone->n<1) return 0; 00407 if (row==0) return 0; 00408 if (row==m-1){ 00409 return 0; 00410 } 00411 info=DSDPVecGetSize(W,&n);DSDPCHKERR(info); 00412 info=DSDPVecGetArray(W,&wd);DSDPCHKERR(info); 00413 wi=(int*)wd; 00414 info=SpRowMatRowNnz(lpcone->A,row-1,wi,n);DSDPCHKERR(info); 00415 info=SpRowIMultAdd(lpcone->A,wi,n,rnnz+1,m-2);DSDPCHKERR(info); 00416 info=DSDPVecRestoreArray(W,&wd);DSDPCHKERR(info); 00417 DSDPFunctionReturn(0); 00418 } 00419 00420 00421 #undef __FUNCT__ 00422 #define __FUNCT__ "LPConeMonitor" 00423 static int LPConeMonitor( void *dcone, int tag){ 00424 DSDPFunctionBegin; 00425 DSDPFunctionReturn(0); 00426 } 00427 00428 #undef __FUNCT__ 00429 #define __FUNCT__ "LPANorm2" 00430 static int LPANorm2( void *dcone, DSDPVec ANorm){ 00431 int i,info; 00432 double dd; 00433 LPCone lpcone=(LPCone)dcone; 00434 DSDPFunctionBegin; 00435 if (lpcone->n<1) return 0; 00436 info=DSDPVecNorm22(lpcone->C,&dd);DSDPCHKERR(info); 00437 info=DSDPVecAddC(ANorm,dd);DSDPCHKERR(info); 00438 for (i=0;i<lpcone->m;i++){ 00439 info=SpRowMatNorm2(lpcone->A,i,&dd);DSDPCHKERR(info); 00440 info=DSDPVecAddElement(ANorm,i+1,dd);DSDPCHKERR(info); 00441 } 00442 info=DSDPVecAddR(ANorm,1.0);DSDPCHKERR(info); 00443 DSDPFunctionReturn(0); 00444 } 00445 00446 00447 static struct DSDPCone_Ops kops; 00448 static const char *lpconename="LP Cone"; 00449 00450 #undef __FUNCT__ 00451 #define __FUNCT__ "LPConeOperationsInitialize" 00452 static int LPConeOperationsInitialize(struct DSDPCone_Ops* coneops){ 00453 int info; 00454 if (coneops==NULL) return 0; 00455 info=DSDPConeOpsInitialize(coneops); DSDPCHKERR(info); 00456 coneops->conehessian=LPConeHessian; 00457 coneops->conerhs=LPConeRHS; 00458 coneops->conesetup=LPConeSetup; 00459 coneops->conesetup2=LPConeSetup2; 00460 coneops->conedestroy=LPConeDestroy; 00461 coneops->conecomputes=LPConeS; 00462 coneops->coneinverts=LPConeInvertS; 00463 coneops->conesetxmaker=LPConeSetX; 00464 coneops->conecomputex=LPConeX; 00465 coneops->conemaxsteplength=LPConeComputeMaxStepLength; 00466 coneops->conelogpotential=LPConePotential; 00467 coneops->conesize=LPConeSize; 00468 coneops->conesparsity=LPConeSparsity; 00469 coneops->conehmultiplyadd=LPConeMultiply; 00470 coneops->conemonitor=LPConeMonitor; 00471 coneops->coneanorm2=LPANorm2; 00472 coneops->id=2; 00473 coneops->name=lpconename; 00474 return 0; 00475 } 00476 00477 #undef __FUNCT__ 00478 #define __FUNCT__ "DSDPAddLP" 00479 int DSDPAddLP(DSDP dsdp,LPCone lpcone){ 00480 int info; 00481 DSDPFunctionBegin; 00482 info=LPConeOperationsInitialize(&kops); DSDPCHKERR(info); 00483 info=DSDPAddCone(dsdp,&kops,(void*)lpcone); DSDPCHKERR(info); 00484 DSDPFunctionReturn(0); 00485 } 00486 00487 #undef __FUNCT__ 00488 #define __FUNCT__ "DSDPCreateLPCone" 00489 00509 int DSDPCreateLPCone(DSDP dsdp, LPCone *dspcone){ 00510 int m,info; 00511 struct LPCone_C *lpcone; 00512 DSDPFunctionBegin; 00513 DSDPCALLOC1(&lpcone,struct LPCone_C,&info);DSDPCHKERR(info); 00514 *dspcone=lpcone; 00515 /* 00516 info=DSDPAddLP(dsdp,lpcone);DSDPCHKERR(info); 00517 */ 00518 info=LPConeOperationsInitialize(&kops); DSDPCHKERR(info); 00519 info=DSDPAddCone(dsdp,&kops,(void*)lpcone); DSDPCHKERR(info); 00520 info=DSDPGetNumberOfVariables(dsdp,&m);DSDPCHKERR(info); 00521 lpcone->m=m; 00522 lpcone->muscale=1.0; 00523 lpcone->n=0; 00524 lpcone->xout=0; 00525 lpcone->r=1.0; 00526 info=DSDPVecCreateSeq(0,&lpcone->C);DSDPCHKERR(info); 00527 info=DSDPVecCreateSeq(0,&lpcone->WY);DSDPCHKERR(info); 00528 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX);DSDPCHKERR(info); 00529 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX2);DSDPCHKERR(info); 00530 info=DSDPVecDuplicate(lpcone->C,&lpcone->PS);DSDPCHKERR(info); 00531 info=DSDPVecDuplicate(lpcone->C,&lpcone->DS);DSDPCHKERR(info); 00532 info=DSDPVecDuplicate(lpcone->C,&lpcone->X);DSDPCHKERR(info); 00533 DSDPFunctionReturn(0); 00534 } 00535 00536 00537 #undef __FUNCT__ 00538 #define __FUNCT__ "LPConeGetXArray" 00539 00556 int LPConeGetXArray(LPCone lpcone,double *x[], int *n){ 00557 int info; 00558 DSDPFunctionBegin; 00559 info=DSDPVecGetArray(lpcone->X,x);DSDPCHKERR(info); 00560 info=DSDPVecGetSize(lpcone->X,n);DSDPCHKERR(info); 00561 DSDPFunctionReturn(0); 00562 } 00563 00564 #undef __FUNCT__ 00565 #define __FUNCT__ "LPConeGetSArray" 00566 /* 00567 \fn int LPConeGetSArray(LPCone lpcone, double *s[], int *n); 00568 \brief Get the array used to store the s variables 00569 \ingroup LPRoutines 00570 \param lpcone LP Cone 00571 \param s array of variables 00572 \param n the dimension of the cone and length of the array 00573 \sa DSDPCreateLPCone() 00574 */ 00575 int LPConeGetSArray(LPCone lpcone,double *s[], int *n){ 00576 int info; 00577 DSDPFunctionBegin; 00578 info=DSDPVecGetArray(lpcone->DS,s);DSDPCHKERR(info); 00579 info=DSDPVecGetSize(lpcone->DS,n);DSDPCHKERR(info); 00580 DSDPFunctionReturn(0); 00581 } 00582 00583 #undef __FUNCT__ 00584 #define __FUNCT__ "LPConeCopyS" 00585 00595 int LPConeCopyS(LPCone lpcone,double s[], int n){ 00596 int i,info; 00597 double *ss,sscale=lpcone->sscale; 00598 DSDPTruth psdefinite; 00599 DSDPFunctionBegin; 00600 info=LPConeS((void*)lpcone,lpcone->Y,DUAL_FACTOR ,&psdefinite);DSDPCHKERR(info); 00601 info=DSDPVecGetArray(lpcone->DS,&ss);DSDPCHKERR(info); 00602 for (i=0;i<n;i++) s[i]=ss[i]/fabs(sscale); 00603 DSDPFunctionReturn(0); 00604 } 00605 00606 #undef __FUNCT__ 00607 #define __FUNCT__ "LPConeGetDimension" 00608 00616 int LPConeGetDimension(LPCone lpcone,int *n){ 00617 DSDPFunctionBegin; 00618 *n=(int)(lpcone->n*lpcone->muscale); 00619 DSDPFunctionReturn(0); 00620 } 00621 00622 00623 #undef __FUNCT__ 00624 #define __FUNCT__ "LPConeScaleBarrier" 00625 int LPConeScaleBarrier(LPCone lpcone,double muscale){ 00626 DSDPFunctionBegin; 00627 if (muscale>0){ 00628 lpcone->muscale=muscale; 00629 } 00630 DSDPFunctionReturn(0); 00631 } 00632 00633 #undef __FUNCT__ 00634 #define __FUNCT__ "LPConeSetData" 00635 00666 int LPConeSetData(LPCone lpcone,int n, const int ik[],const int cols[],const double vals[]){ 00667 int info,i,spot,m=lpcone->m; 00668 DSDPVec C; 00669 DSDPFunctionBegin; 00670 lpcone->n=n; 00671 info=DSDPVecCreateSeq(n,&C);DSDPCHKERR(info); 00672 lpcone->C=C; 00673 info=DSDPVecZero(C);DSDPCHKERR(info); 00674 lpcone->muscale=1.0; 00675 if (n<100) lpcone->muscale=1.0; 00676 if (n<10) lpcone->muscale=1.0; 00677 for (i=ik[0];i<ik[1];i++){ 00678 info=DSDPVecSetElement(C,cols[i],vals[i]); 00679 } 00680 spot=ik[0]; 00681 info=CreateSpRowMatWdata(m,n,vals+spot,cols+spot,ik+1,&lpcone->A);DSDPCHKERR(info); 00682 DSDPFunctionReturn(0); 00683 } 00684 00685 #undef __FUNCT__ 00686 #define __FUNCT__ "LPConeSetData2" 00687 00717 int LPConeSetData2(LPCone lpcone,int n, const int ik[],const int cols[],const double vals[]){ 00718 int info,i,spot,m=lpcone->m; 00719 DSDPVec C; 00720 DSDPFunctionBegin; 00721 lpcone->n=n; 00722 info=DSDPVecCreateSeq(n,&C);DSDPCHKERR(info); 00723 lpcone->C=C; 00724 info=DSDPVecZero(C);DSDPCHKERR(info); 00725 lpcone->muscale=1.0; 00726 if (n<100) lpcone->muscale=1.0; 00727 if (n<10) lpcone->muscale=1.0; 00728 for (i=ik[m];i<ik[m+1];i++){ 00729 info=DSDPVecSetElement(C,cols[i],vals[i]); 00730 } 00731 spot=ik[0]; 00732 info=CreateSpRowMatWdata(m,n,vals+spot,cols+spot,ik,&lpcone->A);DSDPCHKERR(info); 00733 DSDPFunctionReturn(0); 00734 } 00735 00736 #undef __FUNCT__ 00737 #define __FUNCT__ "LPConeView2" 00738 00744 int LPConeView2(LPCone lpcone){ 00745 int info; 00746 DSDPFunctionBegin; 00747 printf("LPCone Constraint Matrix\n"); 00748 info=SpRowMatView(lpcone->A);DSDPCHKERR(info); 00749 printf("LPCone Objective C vector\n"); 00750 info=DSDPVecView(lpcone->C);DSDPCHKERR(info); 00751 DSDPFunctionReturn(0); 00752 } 00753 00754 00755 00756 #undef __FUNCT__ 00757 #define __FUNCT__ "LPConeGetConstraint" 00758 int LPConeGetConstraint(LPCone lpcone,int column,DSDPVec W){ 00759 int n,info; 00760 double *w; 00761 DSDPFunctionBegin; 00762 if (column==0){ 00763 info=DSDPVecCopy(lpcone->C,W);DSDPCHKERR(info); 00764 } else { 00765 info=DSDPVecGetSize(W,&n);DSDPCHKERR(info); 00766 info=DSDPVecGetArray(W,&w);DSDPCHKERR(info); 00767 info=SpRowMatGetRowVector(lpcone->A,column-1,w,n);info=0;DSDPCHKERR(info); 00768 info=DSDPVecRestoreArray(W,&w);DSDPCHKERR(info); 00769 } 00770 DSDPFunctionReturn(0); 00771 } 00772 00773 #undef __FUNCT__ 00774 #define __FUNCT__ "LPConeGetData" 00775 00783 int LPConeGetData(LPCone lpcone,int vari,double vv[], int n){ 00784 int info; 00785 DSDPVec W; 00786 DSDPFunctionBegin; 00787 info=DSDPVecCreateWArray(&W,vv,n);DSDPCHKERR(info); 00788 info=LPConeGetConstraint(lpcone,vari,W);DSDPCHKERR(info); 00789 DSDPFunctionReturn(0); 00790 } 00791 00792 #undef __FUNCT__ 00793 #define __FUNCT__ "LPConeSetXVec" 00794 int LPConeSetXVec(LPCone lpcone,double *xout, int n){ 00795 DSDPFunctionBegin; 00796 if (n==lpcone->n) lpcone->xout=xout; 00797 DSDPFunctionReturn(0); 00798 } 00799 00800 00801 static int vsdot(const int*,const double *,int,const double *, int, double *); 00802 static int checkvsparse(smatx *); 00803 00804 00805 static int CreateSpRowMatWdata(int m,int n,const double vals[], const int cols[], const int ik[], 00806 smatx **A){ 00807 00808 smatx *V; 00809 00810 V=(smatx*) malloc(1*sizeof(smatx)); 00811 if (V==NULL) return 1; 00812 V->nrow=m; 00813 V->ncol=n; 00814 V->owndata=0; 00815 V->an=vals; V->col=cols; V->nnz=ik; 00816 00817 *A=V; 00818 checkvsparse(V); 00819 return 0; 00820 } 00821 00822 static int vsdot(const int ja[], const double an[], int nn0, const double vv[], int n, double *vdot){ 00823 00824 int i; 00825 double tmp=0.0; 00826 00827 for (i=0; i<nn0; i++){ 00828 /* if (ja[i]<n) tmp = tmp + an[i] * vv[ja[i]]; */ 00829 tmp += an[i] * vv[ja[i]]; 00830 } 00831 *vdot=tmp; 00832 return 0; 00833 } 00834 00835 static int checkvsparse(smatx *A){ 00836 int i,k=0,m=A->nrow,tnnz=0; 00837 const int *nnz=A->nnz; 00838 00839 for (i=0;i<m;++i){ 00840 if (nnz[i+1]-nnz[i]>0){ 00841 tnnz++; 00842 } 00843 } 00844 if (tnnz < m/2){ 00845 A->nzrows =(int*)malloc((tnnz)*sizeof(int)); 00846 A->nnzrows=tnnz; 00847 for (i=0;i<m;++i){ 00848 if (nnz[i+1]-nnz[i]>0){ 00849 A->nzrows[k]=i; 00850 k++; 00851 } 00852 } 00853 } else { 00854 A->nzrows = 0; 00855 A->nnzrows = m; 00856 } 00857 return 0; 00858 } 00859 00860 #undef __FUNCT__ 00861 #define __FUNCT__ "SpRowMatMult" 00862 static int SpRowMatMult(smatx* A, const double x[], int n, double y[], int m){ 00863 00864 int i,k1,k2,nrow=A->nrow; 00865 const int *nnz=A->nnz,*col=A->col; 00866 const double *an=A->an; 00867 00868 if (A->ncol != n) return 1; 00869 if (A->nrow != m) return 2; 00870 if (x==0 && n>0) return 3; 00871 if (y==0 && m>0) return 3; 00872 00873 if (m>0){ 00874 memset((void*)y,0,m*sizeof(double)); 00875 for (i=0; i<nrow; i++){ 00876 k1=*(nnz+i); k2=*(nnz+i+1); 00877 vsdot(col+k1,an+k1,k2-k1,x,n,y+i); 00878 } 00879 } 00880 return 0; 00881 } 00882 00883 #undef __FUNCT__ 00884 #define __FUNCT__ "SpRowMatMultTrans" 00885 static int SpRowMatMultTrans(smatx * A, const double x[], int m, double y[], int n){ 00886 00887 int i,j,k1,k2,nrow=A->nrow; 00888 const int *col=A->col,*nnz=A->nnz; 00889 const double *an=A->an; 00890 if (A->ncol != n) return 1; 00891 if (A->nrow != m) return 2; 00892 if (x==0 && m>0) return 3; 00893 if (y==0 && n>0) return 3; 00894 00895 memset((void*)y,0,n*sizeof(double)); 00896 for (i=0; i<nrow; i++){ 00897 k1=nnz[i]; k2=nnz[i+1]; 00898 for (j=k1; j<k2; j++){ 00899 y[col[j]] += an[j]*x[i]; 00900 } 00901 } 00902 00903 return 0; 00904 } 00905 00906 00907 static int SpRowMatRowNnz(smatx *A, int row, int* wi,int n){ 00908 int i,k1,k2; 00909 const int *col=A->col; 00910 DSDPFunctionBegin; 00911 memset((void*)wi,0,n*sizeof(double)); 00912 k1=A->nnz[row]; 00913 k2=A->nnz[row+1]; 00914 for (i=k1; i<k2; i++){ 00915 wi[col[i]]=1; 00916 } 00917 DSDPFunctionReturn(0); 00918 } 00919 00920 static int SpRowIMultAdd(smatx *A,int *wi,int n,int *rnnz,int m){ 00921 00922 int i,j,k1,k2,nrow=A->nrow; 00923 const int *nnz=A->nnz,*col=A->col; 00924 DSDPFunctionBegin; 00925 for (i=0; i<nrow; i++){ 00926 k1=nnz[i]; 00927 k2=nnz[i+1]; 00928 for (j=k1; j<k2; j++){ 00929 if (wi[col[j]]){ 00930 rnnz[i]++; 00931 } 00932 } 00933 } 00934 DSDPFunctionReturn(0); 00935 } 00936 /* 00937 static int SpRowMatAddRowMultiple(smatx* A, int nrow, double ytmp, double row[], int n){ 00938 int k; 00939 int *col=A->col, *nnz=A->nnz; 00940 double *an=A->an; 00941 00942 for (k=nnz[nrow]; k<nnz[nrow+1]; k++){ 00943 row[col[k]] += ytmp * an[k]; 00944 } 00945 00946 return 0; 00947 } 00948 */ 00949 static int SpRowMatNorm2(smatx* A, int nrow, double *norm22){ 00950 int k; 00951 const int *nnz=A->nnz; 00952 double tt=0; 00953 const double *an=A->an; 00954 00955 for (k=nnz[nrow]; k<nnz[nrow+1]; k++){ 00956 tt+=an[k]*an[k]; 00957 } 00958 *norm22=tt; 00959 return 0; 00960 } 00961 00962 00963 #undef __FUNCT__ 00964 #define __FUNCT__ "SpRowMatGetRowVector" 00965 static int SpRowMatGetRowVector(smatx* M, int row, double r[], int m){ 00966 00967 int i,k1,k2; 00968 const int *col=M->col; 00969 const double *an=M->an; 00970 /* 00971 if (M->ncol != m) return 1; 00972 if (row<0 || row>= M->nrow) return 2; 00973 if (r==0) return 3; 00974 */ 00975 memset((void*)r,0,m*sizeof(double)); 00976 k1=M->nnz[row]; 00977 k2=M->nnz[row+1]; 00978 00979 for (i=k1; i<k2; i++){ 00980 r[col[i]]=an[i]; 00981 } 00982 00983 return 0; 00984 } 00985 #undef __FUNCT__ 00986 #define __FUNCT__ "SpRowMatGetScaledRowVector" 00987 static int SpRowMatGetScaledRowVector(smatx* M, int row, const double ss[], double r[], int m){ 00988 00989 int i,k1,k2; 00990 const int *col=M->col; 00991 const double *an=M->an; 00992 /* 00993 if (M->ncol != m) return 1; 00994 if (row<0 || row>= M->nrow) return 2; 00995 if (r==0) return 3; 00996 */ 00997 memset((void*)r,0,m*sizeof(double)); 00998 k1=M->nnz[row]; 00999 k2=M->nnz[row+1]; 01000 01001 for (i=k1; i<k2; i++){ 01002 r[col[i]]=ss[col[i]]*an[i]; 01003 } 01004 01005 return 0; 01006 } 01007 01008 01009 /* 01010 #undef __FUNCT__ 01011 #define __FUNCT__ "SpRowMatZero" 01012 static int SpRowMatZero(smatx* M){ 01013 01014 int nnz=M->nnz[M->nrow]; 01015 memset(M->an,0,(nnz)*sizeof(double)); 01016 01017 return 0; 01018 } 01019 01020 01021 01022 #undef __FUNCT__ 01023 #define __FUNCT__ "SpRowGetSize" 01024 static int SpRowMatGetSize(smatx *M, int *m, int *n){ 01025 *m=M->nrow; 01026 *n=M->ncol; 01027 01028 return 0; 01029 } 01030 */ 01031 #undef __FUNCT__ 01032 #define __FUNCT__ "SpRowMatDestroy" 01033 static int SpRowMatDestroy(smatx* A){ 01034 01035 if (A->owndata){ 01036 printf("Can't free array"); 01037 return 1; 01038 /* 01039 if (A->an) free(A->an); 01040 if (A->col) free(A->col); 01041 if (A->nnz) free(A->nnz); 01042 */ 01043 } 01044 if (A->nzrows) free(A->nzrows); 01045 free(A); 01046 01047 return 0; 01048 } 01049 01050 #undef __FUNCT__ 01051 #define __FUNCT__ "SpRowMatView" 01052 static int SpRowMatView(smatx* M){ 01053 01054 int i,j,k1,k2; 01055 01056 for (i=0; i<M->nrow; i++){ 01057 k1=M->nnz[i]; k2=M->nnz[i+1]; 01058 01059 if (k2-k1 >0){ 01060 printf("Row %d, (Variable y%d) : ",i,i+1); 01061 for (j=k1; j<k2; j++) 01062 printf(" %4.2e x%d + ",M->an[j],M->col[j]); 01063 printf("= dobj%d \n",i+1); 01064 } 01065 } 01066 01067 return 0; 01068 } 01069 01070 #undef __FUNCT__ 01071 #define __FUNCT__ "LPConeView" 01072 01078 int LPConeView(LPCone lpcone){ 01079 01080 smatx* A=lpcone->A; 01081 int i,j,jj,info; 01082 const int *row=A->col,*nnz=A->nnz; 01083 int n=A->ncol,m=A->nrow; 01084 const double *an=A->an; 01085 double cc; 01086 DSDPVec C=lpcone->C; 01087 DSDPFunctionBegin; 01088 printf("LPCone Constraint Matrix\n"); 01089 printf("Number y variables 1 through %d\n",m); 01090 for (i=0; i<n; i++){ 01091 printf("Inequality %d: ",i); 01092 for (j=0;j<m;j++){ 01093 for (jj=nnz[j];jj<nnz[j+1];jj++){ 01094 if (row[jj]==i){ 01095 printf("%4.2e y%d + ",an[jj],j+1); 01096 } 01097 } 01098 } 01099 info=DSDPVecGetElement(C,i,&cc);DSDPCHKERR(info); 01100 printf(" <= %4.2e\n",cc); 01101 } 01102 DSDPFunctionReturn(0); 01103 } 01104 01105