DSDP
sdpconesetup.c
Go to the documentation of this file.
1 #include "dsdpsdp.h"
2 #include "dsdpsys.h"
8 #undef __FUNCT__
9 #define __FUNCT__ "DSDPDataTransposeInitialize"
10 
16  DSDPFunctionBegin;
17  ATranspose->nnzblocks=0;
18  ATranspose->nzblocks=0;
19  ATranspose->idA=0;
20  ATranspose->idAP=0;
21  ATranspose->ttnzmat=0;
22  ATranspose->nnzblocks=0;
23  DSDPFunctionReturn(0);
24 }
25 
26 #undef __FUNCT__
27 #define __FUNCT__ "DSDPDataTransposeSetup"
28 
36 int DSDPDataTransposeSetup(DSDPDataTranspose *ATranspose, SDPblk *blk, int nblocks, int m){
37 
38  int i,ii,kk,vvar,info;
39  int nnzmats,tnzmats=0;
40  DSDPFunctionBegin;
41 
42  info=DSDPDataTransposeTakeDown(ATranspose);DSDPCHKERR(info);
43  /* Determine sparsity pattern of SDP Data Matrices */
44 
45  DSDPCALLOC2(&ATranspose->nnzblocks,int,(m),&info);DSDPCHKERR(info);
46  DSDPCALLOC2(&ATranspose->nzblocks,int*,(m),&info);DSDPCHKERR(info);
47  DSDPCALLOC2(&ATranspose->idA,int*,(m),&info);DSDPCHKERR(info);
48  ATranspose->m=m;
49  for (i=0;i<m;i++){ ATranspose->nnzblocks[i]=0; }
50  for (kk=0; kk<nblocks; kk++){
51  info=DSDPBlockDataMarkNonzeroMatrices(&blk[kk].ADATA,ATranspose->nnzblocks);DSDPCHKERR(info);
52  }
53  for (tnzmats=0,i=0;i<m;i++){ tnzmats += ATranspose->nnzblocks[i];}
54 
55  DSDPCALLOC2(&ATranspose->ttnzmat,int,tnzmats,&info);DSDPCHKERR(info);
56  ATranspose->nzblocks[0]=ATranspose->ttnzmat;
57  for (i=1;i<m;i++){
58  ATranspose->nzblocks[i]=ATranspose->nzblocks[i-1]+ATranspose->nnzblocks[i-1];
59  }
60 
61  DSDPCALLOC2(&ATranspose->idAP,int,tnzmats,&info);DSDPCHKERR(info);
62  ATranspose->idA[0]=ATranspose->idAP;
63  for (i=1;i<m;i++){
64  ATranspose->idA[i]=ATranspose->idA[i-1]+ATranspose->nnzblocks[i-1];
65  }
66 
67  for (i=0;i<m;i++){ATranspose->nnzblocks[i]=0;}
68  for (kk=0; kk<nblocks; kk++){
69  info=DSDPBlockCountNonzeroMatrices(&blk[kk].ADATA,&nnzmats);DSDPCHKERR(info);
70  for (i=0;i<nnzmats;i++){
71  info=DSDPBlockGetMatrix(&blk[kk].ADATA,i,&ii,0,0);DSDPCHKERR(info);
72  vvar=ATranspose->nnzblocks[ii];
73  ATranspose->nzblocks[ii][vvar]=kk;
74  ATranspose->idA[ii][vvar]=i;
75  ATranspose->nnzblocks[ii]++;
76  }
77  }
78 
79  DSDPFunctionReturn(0);
80 }
81 
82 #undef __FUNCT__
83 #define __FUNCT__ "DSDPDataTransposeTakeDown"
84 
90  int info;
91  DSDPFunctionBegin;
92  DSDPFREE(&ATranspose->ttnzmat,&info);DSDPCHKERR(info);
93  DSDPFREE(&ATranspose->idAP,&info);DSDPCHKERR(info);
94  DSDPFREE(&ATranspose->nzblocks,&info);DSDPCHKERR(info);
95  DSDPFREE(&ATranspose->nnzblocks,&info);DSDPCHKERR(info);
96  DSDPFREE(&ATranspose->idA,&info);DSDPCHKERR(info);
97  info=DSDPDataTransposeInitialize(ATranspose);DSDPCHKERR(info);
98  DSDPFunctionReturn(0);
99 }
100 
101 #undef __FUNCT__
102 #define __FUNCT__ "DSDPCreateSDPCone"
103 
113 int DSDPCreateSDPCone(DSDP dsdp, int blocks, SDPCone* dspcone){
114  int i,info;
115  SDPCone sdpcone;
116 
117  DSDPFunctionBegin;
118  DSDPCALLOC1(&sdpcone,struct SDPCone_C,&info);DSDPCHKERR(info);
119  *dspcone=sdpcone;
120  sdpcone->keyid=SDPCONEKEY;
121  info=DSDPAddSDP(dsdp,sdpcone);DSDPCHKERR(info);
122 
123  info=DSDPGetNumberOfVariables(dsdp,&sdpcone->m);DSDPCHKERR(info);
124  DSDPCALLOC2(&sdpcone->blk,SDPblk,blocks,&info); DSDPCHKERR(info);
125  for (i=0;i<blocks; i++){
126  info=DSDPBlockInitialize(&sdpcone->blk[i]); DSDPCHKBLOCKERR(i,info);
127  }
128 
129  sdpcone->nblocks=blocks;
130  sdpcone->optype=3;
131  info=DSDPUseDefaultDualMatrix(sdpcone); DSDPCHKERR(info);
132 
133  sdpcone->nn=0;
134  sdpcone->dsdp=dsdp;
135  info=DSDPDataTransposeInitialize(&sdpcone->ATR); DSDPCHKERR(info);
136  info=DSDPBlockEventZero();DSDPCHKERR(info);
137  info=DSDPDualMatEventZero();DSDPCHKERR(info);
138  info=DSDPVMatEventZero();DSDPCHKERR(info);
139  DSDPFunctionReturn(0);
140 }
141 
142 
144 
145 #undef __FUNCT__
146 #define __FUNCT__ "DSDPBlockSetup"
147 
154 int DSDPBlockSetup(SDPblk *blk, int blockj, DSDPVec WY){
155  int n,info,trank,flag;
156  DSDPFunctionBegin;
157  /*
158  info=DSDPBlockTakeDown(blk);DSDPCHKERR(info);
159  */
160  n=blk->n;
161  info=DSDPVMatExist(blk->T,&flag);DSDPCHKERR(info);
162  if (flag==0){
163  info=DSDPMakeVMat(blk->format,n,&blk->T);DSDPCHKERR(info);
164  }
165 
166  info = DSDPIndexCreate(blk->n,&blk->IS);DSDPCHKERR(info);
167  info = SDPConeVecCreate(blk->n,&blk->W);DSDPCHKERR(info);
168  info = SDPConeVecDuplicate(blk->W,&blk->W2);DSDPCHKERR(info);
169 
170  /* Build Lanczos structures */
171  info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,20); DSDPCHKERR(info);
172  if (n>30){info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,20); DSDPCHKERR(info);}
173  if (n>300){info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,40); DSDPCHKERR(info);}
174  if (n>1000){info=DSDPSetMaximumLanczosIterations(&blk->Lanczos,50); DSDPCHKERR(info);}
175  if (1){
176  info=DSDPFastLanczosSetup(&blk->Lanczos,blk->W);DSDPCHKERR(info);
177  DSDPLogInfo(0,19,"SDP Block %d using Fast Lanczos\n",blockj);
178  } else {
179  info=DSDPRobustLanczosSetup(&blk->Lanczos,blk->W);DSDPCHKERR(info);
180  DSDPLogInfo(0,19,"SDP Block %d using Full Lanczos\n",blockj);
181  }
182 
183  /* Eigenvalues and Eigenvectors */
184  info=DSDPBlockFactorData(&blk->ADATA,blk->T,blk->W);DSDPCHKERR(info);
185  info=DSDPBlockDataRank(&blk->ADATA,&trank,n);DSDPCHKERR(info);
186 
187  info=DSDPCreateS(&blk->ADATA,blk->format,trank,WY,blk->T,blk->W,blk->W2,&blk->S,&blk->SS,&blk->DS,0);DSDPCHKERR(info);
188 
189  DSDPFunctionReturn(0);
190 }
191 
192 #undef __FUNCT__
193 #define __FUNCT__ "SDPConeBlockNNZ"
194 int SDPConeBlockNNZ(SDPblk *blk,int m){
195  int i,ii,n,info,nnz,nnzmats,tnnzmats,tnnz=0;
196  double scl;
197  DSDPDataMat AA;
198  DSDPFunctionBegin;
199  nnzmats=blk->ADATA.nnzmats;tnnzmats=nnzmats;
200  n=blk->n;
201 
202  for (i=0;i<nnzmats;i++){
203  info=DSDPBlockGetMatrix(&blk->ADATA,i,&ii,&scl,&AA);DSDPCHKERR(info);
204  if (ii==0){tnnzmats--; continue;}
205  if (ii==m-1){continue;}
206  info = DSDPDataMatCountNonzeros(AA,&nnz,n); DSDPCHKERR(info);
207  tnnz+= (nnz*(tnnzmats-i));
208  }
209  if (tnnzmats>1){ tnnz=tnnz/((tnnzmats)*(tnnzmats+1)/2); }
210  if (tnnz<1) tnnz = 1;
211  blk->nnz=tnnz;
212  DSDPFunctionReturn(0);
213 }
214 
215 #undef __FUNCT__
216 #define __FUNCT__ "SDPConeSetup2"
217 
225  int kk,n,m,info;
226  double nn=0;
227  SDPblk *blk;
228  DSDPFunctionBegin;
229  info=DSDPVecGetSize(yy0,&m);DSDPCHKERR(info);
230  for (kk=0; kk<sdpcone->nblocks; kk++){
231  blk=&sdpcone->blk[kk];
232  n=blk->n;
233  info=SDPConeBlockNNZ(blk,m);DSDPCHKERR(info);
234  info=DSDPBlockSetup(blk,kk,sdpcone->Work);DSDPCHKERR(info);
235  nn+=n*blk->gammamu;
236  }
237  sdpcone->nn=(int)nn;
238  DSDPFunctionReturn(0);
239 }
240 
241 #undef __FUNCT__
242 #define __FUNCT__ "SDPConeSetup"
243 
249 int SDPConeSetup(SDPCone sdpcone, DSDPVec yy0){
250  int kk,n,m,info;
251  DSDPFunctionBegin;
252 
253  info = DSDPVecGetSize(yy0,&m);DSDPCHKERR(info);
254  if (m!=sdpcone->m+2){DSDPSETERR(8,"CHECK DIMENSION\n");}
255  info = DSDPVecDuplicate(yy0,&sdpcone->Work);DSDPCHKERR(info);
256  info = DSDPVecDuplicate(yy0,&sdpcone->Work2);DSDPCHKERR(info);
257  info = DSDPVecDuplicate(yy0,&sdpcone->YY);DSDPCHKERR(info);
258  info = DSDPVecDuplicate(yy0,&sdpcone->YX);DSDPCHKERR(info);
259  info = DSDPVecDuplicate(yy0,&sdpcone->DYX);DSDPCHKERR(info);
260  for (kk=0; kk<sdpcone->nblocks; kk++){
261  n=sdpcone->blk[kk].n;
262  info=SDPConeSetRIdentity(sdpcone,kk,n,1.0);DSDPCHKERR(info);
263  }
264 
265  info=DSDPDataTransposeSetup(&sdpcone->ATR,sdpcone->blk,sdpcone->nblocks,m); DSDPCHKERR(info);
266  info=DSDPBlockEventInitialize();DSDPCHKERR(info);
267  info=DSDPDualMatEventInitialize();DSDPCHKERR(info);
268  info=DSDPVMatEventInitialize();DSDPCHKERR(info);
269  DSDPFunctionReturn(0);
270 }
271 
272 #undef __FUNCT__
273 #define __FUNCT__ "DSDPBlockInitialize"
274 
280  int info;
281  DSDPFunctionBegin;
282  blk->n=0;
283  blk->format='N';
284  blk->gammamu=1.0;
285  blk->bmu=0.0;
286  blk->nnz=10000000;
287 
288  info = DSDPDualMatInitialize(&blk->S); DSDPCHKERR(info);
289  info = DSDPDualMatInitialize(&blk->SS); DSDPCHKERR(info);
290  info = DSDPDSMatInitialize(&blk->DS); DSDPCHKERR(info);
291  info = DSDPVMatInitialize(&blk->T); DSDPCHKERR(info);
292  info = DSDPLanczosInitialize(&blk->Lanczos); DSDPCHKERR(info);
293  info = DSDPBlockDataInitialize(&blk->ADATA); DSDPCHKERR(info);
294  info = DSDPIndexInitialize(&blk->IS); DSDPCHKERR(info);
295  DSDPFunctionReturn(0);
296 }
297 
298 #undef __FUNCT__
299 #define __FUNCT__ "DSDPBlockTakeDown"
300 
306  int info;
307  DSDPFunctionBegin;
308  if (!blk){DSDPFunctionReturn(0);}
309  info=DSDPBlockTakeDownData(&blk->ADATA);DSDPCHKERR(info);
310  info=SDPConeVecDestroy(&blk->W);DSDPCHKERR(info);
311  info=SDPConeVecDestroy(&blk->W2);DSDPCHKERR(info);
312  info=DSDPIndexDestroy(&blk->IS);DSDPCHKERR(info);
313  info=DSDPLanczosDestroy(&blk->Lanczos);DSDPCHKERR(info);
314  info=DSDPDualMatDestroy(&blk->SS);DSDPCHKERR(info);
315  info=DSDPDualMatDestroy(&blk->S);DSDPCHKERR(info);
316  info=DSDPDSMatDestroy(&blk->DS);DSDPCHKERR(info);
317  info=DSDPVMatDestroy(&blk->T);DSDPCHKERR(info);
318  DSDPFunctionReturn(0);
319 }
320 
321 #undef __FUNCT__
322 #define __FUNCT__ "DSDPConeTakeDown"
323 
329  int blockj,info;
330  DSDPFunctionBegin;
331  for (blockj=0; blockj<sdpcone->nblocks; blockj++){
332  info=DSDPBlockTakeDown(&sdpcone->blk[blockj]);DSDPCHKERR(info);
333  }
334  info=DSDPVecDestroy(&sdpcone->Work);DSDPCHKERR(info);
335  info=DSDPVecDestroy(&sdpcone->Work2);DSDPCHKERR(info);
336  info=DSDPVecDestroy(&sdpcone->YY);DSDPCHKERR(info);
337  info=DSDPVecDestroy(&sdpcone->YX);DSDPCHKERR(info);
338  info=DSDPVecDestroy(&sdpcone->DYX);DSDPCHKERR(info);
339  info=DSDPDataTransposeTakeDown(&sdpcone->ATR);DSDPCHKERR(info);
340  DSDPFunctionReturn(0);
341 }
342 
343 #undef __FUNCT__
344 #define __FUNCT__ "SDPConeDestroy"
345 
350 int SDPConeDestroy(SDPCone sdpcone){
351  int blockj,info;
352  DSDPFunctionBegin;
353  info=DSDPConeTakeDown(sdpcone);DSDPCHKERR(info);
354  for (blockj=0; blockj<sdpcone->nblocks; blockj++){
355  info=DSDPBlockDataDestroy(&sdpcone->blk[blockj].ADATA);DSDPCHKERR(info);
356  }
357  DSDPFREE(&sdpcone->blk,&info);DSDPCHKERR(info);
358  DSDPFREE(&sdpcone,&info);DSDPCHKERR(info);
359  info=DSDPBlockEventZero();DSDPCHKERR(info);
360  info=DSDPDualMatEventZero();DSDPCHKERR(info);
361  info=DSDPVMatEventZero();DSDPCHKERR(info);
362  DSDPFunctionReturn(0);
363 }
364