DSDP
|
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 }