sdpsss.c

Go to the documentation of this file.
00001 #include "dsdpsdp.h"
00002 #include "dsdpsys.h"
00008 int DSDPCreateS(DSDPBlockData*,char,int,DSDPVec,DSDPVMat,SDPConeVec, SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
00009 
00010 static int DSDPCreateDS(DSDPBlockData*, DSDPVMat,int*,int,int,int,int*,int*,DSDPDSMat*);
00011 static int DSDPCreateDS2(DSDPBlockData*, DSDPVMat,int*,int,int,int,int*,int*,DSDPDSMat*);
00012 
00013 static int DSDPCreateS1(DSDPBlockData*,int,DSDPVec,DSDPVMat,SDPConeVec, SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
00014 static int DSDPCreateS2(DSDPBlockData*,int,DSDPVec,DSDPVMat,SDPConeVec, SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
00015 
00016 extern int DSDPXMatPCreate(int,struct DSDPVMat_Ops**,void**);
00017 extern int DSDPXMatPCreateWithData(int,double[],int,struct DSDPVMat_Ops**,void**);
00018 extern int DSDPXMatUCreate(int,struct DSDPVMat_Ops**,void**);
00019 extern int DSDPXMatUCreateWithData(int,double[],int,struct DSDPVMat_Ops**,void**);
00020 
00021 extern int DSDPLAPACKSUDualMatCreate2(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
00022 extern int DSDPLAPACKSUDualMatCreate2P(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
00023 extern int DSDPLAPACKPUDualMatCreate2(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
00024 extern int DSDPDiagDualMatCreateP(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
00025 extern int DSDPDiagDualMatCreateU(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
00026 extern int DSDPDenseDualMatCreate(int, char,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
00027 extern int DSDPSparseDualMatCreate(int, int*, int*, int,char,int*,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
00028 
00029 extern int DSDPSparseMatCreatePattern2P(int,int[],int[],int,struct DSDPDSMat_Ops**,void**);
00030 extern int DSDPSparseMatCreatePattern2U(int,int[],int[],int,struct DSDPDSMat_Ops**,void**);
00031 
00032 extern int DSDPCreateDiagDSMatP(int,struct DSDPDSMat_Ops**,void**);
00033 extern int DSDPCreateDiagDSMatU(int,struct DSDPDSMat_Ops**,void**);
00034 
00035 extern int DSDPCreateDSMatWithArray(int,double[],int,struct DSDPDSMat_Ops**,void**);
00036 extern int DSDPCreateDSMatWithArray2(int,double[],int, struct DSDPDSMat_Ops**,void**);
00037 
00038 extern int DSDPCreateURDSMat(int,struct DSDPDSMat_Ops**,void**);
00039 
00040 /*
00041 #undef __FUNCT__
00042 #define __FUNCT__ "DSDPSetDualMatrix"
00043 int DSDPSetDualMatrix(SDPCone sdpcone,int (*createdualmatrix)(DSDPBlockData*,DSDPVec,DSDPVMat,SDPConeVec,SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*),void*ctx){
00044   DSDPFunctionBegin;
00045   sdpcone->createdualmatrix=createdualmatrix;
00046   sdpcone->dualmatctx=ctx;
00047   DSDPFunctionReturn(0);
00048 }
00049 */
00050 #undef __FUNCT__
00051 #define __FUNCT__ "CountNonzeros"
00052 static int CountNonzeros(DSDPBlockData *ADATA,int m, int rnnz[], int innz[],int n,int *nnz1, int *nnz2)
00053 {
00054   int i,j,info,totalnnz1=0,totalnnz2=0;
00055   
00056   for (i=0;i<n;i++){
00057     memset(rnnz,0,n*sizeof(int));
00058     /* SParsity pattern for DS only requires the constraint matrices A and residual */
00059     for (j=0;j<m;j++) innz[j]=1;innz[0]=0;
00060     info=DSDPBlockDataRowSparsity(ADATA,i,innz,rnnz,n);DSDPCHKERR(info);
00061     for (j=0; j<i; j++){ if (rnnz[j]>0) totalnnz1++;}
00062     /* Adjacency pattern for S also requires the objective matrix C */
00063     for (j=0;j<m;j++) innz[j]=0;innz[0]=1;
00064     info=DSDPBlockDataRowSparsity(ADATA,i,innz,rnnz,n);DSDPCHKERR(info);
00065     for (j=0; j<i; j++){ if (rnnz[j]>0) totalnnz2++;}
00066   }
00067   *nnz1=totalnnz1;
00068   *nnz2=totalnnz2;
00069   return 0;
00070 }
00071 
00072 #undef __FUNCT__
00073 #define __FUNCT__ "CreateS1b"
00074 static int CreateS1b(DSDPBlockData *ADATA, int innz[], int m, int n, int tnnz[], int rnnz[], int snnz[])
00075 {
00076   int i,j,info;
00077   if (ADATA->nnzmats<=0){return 0;};
00078 
00079   memset(rnnz,0,n*sizeof(int));
00080   for (i=0;i<m;i++) innz[i]=1;
00081   innz[0]=0;
00082 
00083   /* Check matrices A for nonzero entries. */
00084   for (i=0;i<n;i++){
00085     memset(tnnz,0,n*sizeof(int));
00086     info=DSDPBlockDataRowSparsity(ADATA,i,innz,tnnz,n);DSDPCHKERR(info);
00087     for (j=0; j<=i; j++){
00088       if (tnnz[j]>0){ *snnz=j; snnz++; rnnz[i]++;}
00089     }
00090   }
00091   return 0;
00092 }
00093 
00094 
00095 #undef __FUNCT__
00096 #define __FUNCT__ "DSDPCreateDS"
00097 int DSDPCreateDS(DSDPBlockData *ADATA, DSDPVMat T, int iworkm[], int m, int n, int nnz, int rnnz[], int tnnz[], DSDPDSMat *B){
00098   int i,n1,*cols,allnnz,info;
00099   double *ss;
00100   void *dsmat;
00101   struct DSDPDSMat_Ops* dsops;
00102   DSDPFunctionBegin;
00103   DSDPLogInfo(0,19,"DS Matrix has %d nonzeros of %d\n",nnz,n*(n-1)/2);
00104   if (nnz==0){
00105     info=DSDPCreateDiagDSMatP(n,&dsops,&dsmat); DSDPCHKERR(info);
00106     DSDPLogInfo(0,19,"Using Diagonal Delta S matrix\n");
00107   } else if (2*nnz +n < n*n/5){
00108     info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
00109     cols=(int*)ss;
00110     info = CreateS1b(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
00111     for (allnnz=0,i=0;i<n;i++){allnnz+=rnnz[i];}
00112     info = DSDPSparseMatCreatePattern2P(n,rnnz,cols,allnnz,&dsops,&dsmat); DSDPCHKERR(info);
00113     info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
00114     DSDPLogInfo(0,19,"Using Sparse Delta S matrix\n");
00115   } else {
00116     info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
00117     info=DSDPCreateDSMatWithArray(n,ss,n1,&dsops,&dsmat); DSDPCHKERR(info);
00118     info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
00119     DSDPLogInfo(0,19,"Using Full Delta S matrix\n");
00120   }
00121   info=DSDPDSMatSetData(B,dsops,dsmat); DSDPCHKERR(info);
00122   DSDPFunctionReturn(0);
00123 }
00124 
00125 
00126 #undef __FUNCT__
00127 #define __FUNCT__ "CreateS1c"
00128 static int CreateS1c(DSDPBlockData *ADATA, int innz[], int m, int n, int tnnz[], int rnnz[], int snnz[])
00129 {
00130   int i,j,info;
00131   memset(rnnz,0,n*sizeof(int));
00132   for (i=0;i<m;i++) innz[i]=1;
00133   /* Check matrix C and A for nonzero entries. */
00134   for (i=0;i<n;i++){
00135     memset(tnnz,0,n*sizeof(int));
00136     info=DSDPBlockDataRowSparsity(ADATA,i,innz,tnnz,n);DSDPCHKERR(info);
00137     for (j=i+1; j<n; j++){
00138       if (tnnz[j]>0){ *snnz=j; snnz++; rnnz[i]++;}
00139     } 
00140   }
00141   return 0;
00142 }
00143 
00144 static int dsdpuselapack=1;
00145 #undef __FUNCT__  
00146 #define __FUNCT__ "SDPConeUseLAPACKForDualMatrix"
00147 int SDPConeUseLAPACKForDualMatrix(SDPCone sdpcone,int flag){ 
00148   DSDPFunctionBegin;
00149   dsdpuselapack = flag;
00150   DSDPFunctionReturn(0);
00151 }
00152 
00153 
00154 #undef __FUNCT__
00155 #define __FUNCT__ "DSDPCreateS"
00156 static int DSDPCreateS1(DSDPBlockData *ADATA, int trank, DSDPVec WY, DSDPVMat T, SDPConeVec W1, SDPConeVec W2, DSDPDualMat *S, DSDPDualMat *SS, DSDPDSMat *DS, void*ctx){
00157 
00158   int nnz,n,n1,*cols,*rnnz,*tnnz,*iworkm,m,info;
00159   int dsnnz,snnz,sfnnz;
00160   double *pss;
00161   void *smat1,*smat2;
00162   struct DSDPDualMat_Ops *sops1,*sops2;
00163 
00164   DSDPFunctionBegin;
00165   info=DSDPVecGetSize(WY,&m);DSDPCHKERR(info); 
00166   info=DSDPVecGetArray(WY,&pss);DSDPCHKERR(info); 
00167   iworkm=(int*)pss;
00168 
00169   info=SDPConeVecGetSize(W1,&n);DSDPCHKERR(info); 
00170   info=SDPConeVecGetArray(W1,&pss);DSDPCHKERR(info); 
00171   tnnz=(int*)pss;
00172   info=SDPConeVecGetArray(W2,&pss);DSDPCHKERR(info); 
00173   rnnz=(int*)pss;
00174 
00175   DSDPLogInfo(0,19,"Compute Sparsity\n");
00176   info = CountNonzeros(ADATA, m, rnnz, iworkm, n, &dsnnz,&snnz); DSDPCHKERR(info);
00177   nnz=snnz;
00178   /* printf("Nonzeros: %d %d of %d\n",dsnnz,snnz,n*(n-1)/2); */
00179   /* TT and DS could share memory */
00180   info=DSDPCreateDS(ADATA,T,iworkm,m,n,dsnnz,rnnz,tnnz,DS);DSDPCHKERR(info);
00181 
00182   if (nnz==0){
00183     info=DSDPDiagDualMatCreateP(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00184     DSDPLogInfo(0,19,"Using Diagonal S matrix\n");
00185   } else if (2*snnz+n+2<n*n/8){
00186     info=DSDPVMatGetArray(T,&pss,&n1); DSDPCHKERR(info);
00187     cols=(int*)pss;
00188     info = CreateS1c(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
00189     info=DSDPSparseDualMatCreate(n,rnnz,cols,trank,'P',&sfnnz,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00190     info=DSDPVMatRestoreArray(T,&pss,&n1); DSDPCHKERR(info);
00191     /*   printf("NNZ: %d %d of %d\n",2*snnz+n,sfnnz*2+n,n*n); */
00192     DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Sparse S matrix\n",nnz,n*(n-1)/2);
00193     DSDPLogInfo(0,19,"Total rank of block: %d, n= %d\n",trank,n);
00194   } else if (n>20 && dsdpuselapack){
00195     info=DSDPLAPACKSUDualMatCreate2P(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00196     DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Full Dense LAPACK S matrix\n",nnz,n*(n-1)/2);
00197   } else if (dsdpuselapack){
00198     info=DSDPLAPACKPUDualMatCreate2(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00199     DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Packed Dense LAPACK S matrix\n",nnz,n*(n-1)/2);
00200   } else {
00201     info=DSDPDenseDualMatCreate(n,'P',&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00202     DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Dense S matrix\n",nnz,n*(n-1)/2);
00203   }
00204   info=DSDPDualMatSetData(S,sops1,smat1); DSDPCHKERR(info);
00205   info=DSDPDualMatSetData(SS,sops2,smat2); DSDPCHKERR(info);
00206 
00207   DSDPFunctionReturn(0);
00208 }
00209 
00210 
00211 
00212 #undef __FUNCT__
00213 #define __FUNCT__ "DSDPCreateDS2"
00214 int DSDPCreateDS2(DSDPBlockData *ADATA, DSDPVMat T, int iworkm[], int m, int n, int nnz, int rnnz[], int tnnz[], DSDPDSMat *B){
00215   int i,n1,*cols,allnnz,info;
00216   double *ss;
00217   void *dsmat;
00218   struct DSDPDSMat_Ops* dsops;
00219   DSDPFunctionBegin;
00220   DSDPLogInfo(0,19,"DS Matrix has %d nonzeros of %d\n",nnz,n*(n-1)/2);
00221   if (nnz==0){
00222     info=DSDPCreateDiagDSMatU(n,&dsops,&dsmat); DSDPCHKERR(info);
00223     DSDPLogInfo(0,19,"Using Diagonal Delta S matrix\n");
00224   } else if (2*nnz +n < n*n/4){
00225     info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
00226     cols=(int*)ss;
00227     info = CreateS1b(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
00228     for (allnnz=0,i=0;i<n;i++){allnnz+=rnnz[i];}
00229     info = DSDPSparseMatCreatePattern2U(n,rnnz,cols,allnnz,&dsops,&dsmat); DSDPCHKERR(info);
00230     info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
00231     DSDPLogInfo(0,19,"Using Sparse Delta S matrix\n");
00232   } else {
00233     info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
00234     info=DSDPCreateDSMatWithArray2(n,ss,n1,&dsops,&dsmat); DSDPCHKERR(info);
00235     info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
00236     DSDPLogInfo(0,19,"Using Full Delta S matrix\n");
00237   }
00238   info=DSDPDSMatSetData(B,dsops,dsmat); DSDPCHKERR(info);
00239   DSDPFunctionReturn(0);
00240 }
00241 
00242 #undef __FUNCT__
00243 #define __FUNCT__ "DSDPCreateS2"
00244 static int DSDPCreateS2(DSDPBlockData *ADATA, int trank, DSDPVec WY, DSDPVMat T, SDPConeVec W1, SDPConeVec W2, DSDPDualMat *S, DSDPDualMat *SS, DSDPDSMat *DS, void*ctx){
00245 
00246   int nnz,n,n1,*cols,*rnnz,*tnnz,*iworkm,m,info;
00247   int dsnnz,snnz,sfnnz;
00248   double *pss;
00249   void *smat1,*smat2;
00250   struct DSDPDualMat_Ops *sops1,*sops2;
00251   DSDPFunctionBegin;
00252 
00253   info=DSDPVecGetSize(WY,&m);DSDPCHKERR(info); 
00254   info=DSDPVecGetArray(WY,&pss);DSDPCHKERR(info); 
00255   iworkm=(int*)pss;
00256 
00257   info=SDPConeVecGetSize(W1,&n);DSDPCHKERR(info); 
00258   info=SDPConeVecGetArray(W1,&pss);DSDPCHKERR(info); 
00259   tnnz=(int*)pss;
00260   info=SDPConeVecGetArray(W2,&pss);DSDPCHKERR(info); 
00261   rnnz=(int*)pss;
00262 
00263   DSDPLogInfo(0,19,"Compute Sparsity\n");
00264   info = CountNonzeros(ADATA, m, rnnz, iworkm, n, &dsnnz,&snnz); DSDPCHKERR(info);
00265   nnz=snnz;
00266   /* printf("Nonzeros: %d %d of %d\n",dsnnz,snnz,n*(n-1)/2); */
00267   /* TT and DS could share memory */
00268   info=DSDPCreateDS2(ADATA,T,iworkm,m,n,dsnnz,rnnz,tnnz,DS);DSDPCHKERR(info);
00269 
00270   if (nnz==0){
00271     info=DSDPDiagDualMatCreateU(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00272     DSDPLogInfo(0,19,"Using Diagonal S matrix\n");
00273   } else if (2*snnz+n+2<n*n/10){
00274     info=DSDPVMatGetArray(T,&pss,&n1); DSDPCHKERR(info);
00275     cols=(int*)pss;
00276     info = CreateS1c(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
00277     info=DSDPSparseDualMatCreate(n,rnnz,cols,trank,'U',&sfnnz,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00278     info=DSDPVMatRestoreArray(T,&pss,&n1); DSDPCHKERR(info);
00279     DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Sparse S matrix\n",nnz,n*(n-1)/2);
00280   } else if (dsdpuselapack){
00281     info=DSDPLAPACKSUDualMatCreate2(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00282     DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Full Dense LAPACK S matrix\n",nnz,n*(n-1)/2);
00283   } else {
00284     info=DSDPDenseDualMatCreate(n,'U',&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
00285     DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Packed Dense S matrix\n",nnz,n*(n-1)/2);
00286   }
00287   info=DSDPDualMatSetData(S,sops1,smat1); DSDPCHKERR(info);
00288   info=DSDPDualMatSetData(SS,sops2,smat2); DSDPCHKERR(info);
00289 
00290   DSDPFunctionReturn(0);
00291 }
00292 
00293 
00294 int DSDPCreateS(DSDPBlockData*,char,int,DSDPVec,DSDPVMat,SDPConeVec,SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
00295 
00296 #undef __FUNCT__
00297 #define __FUNCT__ "DSDPCreateS"
00298 
00314 int DSDPCreateS(DSDPBlockData *ADATA, char UPLQ, int trank, DSDPVec WY, DSDPVMat T, SDPConeVec W1, SDPConeVec W2, DSDPDualMat *S, DSDPDualMat *SS, DSDPDSMat *DS, void*ctx){
00315   int info;
00316   DSDPFunctionBegin;
00317   switch (UPLQ){
00318   case 'P':
00319     info=DSDPCreateS1(ADATA,trank,WY,T,W1,W2,S,SS,DS,ctx);DSDPCHKERR(info); 
00320     break;
00321   case 'U':
00322     info=DSDPCreateS2(ADATA,trank,WY,T,W1,W2,S,SS,DS,ctx);DSDPCHKERR(info); 
00323     break;
00324   }
00325   DSDPFunctionReturn(0);
00326 }
00327 
00328 
00329 
00330 #undef __FUNCT__
00331 #define __FUNCT__ "DSDPCreateS"
00332 int DSDPUseDefaultDualMatrix(SDPCone sdpcone){
00333   DSDPFunctionBegin;
00334   /*
00335   int info;
00336   info=DSDPSetDualMatrix(sdpcone,DSDPCreateS2,0); DSDPCHKERR(info);
00337   */
00338   DSDPFunctionReturn(0);
00339 }
00340 
00341 #undef __FUNCT__
00342 #define __FUNCT__ "DSDPMakeVMat"
00343 
00351 int DSDPMakeVMat(char UPLQ, int n, DSDPVMat *X){
00352   int info;
00353   void *xmat;
00354   struct DSDPVMat_Ops* xops;
00355   DSDPFunctionBegin;
00356   switch (UPLQ){
00357   case 'P':
00358     info=DSDPXMatPCreate(n,&xops,&xmat);DSDPCHKERR(info); 
00359     break;
00360   case 'U':
00361     info=DSDPXMatUCreate(n,&xops,&xmat);DSDPCHKERR(info); 
00362     break;
00363   }
00364   info=DSDPVMatSetData(X,xops,xmat); DSDPCHKERR(info);
00365   DSDPFunctionReturn(0);
00366 }
00367 
00368 #undef __FUNCT__
00369 #define __FUNCT__ "DSDPMakeVMatWithArray"
00370 
00381 int DSDPMakeVMatWithArray(char UPLQ, double xx[], int nnz, int n, DSDPVMat *X){
00382   int info;
00383   void *xmat;
00384   struct DSDPVMat_Ops* xops;
00385   DSDPFunctionBegin;
00386   switch (UPLQ){
00387   case 'P':
00388     info=DSDPXMatPCreateWithData(n,xx,nnz,&xops,&xmat);DSDPCHKERR(info); 
00389     break;
00390   case 'U':
00391     info=DSDPXMatUCreateWithData(n,xx,nnz,&xops,&xmat);DSDPCHKERR(info); 
00392     break;
00393   }
00394   info=DSDPVMatSetData(X,xops,xmat); DSDPCHKERR(info);
00395   DSDPFunctionReturn(0);
00396 }

Generated on Sun Mar 23 07:30:49 2008 for DSDP by  doxygen 1.5.5