DSDP
sdpsss.c
Go to the documentation of this file.
1 #include "dsdpsdp.h"
2 #include "dsdpsys.h"
9 
10 static int DSDPCreateDS(DSDPBlockData*, DSDPVMat,int*,int,int,int,int*,int*,DSDPDSMat*);
11 static int DSDPCreateDS2(DSDPBlockData*, DSDPVMat,int*,int,int,int,int*,int*,DSDPDSMat*);
12 
13 static int DSDPCreateS1(DSDPBlockData*,int,DSDPVec,DSDPVMat,SDPConeVec, SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
14 static int DSDPCreateS2(DSDPBlockData*,int,DSDPVec,DSDPVMat,SDPConeVec, SDPConeVec,DSDPDualMat*, DSDPDualMat*, DSDPDSMat*, void*);
15 
16 extern int DSDPXMatPCreate(int,struct DSDPVMat_Ops**,void**);
17 extern int DSDPXMatPCreateWithData(int,double[],int,struct DSDPVMat_Ops**,void**);
18 extern int DSDPXMatUCreate(int,struct DSDPVMat_Ops**,void**);
19 extern int DSDPXMatUCreateWithData(int,double[],int,struct DSDPVMat_Ops**,void**);
20 
21 extern int DSDPLAPACKSUDualMatCreate2(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
22 extern int DSDPLAPACKSUDualMatCreate2P(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
23 extern int DSDPLAPACKPUDualMatCreate2(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
24 extern int DSDPDiagDualMatCreateP(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
25 extern int DSDPDiagDualMatCreateU(int,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
26 extern int DSDPDenseDualMatCreate(int, char,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
27 extern int DSDPSparseDualMatCreate(int, int*, int*, int,char,int*,struct DSDPDualMat_Ops**,void**, struct DSDPDualMat_Ops**,void**);
28 
29 extern int DSDPSparseMatCreatePattern2P(int,int[],int[],int,struct DSDPDSMat_Ops**,void**);
30 extern int DSDPSparseMatCreatePattern2U(int,int[],int[],int,struct DSDPDSMat_Ops**,void**);
31 
32 extern int DSDPCreateDiagDSMatP(int,struct DSDPDSMat_Ops**,void**);
33 extern int DSDPCreateDiagDSMatU(int,struct DSDPDSMat_Ops**,void**);
34 
35 extern int DSDPCreateDSMatWithArray(int,double[],int,struct DSDPDSMat_Ops**,void**);
36 extern int DSDPCreateDSMatWithArray2(int,double[],int, struct DSDPDSMat_Ops**,void**);
37 
38 extern int DSDPCreateURDSMat(int,struct DSDPDSMat_Ops**,void**);
39 
40 /*
41 #undef __FUNCT__
42 #define __FUNCT__ "DSDPSetDualMatrix"
43 int 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"
52 static 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"
74 static 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"
97 int 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"
128 static 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 
144 static int dsdpuselapack=1;
145 #undef __FUNCT__
146 #define __FUNCT__ "SDPConeUseLAPACKForDualMatrix"
147 int SDPConeUseLAPACKForDualMatrix(SDPCone sdpcone,int flag){
148  DSDPFunctionBegin;
149  dsdpuselapack = flag;
150  DSDPFunctionReturn(0);
151 }
152 
153 
154 #undef __FUNCT__
155 #define __FUNCT__ "DSDPCreateS"
156 static 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"
214 int 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"
244 static 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"
298 
314 int 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"
332 int 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"
343 
351 int 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"
370 
381 int 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 }
Table of function pointers that operate on the dense matrix.
Definition: dsdpxmat_impl.h:13
int DSDPMakeVMat(char UPLQ, int n, DSDPVMat *X)
Allocate V matrix.
Definition: sdpsss.c:351
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
int DSDPVMatGetArray(DSDPVMat X, double **v, int *nn)
Get the array that stores the matrix.
Definition: dsdpxmat.c:211
Error handling, printing, and profiling.
int DSDPVMatRestoreArray(DSDPVMat X, double **v, int *nn)
Restore the array that stores the matrix.
Definition: dsdpxmat.c:233
Internal structure for data in one block of semidefintie.
Definition: dsdpsdp.h:39
Vector whose length corresponds to dimension of a block in a cone.
Definition: sdpconevec.h:13
Table of function pointers that operate on the S matrix.
int DSDPCreateS(DSDPBlockData *, char, int, DSDPVec, DSDPVMat, SDPConeVec, SDPConeVec, DSDPDualMat *, DSDPDualMat *, DSDPDSMat *, void *)
Create S1, S2, and DS.
Definition: sdpsss.c:314
int DSDPMakeVMatWithArray(char UPLQ, double xx[], int nnz, int n, DSDPVMat *X)
Allocate V matrix using the given array.
Definition: sdpsss.c:381
Represents an S matrix for one block in the semidefinite cone.
Definition: dsdpdualmat.h:18
int DSDPBlockDataRowSparsity(DSDPBlockData *ADATA, int row, int ai[], int rnnz[], int n)
Determine sparsity pattern of data.
Definition: dsdpblock.c:330
Internal structure for semidefinite cone.
Definition: dsdpsdp.h:80
Symmetric Delta S matrix for one block in the semidefinite cone.
Symmetric Delta S matrix for one block in the semidefinite cone.
Definition: dsdpdsmat.h:23
Dense symmetric matrix for one block in the semidefinite cone.
Definition: dsdpxmat.h:17
Internal SDPCone data structures and routines.