DSDP
sdpsss.c
Go to the documentation of this file.
1#include "dsdpsdp.h"
2#include "dsdpsys.h"
9
10static int DSDPCreateDS(DSDPBlockData*, DSDPVMat,int*,int,int,int,int*,int*,DSDPDSMat*);
11static int DSDPCreateDS2(DSDPBlockData*, DSDPVMat,int*,int,int,int,int*,int*,DSDPDSMat*);
12
13static int DSDPCreateS1(DSDPBlockData*,int,DSDPVec,DSDPVMat,SDPConeVec, SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
14static int DSDPCreateS2(DSDPBlockData*,int,DSDPVec,DSDPVMat,SDPConeVec, SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
15
16extern int DSDPXMatPCreate(int,struct DSDPVMat_Ops**,void**);
17extern int DSDPXMatPCreateWithData(int,double[],int,struct DSDPVMat_Ops**,void**);
18extern int DSDPXMatUCreate(int,struct DSDPVMat_Ops**,void**);
19extern int DSDPXMatUCreateWithData(int,double[],int,struct DSDPVMat_Ops**,void**);
20
21extern int DSDPLAPACKSUDualMatCreate2(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
22extern int DSDPLAPACKSUDualMatCreate2P(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
23extern int DSDPLAPACKPUDualMatCreate2(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
24extern int DSDPDiagDualMatCreateP(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
25extern int DSDPDiagDualMatCreateU(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
26extern int DSDPDenseDualMatCreate(int, char,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
27extern int DSDPSparseDualMatCreate(int, int*, int*, int,char,int*,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
28
29extern int DSDPSparseMatCreatePattern2P(int,int[],int[],int,struct DSDPDSMat_Ops**,void**);
30extern int DSDPSparseMatCreatePattern2U(int,int[],int[],int,struct DSDPDSMat_Ops**,void**);
31
32extern int DSDPCreateDiagDSMatP(int,struct DSDPDSMat_Ops**,void**);
33extern int DSDPCreateDiagDSMatU(int,struct DSDPDSMat_Ops**,void**);
34
35extern int DSDPCreateDSMatWithArray(int,double[],int,struct DSDPDSMat_Ops**,void**);
36extern int DSDPCreateDSMatWithArray2(int,double[],int, struct DSDPDSMat_Ops**,void**);
37
38extern int DSDPCreateURDSMat(int,struct DSDPDSMat_Ops**,void**);
39
40/*
41#undef __FUNCT__
42#define __FUNCT__ "DSDPSetDualMatrix"
43int DSDPSetDualMatrix(SDPCone sdpcone,int (*createdualmatrix)(DSDPBlockData*,DSDPVec,DSDPVMat,SDPConeVec,SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*),void*ctx){
44 DSDPFunctionBegin;
45 sdpcone->createdualmatrix=createdualmatrix;
46 sdpcone->dualmatctx=ctx;
47 DSDPFunctionReturn(0);
48}
49*/
50#undef __FUNCT__
51#define __FUNCT__ "CountNonzeros"
52static int CountNonzeros(DSDPBlockData *ADATA,int m, int rnnz[], int innz[],int n,int *nnz1, int *nnz2)
53{
54 int i,j,info,totalnnz1=0,totalnnz2=0;
55
56 for (i=0;i<n;i++){
57 memset(rnnz,0,n*sizeof(int));
58 /* SParsity pattern for DS only requires the constraint matrices A and residual */
59 for (j=0;j<m;j++) innz[j]=1;innz[0]=0;
60 info=DSDPBlockDataRowSparsity(ADATA,i,innz,rnnz,n);DSDPCHKERR(info);
61 for (j=0; j<i; j++){ if (rnnz[j]>0) totalnnz1++;}
62 /* Adjacency pattern for S also requires the objective matrix C */
63 for (j=0;j<m;j++) innz[j]=0;innz[0]=1;
64 info=DSDPBlockDataRowSparsity(ADATA,i,innz,rnnz,n);DSDPCHKERR(info);
65 for (j=0; j<i; j++){ if (rnnz[j]>0) totalnnz2++;}
66 }
67 *nnz1=totalnnz1;
68 *nnz2=totalnnz2;
69 return 0;
70}
71
72#undef __FUNCT__
73#define __FUNCT__ "CreateS1b"
74static int CreateS1b(DSDPBlockData *ADATA, int innz[], int m, int n, int tnnz[], int rnnz[], int snnz[])
75{
76 int i,j,info;
77 if (ADATA->nnzmats<=0){return 0;};
78
79 memset(rnnz,0,n*sizeof(int));
80 for (i=0;i<m;i++) innz[i]=1;
81 innz[0]=0;
82
83 /* Check matrices A for nonzero entries. */
84 for (i=0;i<n;i++){
85 memset(tnnz,0,n*sizeof(int));
86 info=DSDPBlockDataRowSparsity(ADATA,i,innz,tnnz,n);DSDPCHKERR(info);
87 for (j=0; j<=i; j++){
88 if (tnnz[j]>0){ *snnz=j; snnz++; rnnz[i]++;}
89 }
90 }
91 return 0;
92}
93
94
95#undef __FUNCT__
96#define __FUNCT__ "DSDPCreateDS"
97int DSDPCreateDS(DSDPBlockData *ADATA, DSDPVMat T, int iworkm[], int m, int n, int nnz, int rnnz[], int tnnz[], DSDPDSMat *B){
98 int i,n1,*cols,allnnz,info;
99 double *ss;
100 void *dsmat;
101 struct DSDPDSMat_Ops* dsops;
102 DSDPFunctionBegin;
103 DSDPLogInfo(0,19,"DS Matrix has %d nonzeros of %d\n",nnz,n*(n-1)/2);
104 if (nnz==0){
105 info=DSDPCreateDiagDSMatP(n,&dsops,&dsmat); DSDPCHKERR(info);
106 DSDPLogInfo(0,19,"Using Diagonal Delta S matrix\n");
107 } else if (2*nnz +n < n*n/5){
108 info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
109 cols=(int*)ss;
110 info = CreateS1b(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
111 for (allnnz=0,i=0;i<n;i++){allnnz+=rnnz[i];}
112 info = DSDPSparseMatCreatePattern2P(n,rnnz,cols,allnnz,&dsops,&dsmat); DSDPCHKERR(info);
113 info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
114 DSDPLogInfo(0,19,"Using Sparse Delta S matrix\n");
115 } else {
116 info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
117 info=DSDPCreateDSMatWithArray(n,ss,n1,&dsops,&dsmat); DSDPCHKERR(info);
118 info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
119 DSDPLogInfo(0,19,"Using Full Delta S matrix\n");
120 }
121 info=DSDPDSMatSetData(B,dsops,dsmat); DSDPCHKERR(info);
122 DSDPFunctionReturn(0);
123}
124
125
126#undef __FUNCT__
127#define __FUNCT__ "CreateS1c"
128static int CreateS1c(DSDPBlockData *ADATA, int innz[], int m, int n, int tnnz[], int rnnz[], int snnz[])
129{
130 int i,j,info;
131 memset(rnnz,0,n*sizeof(int));
132 for (i=0;i<m;i++) innz[i]=1;
133 /* Check matrix C and A for nonzero entries. */
134 for (i=0;i<n;i++){
135 memset(tnnz,0,n*sizeof(int));
136 info=DSDPBlockDataRowSparsity(ADATA,i,innz,tnnz,n);DSDPCHKERR(info);
137 for (j=i+1; j<n; j++){
138 if (tnnz[j]>0){ *snnz=j; snnz++; rnnz[i]++;}
139 }
140 }
141 return 0;
142}
143
144static int dsdpuselapack=1;
145#undef __FUNCT__
146#define __FUNCT__ "SDPConeUseLAPACKForDualMatrix"
147int SDPConeUseLAPACKForDualMatrix(SDPCone sdpcone,int flag){
148 DSDPFunctionBegin;
149 dsdpuselapack = flag;
150 DSDPFunctionReturn(0);
151}
152
153
154#undef __FUNCT__
155#define __FUNCT__ "DSDPCreateS"
156static int DSDPCreateS1(DSDPBlockData *ADATA, int trank, DSDPVec WY, DSDPVMat T, SDPConeVec W1, SDPConeVec W2, DSDPDualMat *S, DSDPDualMat *SS, DSDPDSMat *DS, void*ctx){
157
158 int nnz,n,n1,*cols,*rnnz,*tnnz,*iworkm,m,info;
159 int dsnnz,snnz,sfnnz;
160 double *pss;
161 void *smat1,*smat2;
162 struct DSDPDualMat_Ops *sops1,*sops2;
163
164 DSDPFunctionBegin;
165 info=DSDPVecGetSize(WY,&m);DSDPCHKERR(info);
166 info=DSDPVecGetArray(WY,&pss);DSDPCHKERR(info);
167 iworkm=(int*)pss;
168
169 info=SDPConeVecGetSize(W1,&n);DSDPCHKERR(info);
170 info=SDPConeVecGetArray(W1,&pss);DSDPCHKERR(info);
171 tnnz=(int*)pss;
172 info=SDPConeVecGetArray(W2,&pss);DSDPCHKERR(info);
173 rnnz=(int*)pss;
174
175 DSDPLogInfo(0,19,"Compute Sparsity\n");
176 info = CountNonzeros(ADATA, m, rnnz, iworkm, n, &dsnnz,&snnz); DSDPCHKERR(info);
177 nnz=snnz;
178 /* printf("Nonzeros: %d %d of %d\n",dsnnz,snnz,n*(n-1)/2); */
179 /* TT and DS could share memory */
180 info=DSDPCreateDS(ADATA,T,iworkm,m,n,dsnnz,rnnz,tnnz,DS);DSDPCHKERR(info);
181
182 if (nnz==0){
183 info=DSDPDiagDualMatCreateP(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
184 DSDPLogInfo(0,19,"Using Diagonal S matrix\n");
185 } else if (2*snnz+n+2<n*n/8){
186 info=DSDPVMatGetArray(T,&pss,&n1); DSDPCHKERR(info);
187 cols=(int*)pss;
188 info = CreateS1c(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
189 info=DSDPSparseDualMatCreate(n,rnnz,cols,trank,'P',&sfnnz,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
190 info=DSDPVMatRestoreArray(T,&pss,&n1); DSDPCHKERR(info);
191 /* printf("NNZ: %d %d of %d\n",2*snnz+n,sfnnz*2+n,n*n); */
192 DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Sparse S matrix\n",nnz,n*(n-1)/2);
193 DSDPLogInfo(0,19,"Total rank of block: %d, n= %d\n",trank,n);
194 } else if (n>20 && dsdpuselapack){
195 info=DSDPLAPACKSUDualMatCreate2P(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
196 DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Full Dense LAPACK S matrix\n",nnz,n*(n-1)/2);
197 } else if (dsdpuselapack){
198 info=DSDPLAPACKPUDualMatCreate2(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
199 DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Packed Dense LAPACK S matrix\n",nnz,n*(n-1)/2);
200 } else {
201 info=DSDPDenseDualMatCreate(n,'P',&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
202 DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Dense S matrix\n",nnz,n*(n-1)/2);
203 }
204 info=DSDPDualMatSetData(S,sops1,smat1); DSDPCHKERR(info);
205 info=DSDPDualMatSetData(SS,sops2,smat2); DSDPCHKERR(info);
206
207 DSDPFunctionReturn(0);
208}
209
210
211
212#undef __FUNCT__
213#define __FUNCT__ "DSDPCreateDS2"
214int DSDPCreateDS2(DSDPBlockData *ADATA, DSDPVMat T, int iworkm[], int m, int n, int nnz, int rnnz[], int tnnz[], DSDPDSMat *B){
215 int i,n1,*cols,allnnz,info;
216 double *ss;
217 void *dsmat;
218 struct DSDPDSMat_Ops* dsops;
219 DSDPFunctionBegin;
220 DSDPLogInfo(0,19,"DS Matrix has %d nonzeros of %d\n",nnz,n*(n-1)/2);
221 if (nnz==0){
222 info=DSDPCreateDiagDSMatU(n,&dsops,&dsmat); DSDPCHKERR(info);
223 DSDPLogInfo(0,19,"Using Diagonal Delta S matrix\n");
224 } else if (2*nnz +n < n*n/4){
225 info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
226 cols=(int*)ss;
227 info = CreateS1b(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
228 for (allnnz=0,i=0;i<n;i++){allnnz+=rnnz[i];}
229 info = DSDPSparseMatCreatePattern2U(n,rnnz,cols,allnnz,&dsops,&dsmat); DSDPCHKERR(info);
230 info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
231 DSDPLogInfo(0,19,"Using Sparse Delta S matrix\n");
232 } else {
233 info=DSDPVMatGetArray(T,&ss,&n1); DSDPCHKERR(info);
234 info=DSDPCreateDSMatWithArray2(n,ss,n1,&dsops,&dsmat); DSDPCHKERR(info);
235 info=DSDPVMatRestoreArray(T,&ss,&n1); DSDPCHKERR(info);
236 DSDPLogInfo(0,19,"Using Full Delta S matrix\n");
237 }
238 info=DSDPDSMatSetData(B,dsops,dsmat); DSDPCHKERR(info);
239 DSDPFunctionReturn(0);
240}
241
242#undef __FUNCT__
243#define __FUNCT__ "DSDPCreateS2"
244static int DSDPCreateS2(DSDPBlockData *ADATA, int trank, DSDPVec WY, DSDPVMat T, SDPConeVec W1, SDPConeVec W2, DSDPDualMat *S, DSDPDualMat *SS, DSDPDSMat *DS, void*ctx){
245
246 int nnz,n,n1,*cols,*rnnz,*tnnz,*iworkm,m,info;
247 int dsnnz,snnz,sfnnz;
248 double *pss;
249 void *smat1,*smat2;
250 struct DSDPDualMat_Ops *sops1,*sops2;
251 DSDPFunctionBegin;
252
253 info=DSDPVecGetSize(WY,&m);DSDPCHKERR(info);
254 info=DSDPVecGetArray(WY,&pss);DSDPCHKERR(info);
255 iworkm=(int*)pss;
256
257 info=SDPConeVecGetSize(W1,&n);DSDPCHKERR(info);
258 info=SDPConeVecGetArray(W1,&pss);DSDPCHKERR(info);
259 tnnz=(int*)pss;
260 info=SDPConeVecGetArray(W2,&pss);DSDPCHKERR(info);
261 rnnz=(int*)pss;
262
263 DSDPLogInfo(0,19,"Compute Sparsity\n");
264 info = CountNonzeros(ADATA, m, rnnz, iworkm, n, &dsnnz,&snnz); DSDPCHKERR(info);
265 nnz=snnz;
266 /* printf("Nonzeros: %d %d of %d\n",dsnnz,snnz,n*(n-1)/2); */
267 /* TT and DS could share memory */
268 info=DSDPCreateDS2(ADATA,T,iworkm,m,n,dsnnz,rnnz,tnnz,DS);DSDPCHKERR(info);
269
270 if (nnz==0){
271 info=DSDPDiagDualMatCreateU(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
272 DSDPLogInfo(0,19,"Using Diagonal S matrix\n");
273 } else if (2*snnz+n+2<n*n/10){
274 info=DSDPVMatGetArray(T,&pss,&n1); DSDPCHKERR(info);
275 cols=(int*)pss;
276 info = CreateS1c(ADATA, iworkm, m, n, tnnz, rnnz, cols); DSDPCHKERR(info);
277 info=DSDPSparseDualMatCreate(n,rnnz,cols,trank,'U',&sfnnz,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
278 info=DSDPVMatRestoreArray(T,&pss,&n1); DSDPCHKERR(info);
279 DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Sparse S matrix\n",nnz,n*(n-1)/2);
280 } else if (dsdpuselapack){
281 info=DSDPLAPACKSUDualMatCreate2(n,&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
282 DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Full Dense LAPACK S matrix\n",nnz,n*(n-1)/2);
283 } else {
284 info=DSDPDenseDualMatCreate(n,'U',&sops1,&smat1,&sops2,&smat2); DSDPCHKERR(info);
285 DSDPLogInfo(0,19,"Count %d of %d nonzeros: Using Packed Dense S matrix\n",nnz,n*(n-1)/2);
286 }
287 info=DSDPDualMatSetData(S,sops1,smat1); DSDPCHKERR(info);
288 info=DSDPDualMatSetData(SS,sops2,smat2); DSDPCHKERR(info);
289
290 DSDPFunctionReturn(0);
291}
292
293
295
296#undef __FUNCT__
297#define __FUNCT__ "DSDPCreateS"
314int DSDPCreateS(DSDPBlockData *ADATA, char UPLQ, int trank, DSDPVec WY, DSDPVMat T, SDPConeVec W1, SDPConeVec W2, DSDPDualMat *S, DSDPDualMat *SS, DSDPDSMat *DS, void*ctx){
315 int info;
316 DSDPFunctionBegin;
317 switch (UPLQ){
318 case 'P':
319 info=DSDPCreateS1(ADATA,trank,WY,T,W1,W2,S,SS,DS,ctx);DSDPCHKERR(info);
320 break;
321 case 'U':
322 info=DSDPCreateS2(ADATA,trank,WY,T,W1,W2,S,SS,DS,ctx);DSDPCHKERR(info);
323 break;
324 }
325 DSDPFunctionReturn(0);
326}
327
328
329
330#undef __FUNCT__
331#define __FUNCT__ "DSDPCreateS"
332int DSDPUseDefaultDualMatrix(SDPCone sdpcone){
333 DSDPFunctionBegin;
334 /*
335 int info;
336 info=DSDPSetDualMatrix(sdpcone,DSDPCreateS2,0); DSDPCHKERR(info);
337 */
338 DSDPFunctionReturn(0);
339}
340
341#undef __FUNCT__
342#define __FUNCT__ "DSDPMakeVMat"
351int DSDPMakeVMat(char UPLQ, int n, DSDPVMat *X){
352 int info;
353 void *xmat;
354 struct DSDPVMat_Ops* xops;
355 DSDPFunctionBegin;
356 switch (UPLQ){
357 case 'P':
358 info=DSDPXMatPCreate(n,&xops,&xmat);DSDPCHKERR(info);
359 break;
360 case 'U':
361 info=DSDPXMatUCreate(n,&xops,&xmat);DSDPCHKERR(info);
362 break;
363 }
364 info=DSDPVMatSetData(X,xops,xmat); DSDPCHKERR(info);
365 DSDPFunctionReturn(0);
366}
367
368#undef __FUNCT__
369#define __FUNCT__ "DSDPMakeVMatWithArray"
381int DSDPMakeVMatWithArray(char UPLQ, double xx[], int nnz, int n, DSDPVMat *X){
382 int info;
383 void *xmat;
384 struct DSDPVMat_Ops* xops;
385 DSDPFunctionBegin;
386 switch (UPLQ){
387 case 'P':
388 info=DSDPXMatPCreateWithData(n,xx,nnz,&xops,&xmat);DSDPCHKERR(info);
389 break;
390 case 'U':
391 info=DSDPXMatUCreateWithData(n,xx,nnz,&xops,&xmat);DSDPCHKERR(info);
392 break;
393 }
394 info=DSDPVMatSetData(X,xops,xmat); DSDPCHKERR(info);
395 DSDPFunctionReturn(0);
396}
int DSDPBlockDataRowSparsity(DSDPBlockData *ADATA, int row, int ai[], int rnnz[], int n)
Determine sparsity pattern of data.
Definition dsdpblock.c:330
int DSDPDSMatSetData(DSDPDSMat *M, struct DSDPDSMat_Ops *ops, void *data)
Set the opaque pointer and function pointers to the matrix.
Definition dsdpdsmat.c:31
int DSDPDualMatSetData(DSDPDualMat *S, struct DSDPDualMat_Ops *ops, void *data)
Set the opaque pointer and function pointers to the matrix.
Definition dsdpdualmat.c:49
Internal SDPCone data structures and routines.
Error handling, printing, and profiling.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition dsdpvec.h:25
int DSDPVMatSetData(DSDPVMat *X, struct DSDPVMat_Ops *ops, void *data)
Set opaque pointer an function pointers.
Definition dsdpxmat.c:39
int DSDPVMatGetArray(DSDPVMat X, double **v, int *nn)
Get the array that stores the matrix.
Definition dsdpxmat.c:211
int DSDPVMatRestoreArray(DSDPVMat X, double **v, int *nn)
Restore the array that stores the matrix.
Definition dsdpxmat.c:233
int DSDPCreateS(DSDPBlockData *, char, int, DSDPVec, DSDPVMat, SDPConeVec, SDPConeVec, DSDPDualMat *, DSDPDualMat *, DSDPDSMat *, void *)
Create S1, S2, and DS.
Definition sdpsss.c:314
int DSDPMakeVMat(char UPLQ, int n, DSDPVMat *X)
Allocate V matrix.
Definition sdpsss.c:351
int DSDPMakeVMatWithArray(char UPLQ, double xx[], int nnz, int n, DSDPVMat *X)
Allocate V matrix using the given array.
Definition sdpsss.c:381
Internal structure for data in one block of semidefintie.
Definition dsdpsdp.h:39
Symmetric Delta S matrix for one block in the semidefinite cone.
Definition dsdpdsmat.h:23
Symmetric Delta S matrix for one block in the semidefinite cone.
Represents an S matrix for one block in the semidefinite cone.
Definition dsdpdualmat.h:18
Table of function pointers that operate on the S matrix.
Dense symmetric matrix for one block in the semidefinite cone.
Definition dsdpxmat.h:17
Table of function pointers that operate on the dense matrix.
Vector whose length corresponds to dimension of a block in a cone.
Definition sdpconevec.h:13
Internal structure for semidefinite cone.
Definition dsdpsdp.h:80