dsdplp.c

Go to the documentation of this file.
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 

Generated on Thu May 22 09:42:47 2008 for DSDP by  doxygen 1.5.5