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
00042
00043
00044
00045
00046
00047
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
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
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
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
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
00179
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
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
00267
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
00336
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 }