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
00053
00054
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
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
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
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
00568
00569
00570
00571
00572
00573
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
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
00938
00939
00940
00941
00942
00943
00944
00945
00946
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
00972
00973
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
00994
00995
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
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
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
01040
01041
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