DSDP
dsdplp.c
Go to the documentation of this file.
1
2#include "dsdpcone_impl.h"
3#include "dsdpsys.h"
4#include "dsdp5.h"
5
6#include <stdio.h>
7#include <stdlib.h>
8#include <math.h>
9#include <string.h>
14typedef struct {
15 int nrow;
16 int ncol;
17 int owndata;
18
19 const double *an;
20 const int *col;
21 const int *nnz;
22 int *nzrows;
23 int nnzrows;
24
25} smatx;
26
31struct LPCone_C{
32 smatx *A,*AT;
33 DSDPVec C;
34 DSDPVec PS,DS,X;
35 double sscale;
36 double r;
37 double muscale;
38 DSDPVec Y,WY,WY2,WX,WX2;
39 double *xout;
40 int n,m;
41};
42
43
44static int CreateSpRowMatWdata(int,int,const double[], const int[], const int[],smatx **);
45static int SpRowMatMult(smatx*,const double[], int , double[], int);
46static int SpRowMatMultTrans(smatx *, const double[],int, double[],int);
47static int SpRowMatGetRowVector(smatx*, int, double*,int);
48static int SpRowMatGetScaledRowVector(smatx*, int, const double[], double*, int);
49static int SpRowMatDestroy(smatx*);
50static int SpRowMatView(smatx*);
51/*
52static int SpRowMatGetSize(smatx *, int *, int *);
53static int SpRowMatZero(smatx*);
54static int SpRowMatAddRowMultiple(smatx*, int, double, double[], int);
55*/
56static int SpRowIMultAdd(smatx*,int*,int,int *,int);
57static int SpRowMatRowNnz(smatx*, int, int*,int);
58static int SpRowMatNorm2(smatx*, int, double*);
59
60
61#undef __FUNCT__
62#define __FUNCT__ "LPConeSetUp"
63static int LPConeSetup(void *dcone,DSDPVec y){
64 int m,info;
65 LPCone lpcone=(LPCone)dcone;
66 DSDPFunctionBegin;
67 if (lpcone->n<1) return 0;
68 m=lpcone->m;
69 info=DSDPVecCreateSeq(m+2,&lpcone->WY);DSDPCHKERR(info);
70 info=DSDPVecDuplicate(lpcone->WY,&lpcone->WY2);DSDPCHKERR(info);
71 info=DSDPVecDuplicate(lpcone->WY,&lpcone->Y);DSDPCHKERR(info);
72 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX);DSDPCHKERR(info);
73 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX2);DSDPCHKERR(info);
74 info=DSDPVecDuplicate(lpcone->C,&lpcone->PS);DSDPCHKERR(info);
75 info=DSDPVecDuplicate(lpcone->C,&lpcone->DS);DSDPCHKERR(info);
76 info=DSDPVecDuplicate(lpcone->C,&lpcone->X);DSDPCHKERR(info);
77 DSDPFunctionReturn(0);
78}
79
80#undef __FUNCT__
81#define __FUNCT__ "LPConeSetUp2"
82static int LPConeSetup2(void *dcone, DSDPVec Y, DSDPSchurMat M){
83 LPCone lpcone=(LPCone)dcone;
84 DSDPFunctionBegin;
85 DSDPLogInfo(0,19,"Setup LP Cone of dimension: %d\n",lpcone->n);
86 DSDPFunctionReturn(0);
87}
88
89
90#undef __FUNCT__
91#define __FUNCT__ "LPConeDestroy"
92static int LPConeDestroy(void *dcone){
93 int info;
94 LPCone lpcone=(LPCone)dcone;
95 DSDPFunctionBegin;
96 if (lpcone->n<1) return 0;
97 info=DSDPVecDestroy(&lpcone->DS);DSDPCHKERR(info);
98 info=DSDPVecDestroy(&lpcone->PS);DSDPCHKERR(info);
99 info=DSDPVecDestroy(&lpcone->C);DSDPCHKERR(info);
100 info=DSDPVecDestroy(&lpcone->X);DSDPCHKERR(info);
101 info=SpRowMatDestroy(lpcone->A);DSDPCHKERR(info);
102 info=DSDPVecDestroy(&lpcone->WX2);DSDPCHKERR(info);
103 info=DSDPVecDestroy(&lpcone->WY2);DSDPCHKERR(info);
104 info=DSDPVecDestroy(&lpcone->WY);DSDPCHKERR(info);
105 info=DSDPVecDestroy(&lpcone->Y);DSDPCHKERR(info);
106 info=DSDPVecDestroy(&lpcone->WX);DSDPCHKERR(info);
107 DSDPFREE(&lpcone,&info);DSDPCHKERR(info);
108 DSDPFunctionReturn(0);
109}
110
111#undef __FUNCT__
112#define __FUNCT__ "LPConeSize"
113static int LPConeSize(void *dcone, double *n){
114 LPCone lpcone=(LPCone)dcone;
115 DSDPFunctionBegin;
116 *n=lpcone->muscale*lpcone->n;
117 DSDPFunctionReturn(0);
118}
119
120
121
122#undef __FUNCT__
123#define __FUNCT__ "LPComputeAX"
124static int LPComputeAX( LPCone lpcone,DSDPVec X, DSDPVec Y){
125 int info,m=lpcone->m,n=lpcone->n;
126 double *x,*y,ppobj;
127 smatx *A=lpcone->A;
128 DSDPFunctionBegin;
129 if (lpcone->n<1) return 0;
130 info=DSDPVecGetSize(X,&n);DSDPCHKERR(info);
131 info=DSDPVecDot(lpcone->C,X,&ppobj);DSDPCHKERR(info);
132 info=DSDPVecSetC(Y,ppobj);
133 info=DSDPVecSum(X,&ppobj);DSDPCHKERR(info);
134 info=DSDPVecSetR(Y,ppobj*lpcone->r);DSDPCHKERR(info);
135 info=DSDPVecGetArray(Y,&y);DSDPCHKERR(info);
136 info=DSDPVecGetArray(X,&x);DSDPCHKERR(info);
137 info=SpRowMatMult(A,x,n,y+1,m);
138 info=DSDPVecRestoreArray(X,&x);DSDPCHKERR(info);
139 info=DSDPVecRestoreArray(Y,&y);DSDPCHKERR(info);
140 DSDPFunctionReturn(0);
141}
142
143#undef __FUNCT__
144#define __FUNCT__ "LPComputeATY"
145static int LPComputeATY(LPCone lpcone,DSDPVec Y, DSDPVec S){
146 int info,m=lpcone->m,n=lpcone->n;
147 double cc,r,*s,*y;
148 DSDPVec C=lpcone->C;
149 smatx *A=lpcone->A;
150 DSDPFunctionBegin;
151 if (lpcone->n<1) return 0;
152 info=DSDPVecGetC(Y,&cc);DSDPCHKERR(info);
153 info=DSDPVecGetR(Y,&r);DSDPCHKERR(info);
154 info=DSDPVecGetSize(S,&n);DSDPCHKERR(info);
155 info=DSDPVecGetArray(S,&s);DSDPCHKERR(info);
156 info=DSDPVecGetArray(Y,&y);DSDPCHKERR(info);
157 info=SpRowMatMultTrans(A,y+1,m,s,n); DSDPCHKERR(info);
158 info=DSDPVecRestoreArray(S,&s);DSDPCHKERR(info);
159 info=DSDPVecRestoreArray(Y,&y);DSDPCHKERR(info);
160 info=DSDPVecAXPY(cc,C,S);DSDPCHKERR(info);
161 info=DSDPVecShift(r*lpcone->r,S);DSDPCHKERR(info);
162 info=DSDPVecScale(-1.0,S);DSDPCHKERR(info);
163 DSDPFunctionReturn(0);
164}
165#undef __FUNCT__
166#define __FUNCT__ "LPConeHessian"
167static int LPConeHessian(void* dcone, double mu, DSDPSchurMat M,
168 DSDPVec vrhs1, DSDPVec vrhs2){
169 int info,i,m,n,ncols;
170 double r=1.0,*wx,*wx2;
171 LPCone lpcone=(LPCone)dcone;
172 DSDPVec WX=lpcone->WX,WX2=lpcone->WX2,WY=lpcone->WY,WY2=lpcone->WY2,S=lpcone->DS;
173 smatx *A=lpcone->A;
174
175 DSDPFunctionBegin;
176 if (lpcone->n<1) return 0;
177 mu*=lpcone->muscale;
178 info=DSDPVecGetSize(vrhs1,&m);DSDPCHKERR(info);
179 info=DSDPVecGetSize(WX,&n);DSDPCHKERR(info);
180 info=DSDPVecSet(mu,WX2);DSDPCHKERR(info);
181 info=DSDPVecPointwiseDivide(WX2,S,WX2);DSDPCHKERR(info);
182 info=DSDPVecPointwiseDivide(WX2,S,WX2);DSDPCHKERR(info);
183 for (i=0;i<m;i++){
184
185 info=DSDPSchurMatRowColumnScaling(M,i,WY2,&ncols); DSDPCHKERR(info);
186 if (ncols==0) continue;
187
188 if (i==0){
189 info=DSDPVecPointwiseMult(lpcone->C,WX2,WX); DSDPCHKERR(info);
190 } else if (i==m-1){
191 info=DSDPVecScaleCopy(WX2,r,WX); DSDPCHKERR(info);
192 } else {
193 info=DSDPVecGetArray(WX,&wx);
194 info=DSDPVecGetArray(WX2,&wx2);DSDPCHKERR(info);
195 info=SpRowMatGetScaledRowVector(A,i-1,wx2,wx,n);
196 info=DSDPVecRestoreArray(WX,&wx);
197 info=DSDPVecRestoreArray(WX2,&wx2);
198 }
199
200 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info);
201
202 info=DSDPVecPointwiseMult(WY2,WY,WY);DSDPCHKERR(info);
203
204 info=DSDPSchurMatAddRow(M,i,1.0,WY);DSDPCHKERR(info);
205 }
206
207 /* Compute AS^{-1} */
208 info=DSDPVecSet(mu,WX);DSDPCHKERR(info);
209 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info);
210 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info);
211
212 info=DSDPSchurMatDiagonalScaling(M, WY2);DSDPCHKERR(info);
213 info=DSDPVecPointwiseMult(WY2,WY,WY);DSDPCHKERR(info);
214 info=DSDPVecAXPY(1.0,WY,vrhs2);DSDPCHKERR(info);
215
216 DSDPFunctionReturn(0);
217}
218
219#undef __FUNCT__
220#define __FUNCT__ "LPConeHessian"
221static int LPConeRHS(void* dcone, double mu, DSDPVec vrow,
222 DSDPVec vrhs1, DSDPVec vrhs2){
223 int info;
224 LPCone lpcone=(LPCone)dcone;
225 DSDPVec WX=lpcone->WX,WY=lpcone->WY,S=lpcone->DS;
226
227 DSDPFunctionBegin;
228 if (lpcone->n<1) return 0;
229 mu*=lpcone->muscale;
230
231 /* Compute AS^{-1} */
232 info=DSDPVecSet(mu,WX);DSDPCHKERR(info);
233 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info);
234 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info);
235
236 info=DSDPVecPointwiseMult(vrow,WY,WY);DSDPCHKERR(info);
237 info=DSDPVecAXPY(1.0,WY,vrhs2);DSDPCHKERR(info);
238
239 DSDPFunctionReturn(0);
240}
241
242#undef __FUNCT__
243#define __FUNCT__ "LPConeMultiply"
244static int LPConeMultiply(void* dcone, double mu, DSDPVec vrow, DSDPVec vin, DSDPVec vout){
245 int info;
246 LPCone lpcone=(LPCone)dcone;
247 DSDPVec WX=lpcone->WX,S=lpcone->DS,WY=lpcone->WY;
248 DSDPFunctionBegin;
249 if (lpcone->n<1) return 0;
250 mu*=lpcone->muscale;
251 info=LPComputeATY(lpcone,vin,WX);DSDPCHKERR(info);
252 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info);
253 info=DSDPVecScale(mu,WX);DSDPCHKERR(info);
254 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info);
255 info=LPComputeAX(lpcone,WX,WY);DSDPCHKERR(info);
256 info=DSDPVecPointwiseMult(WY,vrow,WY);DSDPCHKERR(info);
257 info=DSDPVecAXPY(1.0,WY,vout);DSDPCHKERR(info);
258 DSDPFunctionReturn(0);
259}
260
261#undef __FUNCT__
262#define __FUNCT__ "LPConeSetX"
263static int LPConeSetX(void* dcone,double mu, DSDPVec Y,DSDPVec DY){
264 DSDPFunctionBegin;
265 DSDPFunctionReturn(0);
266}
267
268#undef __FUNCT__
269#define __FUNCT__ "LPConeX"
270static int LPConeX(void* dcone,double mu, DSDPVec Y,DSDPVec DY,
271 DSDPVec AX , double *tracexs){
272 int info;
273 double dtracexs;
274 LPCone lpcone=(LPCone)dcone;
275 DSDPVec S=lpcone->PS,WX=lpcone->WX,X=lpcone->X,DS=lpcone->DS,WY=lpcone->WY;
276 double *xx,*xout=lpcone->xout;
277 DSDPFunctionBegin;
278
279 if (lpcone->n<1) return 0;
280 mu*=lpcone->muscale;
281 info=LPComputeATY(lpcone,Y,S);DSDPCHKERR(info);
282
283 info=DSDPVecSet(1.0,WX);
284 info=DSDPVecPointwiseDivide(WX,S,WX);DSDPCHKERR(info);
285
286 info=LPComputeATY(lpcone,DY,DS);DSDPCHKERR(info);
287 info=DSDPVecPointwiseMult(WX,DS,X);DSDPCHKERR(info);
288
289 info=DSDPVecScale(-mu,WX);DSDPCHKERR(info);
290 info=DSDPVecPointwiseMult(WX,X,X);DSDPCHKERR(info);
291 info=DSDPVecAXPY(-1.0,WX,X);DSDPCHKERR(info);
292 info=DSDPVecGetArray(X,&xx);DSDPCHKERR(info);
293 for (info=0;info<lpcone->n;info++){
294 if (xx[info]<0) xx[info]=0;
295 }
296 info=DSDPVecRestoreArray(X,&xx);DSDPCHKERR(info);
297 info=LPComputeAX(lpcone,X,WY);DSDPCHKERR(info);
298 info=DSDPVecAXPY(1.0,WY,AX);DSDPCHKERR(info);
299 info=DSDPVecDot(S,X,&dtracexs);DSDPCHKERR(info);
300 *tracexs+=dtracexs;
301 info=DSDPVecGetArray(X,&xx);DSDPCHKERR(info);
302 if (xout){
303 for (info=0;info<lpcone->n;info++){
304 if (xout){ xout[info]=xx[info]; }
305 }
306 }
307 info=DSDPVecRestoreArray(X,&xx);DSDPCHKERR(info);
308 DSDPFunctionReturn(0);
309}
310
311
312#undef __FUNCT__
313#define __FUNCT__ "LPConeS"
314static int LPConeS(void* dcone, DSDPVec Y, DSDPDualFactorMatrix flag,
315 DSDPTruth *psdefinite){
316 int i,n,info;
317 double *s;
318 LPCone lpcone=(LPCone)dcone;
319 DSDPVec S;
320
321 DSDPFunctionBegin;
322
323 if (lpcone->n<1) return 0;
324 if (flag==DUAL_FACTOR){
325 S=lpcone->DS;
326 } else {
327 S=lpcone->PS;
328 }
329
330 info=DSDPVecCopy(Y,lpcone->Y);DSDPCHKERR(info);
331 info=LPComputeATY(lpcone,Y,S);DSDPCHKERR(info);
332 info=DSDPVecGetC(Y,&lpcone->sscale);
333 info=DSDPVecGetSize(S,&n);DSDPCHKERR(info);
334 info=DSDPVecGetArray(S,&s);DSDPCHKERR(info);
335 *psdefinite=DSDP_TRUE;
336 for (i=0;i<n;i++){ if (s[i]<=0) *psdefinite=DSDP_FALSE;}
337 info=DSDPVecRestoreArray(S,&s);DSDPCHKERR(info);
338
339 DSDPFunctionReturn(0);
340}
341#undef __FUNCT__
342#define __FUNCT__ "LPConeInvertS"
343static int LPConeInvertS(void* dcone){
344 DSDPFunctionBegin;
345 DSDPFunctionReturn(0);
346}
347
348#undef __FUNCT__
349#define __FUNCT__ "LPConeComputeMaxStepLength"
350static int LPConeComputeMaxStepLength(void* dcone, DSDPVec DY, DSDPDualFactorMatrix flag, double *maxsteplength){
351 int i,n,info;
352 double *s,*ds,mstep=1.0e200;
353 LPCone lpcone=(LPCone)dcone;
354 DSDPVec S,DS=lpcone->WX;
355 DSDPFunctionBegin;
356 if (lpcone->n<1) return 0;
357 if (flag==DUAL_FACTOR){
358 S=lpcone->DS;
359 } else {
360 S=lpcone->PS;
361 }
362
363 info=LPComputeATY(lpcone,DY,DS);DSDPCHKERR(info);
364
365 info=DSDPVecGetSize(DS,&n);DSDPCHKERR(info);
366 info=DSDPVecGetArray(S,&s);DSDPCHKERR(info);
367 info=DSDPVecGetArray(DS,&ds);DSDPCHKERR(info);
368 for (i=0;i<n;i++) if (ds[i]<0){mstep=DSDPMin(mstep,-s[i]/ds[i]);}
369 info=DSDPVecRestoreArray(S,&s);DSDPCHKERR(info);
370 info=DSDPVecRestoreArray(DS,&ds);DSDPCHKERR(info);
371
372 *maxsteplength=mstep;
373
374 DSDPFunctionReturn(0);
375}
376
377
378#undef __FUNCT__
379#define __FUNCT__ "LPConePotential"
380static int LPConePotential(void* dcone, double *logobj, double *logdet){
381 int i,n,info;
382 double *s,mu,sumlog=0;
383 LPCone lpcone=(LPCone)dcone;
384 DSDPFunctionBegin;
385 if (lpcone->n<1) return 0;
386 mu=lpcone->muscale;
387 info=DSDPVecGetArray(lpcone->DS,&s);DSDPCHKERR(info);
388 info=DSDPVecGetSize(lpcone->DS,&n);DSDPCHKERR(info);
389 for (i=0;i<n;i++){
390 sumlog+= mu*log(s[i]);
391 }
392 info=DSDPVecRestoreArray(lpcone->DS,&s);DSDPCHKERR(info);
393 *logdet=sumlog;
394 *logobj=0;
395 DSDPFunctionReturn(0);
396}
397
398#undef __FUNCT__
399#define __FUNCT__ "LPConeSparsity"
400static int LPConeSparsity(void *dcone,int row, int *tnnz, int rnnz[], int m){
401 int info,*wi,n;
402 double *wd;
403 LPCone lpcone=(LPCone)dcone;
404 DSDPVec W=lpcone->WX;
405 DSDPFunctionBegin;
406 if (lpcone->n<1) return 0;
407 if (row==0) return 0;
408 if (row==m-1){
409 return 0;
410 }
411 info=DSDPVecGetSize(W,&n);DSDPCHKERR(info);
412 info=DSDPVecGetArray(W,&wd);DSDPCHKERR(info);
413 wi=(int*)wd;
414 info=SpRowMatRowNnz(lpcone->A,row-1,wi,n);DSDPCHKERR(info);
415 info=SpRowIMultAdd(lpcone->A,wi,n,rnnz+1,m-2);DSDPCHKERR(info);
416 info=DSDPVecRestoreArray(W,&wd);DSDPCHKERR(info);
417 DSDPFunctionReturn(0);
418}
419
420
421#undef __FUNCT__
422#define __FUNCT__ "LPConeMonitor"
423static int LPConeMonitor( void *dcone, int tag){
424 DSDPFunctionBegin;
425 DSDPFunctionReturn(0);
426}
427
428#undef __FUNCT__
429#define __FUNCT__ "LPANorm2"
430static int LPANorm2( void *dcone, DSDPVec ANorm){
431 int i,info;
432 double dd;
433 LPCone lpcone=(LPCone)dcone;
434 DSDPFunctionBegin;
435 if (lpcone->n<1) return 0;
436 info=DSDPVecNorm22(lpcone->C,&dd);DSDPCHKERR(info);
437 info=DSDPVecAddC(ANorm,dd);DSDPCHKERR(info);
438 for (i=0;i<lpcone->m;i++){
439 info=SpRowMatNorm2(lpcone->A,i,&dd);DSDPCHKERR(info);
440 info=DSDPVecAddElement(ANorm,i+1,dd);DSDPCHKERR(info);
441 }
442 info=DSDPVecAddR(ANorm,1.0);DSDPCHKERR(info);
443 DSDPFunctionReturn(0);
444}
445
446
447static struct DSDPCone_Ops kops;
448static const char *lpconename="LP Cone";
449
450#undef __FUNCT__
451#define __FUNCT__ "LPConeOperationsInitialize"
452static int LPConeOperationsInitialize(struct DSDPCone_Ops* coneops){
453 int info;
454 if (coneops==NULL) return 0;
455 info=DSDPConeOpsInitialize(coneops); DSDPCHKERR(info);
456 coneops->conehessian=LPConeHessian;
457 coneops->conerhs=LPConeRHS;
458 coneops->conesetup=LPConeSetup;
459 coneops->conesetup2=LPConeSetup2;
460 coneops->conedestroy=LPConeDestroy;
461 coneops->conecomputes=LPConeS;
462 coneops->coneinverts=LPConeInvertS;
463 coneops->conesetxmaker=LPConeSetX;
464 coneops->conecomputex=LPConeX;
465 coneops->conemaxsteplength=LPConeComputeMaxStepLength;
466 coneops->conelogpotential=LPConePotential;
467 coneops->conesize=LPConeSize;
468 coneops->conesparsity=LPConeSparsity;
469 coneops->conehmultiplyadd=LPConeMultiply;
470 coneops->conemonitor=LPConeMonitor;
471 coneops->coneanorm2=LPANorm2;
472 coneops->id=2;
473 coneops->name=lpconename;
474 return 0;
475}
476
477#undef __FUNCT__
478#define __FUNCT__ "DSDPAddLP"
479int DSDPAddLP(DSDP dsdp,LPCone lpcone){
480 int info;
481 DSDPFunctionBegin;
482 info=LPConeOperationsInitialize(&kops); DSDPCHKERR(info);
483 info=DSDPAddCone(dsdp,&kops,(void*)lpcone); DSDPCHKERR(info);
484 DSDPFunctionReturn(0);
485}
486
487#undef __FUNCT__
488#define __FUNCT__ "DSDPCreateLPCone"
509int DSDPCreateLPCone(DSDP dsdp, LPCone *dspcone){
510 int m,info;
511 struct LPCone_C *lpcone;
512 DSDPFunctionBegin;
513 DSDPCALLOC1(&lpcone,struct LPCone_C,&info);DSDPCHKERR(info);
514 *dspcone=lpcone;
515 /*
516 info=DSDPAddLP(dsdp,lpcone);DSDPCHKERR(info);
517 */
518 info=LPConeOperationsInitialize(&kops); DSDPCHKERR(info);
519 info=DSDPAddCone(dsdp,&kops,(void*)lpcone); DSDPCHKERR(info);
520 info=DSDPGetNumberOfVariables(dsdp,&m);DSDPCHKERR(info);
521 lpcone->m=m;
522 lpcone->muscale=1.0;
523 lpcone->n=0;
524 lpcone->xout=0;
525 lpcone->r=1.0;
526 info=DSDPVecCreateSeq(0,&lpcone->C);DSDPCHKERR(info);
527 info=DSDPVecCreateSeq(0,&lpcone->WY);DSDPCHKERR(info);
528 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX);DSDPCHKERR(info);
529 info=DSDPVecDuplicate(lpcone->C,&lpcone->WX2);DSDPCHKERR(info);
530 info=DSDPVecDuplicate(lpcone->C,&lpcone->PS);DSDPCHKERR(info);
531 info=DSDPVecDuplicate(lpcone->C,&lpcone->DS);DSDPCHKERR(info);
532 info=DSDPVecDuplicate(lpcone->C,&lpcone->X);DSDPCHKERR(info);
533 DSDPFunctionReturn(0);
534}
535
536
537#undef __FUNCT__
538#define __FUNCT__ "LPConeGetXArray"
556int LPConeGetXArray(LPCone lpcone,double *x[], int *n){
557 int info;
558 DSDPFunctionBegin;
559 info=DSDPVecGetArray(lpcone->X,x);DSDPCHKERR(info);
560 info=DSDPVecGetSize(lpcone->X,n);DSDPCHKERR(info);
561 DSDPFunctionReturn(0);
562}
563
564#undef __FUNCT__
565#define __FUNCT__ "LPConeGetSArray"
566/*
567\fn int LPConeGetSArray(LPCone lpcone, double *s[], int *n);
568\brief Get the array used to store the s variables
569\ingroup LPRoutines
570\param lpcone LP Cone
571\param s array of variables
572\param n the dimension of the cone and length of the array
573\sa DSDPCreateLPCone()
574 */
575int LPConeGetSArray(LPCone lpcone,double *s[], int *n){
576 int info;
577 DSDPFunctionBegin;
578 info=DSDPVecGetArray(lpcone->DS,s);DSDPCHKERR(info);
579 info=DSDPVecGetSize(lpcone->DS,n);DSDPCHKERR(info);
580 DSDPFunctionReturn(0);
581}
582
583#undef __FUNCT__
584#define __FUNCT__ "LPConeCopyS"
595int LPConeCopyS(LPCone lpcone,double s[], int n){
596 int i,info;
597 double *ss,sscale=lpcone->sscale;
598 DSDPTruth psdefinite;
599 DSDPFunctionBegin;
600 info=LPConeS((void*)lpcone,lpcone->Y,DUAL_FACTOR ,&psdefinite);DSDPCHKERR(info);
601 info=DSDPVecGetArray(lpcone->DS,&ss);DSDPCHKERR(info);
602 for (i=0;i<n;i++) s[i]=ss[i]/fabs(sscale);
603 DSDPFunctionReturn(0);
604}
605
606#undef __FUNCT__
607#define __FUNCT__ "LPConeGetDimension"
616int LPConeGetDimension(LPCone lpcone,int *n){
617 DSDPFunctionBegin;
618 *n=(int)(lpcone->n*lpcone->muscale);
619 DSDPFunctionReturn(0);
620}
621
622
623#undef __FUNCT__
624#define __FUNCT__ "LPConeScaleBarrier"
625int LPConeScaleBarrier(LPCone lpcone,double muscale){
626 DSDPFunctionBegin;
627 if (muscale>0){
628 lpcone->muscale=muscale;
629 }
630 DSDPFunctionReturn(0);
631}
632
633#undef __FUNCT__
634#define __FUNCT__ "LPConeSetData"
666int LPConeSetData(LPCone lpcone,int n, const int ik[],const int cols[],const double vals[]){
667 int info,i,spot,m=lpcone->m;
668 DSDPVec C;
669 DSDPFunctionBegin;
670 lpcone->n=n;
671 info=DSDPVecCreateSeq(n,&C);DSDPCHKERR(info);
672 lpcone->C=C;
673 info=DSDPVecZero(C);DSDPCHKERR(info);
674 lpcone->muscale=1.0;
675 if (n<100) lpcone->muscale=1.0;
676 if (n<10) lpcone->muscale=1.0;
677 for (i=ik[0];i<ik[1];i++){
678 info=DSDPVecSetElement(C,cols[i],vals[i]);
679 }
680 spot=ik[0];
681 info=CreateSpRowMatWdata(m,n,vals+spot,cols+spot,ik+1,&lpcone->A);DSDPCHKERR(info);
682 DSDPFunctionReturn(0);
683}
684
685#undef __FUNCT__
686#define __FUNCT__ "LPConeSetData2"
717int LPConeSetData2(LPCone lpcone,int n, const int ik[],const int cols[],const double vals[]){
718 int info,i,spot,m=lpcone->m;
719 DSDPVec C;
720 DSDPFunctionBegin;
721 lpcone->n=n;
722 info=DSDPVecCreateSeq(n,&C);DSDPCHKERR(info);
723 lpcone->C=C;
724 info=DSDPVecZero(C);DSDPCHKERR(info);
725 lpcone->muscale=1.0;
726 if (n<100) lpcone->muscale=1.0;
727 if (n<10) lpcone->muscale=1.0;
728 for (i=ik[m];i<ik[m+1];i++){
729 info=DSDPVecSetElement(C,cols[i],vals[i]);
730 }
731 spot=ik[0];
732 info=CreateSpRowMatWdata(m,n,vals+spot,cols+spot,ik,&lpcone->A);DSDPCHKERR(info);
733 DSDPFunctionReturn(0);
734}
735
736#undef __FUNCT__
737#define __FUNCT__ "LPConeView2"
744int LPConeView2(LPCone lpcone){
745 int info;
746 DSDPFunctionBegin;
747 printf("LPCone Constraint Matrix\n");
748 info=SpRowMatView(lpcone->A);DSDPCHKERR(info);
749 printf("LPCone Objective C vector\n");
750 info=DSDPVecView(lpcone->C);DSDPCHKERR(info);
751 DSDPFunctionReturn(0);
752}
753
754
755
756#undef __FUNCT__
757#define __FUNCT__ "LPConeGetConstraint"
758int LPConeGetConstraint(LPCone lpcone,int column,DSDPVec W){
759 int n,info;
760 double *w;
761 DSDPFunctionBegin;
762 if (column==0){
763 info=DSDPVecCopy(lpcone->C,W);DSDPCHKERR(info);
764 } else {
765 info=DSDPVecGetSize(W,&n);DSDPCHKERR(info);
766 info=DSDPVecGetArray(W,&w);DSDPCHKERR(info);
767 info=SpRowMatGetRowVector(lpcone->A,column-1,w,n);info=0;DSDPCHKERR(info);
768 info=DSDPVecRestoreArray(W,&w);DSDPCHKERR(info);
769 }
770 DSDPFunctionReturn(0);
771}
772
773#undef __FUNCT__
774#define __FUNCT__ "LPConeGetData"
783int LPConeGetData(LPCone lpcone,int vari,double vv[], int n){
784 int info;
785 DSDPVec W;
786 DSDPFunctionBegin;
787 info=DSDPVecCreateWArray(&W,vv,n);DSDPCHKERR(info);
788 info=LPConeGetConstraint(lpcone,vari,W);DSDPCHKERR(info);
789 DSDPFunctionReturn(0);
790}
791
792#undef __FUNCT__
793#define __FUNCT__ "LPConeSetXVec"
794int LPConeSetXVec(LPCone lpcone,double *xout, int n){
795 DSDPFunctionBegin;
796 if (n==lpcone->n) lpcone->xout=xout;
797 DSDPFunctionReturn(0);
798}
799
800
801static int vsdot(const int*,const double *,int,const double *, int, double *);
802static int checkvsparse(smatx *);
803
804
805static int CreateSpRowMatWdata(int m,int n,const double vals[], const int cols[], const int ik[],
806 smatx **A){
807
808 smatx *V;
809
810 V=(smatx*) malloc(1*sizeof(smatx));
811 if (V==NULL) return 1;
812 V->nrow=m;
813 V->ncol=n;
814 V->owndata=0;
815 V->an=vals; V->col=cols; V->nnz=ik;
816
817 *A=V;
818 checkvsparse(V);
819 return 0;
820}
821
822static int vsdot(const int ja[], const double an[], int nn0, const double vv[], int n, double *vdot){
823
824 int i;
825 double tmp=0.0;
826
827 for (i=0; i<nn0; i++){
828 /* if (ja[i]<n) tmp = tmp + an[i] * vv[ja[i]]; */
829 tmp += an[i] * vv[ja[i]];
830 }
831 *vdot=tmp;
832 return 0;
833}
834
835static int checkvsparse(smatx *A){
836 int i,k=0,m=A->nrow,tnnz=0;
837 const int *nnz=A->nnz;
838
839 for (i=0;i<m;++i){
840 if (nnz[i+1]-nnz[i]>0){
841 tnnz++;
842 }
843 }
844 if (tnnz < m/2){
845 A->nzrows =(int*)malloc((tnnz)*sizeof(int));
846 A->nnzrows=tnnz;
847 for (i=0;i<m;++i){
848 if (nnz[i+1]-nnz[i]>0){
849 A->nzrows[k]=i;
850 k++;
851 }
852 }
853 } else {
854 A->nzrows = 0;
855 A->nnzrows = m;
856 }
857 return 0;
858}
859
860#undef __FUNCT__
861#define __FUNCT__ "SpRowMatMult"
862static int SpRowMatMult(smatx* A, const double x[], int n, double y[], int m){
863
864 int i,k1,k2,nrow=A->nrow;
865 const int *nnz=A->nnz,*col=A->col;
866 const double *an=A->an;
867
868 if (A->ncol != n) return 1;
869 if (A->nrow != m) return 2;
870 if (x==0 && n>0) return 3;
871 if (y==0 && m>0) return 3;
872
873 if (m>0){
874 memset((void*)y,0,m*sizeof(double));
875 for (i=0; i<nrow; i++){
876 k1=*(nnz+i); k2=*(nnz+i+1);
877 vsdot(col+k1,an+k1,k2-k1,x,n,y+i);
878 }
879 }
880 return 0;
881}
882
883#undef __FUNCT__
884#define __FUNCT__ "SpRowMatMultTrans"
885static int SpRowMatMultTrans(smatx * A, const double x[], int m, double y[], int n){
886
887 int i,j,k1,k2,nrow=A->nrow;
888 const int *col=A->col,*nnz=A->nnz;
889 const double *an=A->an;
890 if (A->ncol != n) return 1;
891 if (A->nrow != m) return 2;
892 if (x==0 && m>0) return 3;
893 if (y==0 && n>0) return 3;
894
895 memset((void*)y,0,n*sizeof(double));
896 for (i=0; i<nrow; i++){
897 k1=nnz[i]; k2=nnz[i+1];
898 for (j=k1; j<k2; j++){
899 y[col[j]] += an[j]*x[i];
900 }
901 }
902
903 return 0;
904}
905
906
907static int SpRowMatRowNnz(smatx *A, int row, int* wi,int n){
908 int i,k1,k2;
909 const int *col=A->col;
910 DSDPFunctionBegin;
911 memset((void*)wi,0,n*sizeof(double));
912 k1=A->nnz[row];
913 k2=A->nnz[row+1];
914 for (i=k1; i<k2; i++){
915 wi[col[i]]=1;
916 }
917 DSDPFunctionReturn(0);
918}
919
920static int SpRowIMultAdd(smatx *A,int *wi,int n,int *rnnz,int m){
921
922 int i,j,k1,k2,nrow=A->nrow;
923 const int *nnz=A->nnz,*col=A->col;
924 DSDPFunctionBegin;
925 for (i=0; i<nrow; i++){
926 k1=nnz[i];
927 k2=nnz[i+1];
928 for (j=k1; j<k2; j++){
929 if (wi[col[j]]){
930 rnnz[i]++;
931 }
932 }
933 }
934 DSDPFunctionReturn(0);
935}
936/*
937static int SpRowMatAddRowMultiple(smatx* A, int nrow, double ytmp, double row[], int n){
938 int k;
939 int *col=A->col, *nnz=A->nnz;
940 double *an=A->an;
941
942 for (k=nnz[nrow]; k<nnz[nrow+1]; k++){
943 row[col[k]] += ytmp * an[k];
944 }
945
946 return 0;
947}
948*/
949static int SpRowMatNorm2(smatx* A, int nrow, double *norm22){
950 int k;
951 const int *nnz=A->nnz;
952 double tt=0;
953 const double *an=A->an;
954
955 for (k=nnz[nrow]; k<nnz[nrow+1]; k++){
956 tt+=an[k]*an[k];
957 }
958 *norm22=tt;
959 return 0;
960}
961
962
963#undef __FUNCT__
964#define __FUNCT__ "SpRowMatGetRowVector"
965static int SpRowMatGetRowVector(smatx* M, int row, double r[], int m){
966
967 int i,k1,k2;
968 const int *col=M->col;
969 const double *an=M->an;
970 /*
971 if (M->ncol != m) return 1;
972 if (row<0 || row>= M->nrow) return 2;
973 if (r==0) return 3;
974 */
975 memset((void*)r,0,m*sizeof(double));
976 k1=M->nnz[row];
977 k2=M->nnz[row+1];
978
979 for (i=k1; i<k2; i++){
980 r[col[i]]=an[i];
981 }
982
983 return 0;
984}
985#undef __FUNCT__
986#define __FUNCT__ "SpRowMatGetScaledRowVector"
987static int SpRowMatGetScaledRowVector(smatx* M, int row, const double ss[], double r[], int m){
988
989 int i,k1,k2;
990 const int *col=M->col;
991 const double *an=M->an;
992 /*
993 if (M->ncol != m) return 1;
994 if (row<0 || row>= M->nrow) return 2;
995 if (r==0) return 3;
996 */
997 memset((void*)r,0,m*sizeof(double));
998 k1=M->nnz[row];
999 k2=M->nnz[row+1];
1000
1001 for (i=k1; i<k2; i++){
1002 r[col[i]]=ss[col[i]]*an[i];
1003 }
1004
1005 return 0;
1006}
1007
1008
1009/*
1010#undef __FUNCT__
1011#define __FUNCT__ "SpRowMatZero"
1012static int SpRowMatZero(smatx* M){
1013
1014 int nnz=M->nnz[M->nrow];
1015 memset(M->an,0,(nnz)*sizeof(double));
1016
1017 return 0;
1018}
1019
1020
1021
1022#undef __FUNCT__
1023#define __FUNCT__ "SpRowGetSize"
1024static int SpRowMatGetSize(smatx *M, int *m, int *n){
1025 *m=M->nrow;
1026 *n=M->ncol;
1027
1028 return 0;
1029}
1030*/
1031#undef __FUNCT__
1032#define __FUNCT__ "SpRowMatDestroy"
1033static int SpRowMatDestroy(smatx* A){
1034
1035 if (A->owndata){
1036 printf("Can't free array");
1037 return 1;
1038 /*
1039 if (A->an) free(A->an);
1040 if (A->col) free(A->col);
1041 if (A->nnz) free(A->nnz);
1042 */
1043 }
1044 if (A->nzrows) free(A->nzrows);
1045 free(A);
1046
1047 return 0;
1048}
1049
1050#undef __FUNCT__
1051#define __FUNCT__ "SpRowMatView"
1052static int SpRowMatView(smatx* M){
1053
1054 int i,j,k1,k2;
1055
1056 for (i=0; i<M->nrow; i++){
1057 k1=M->nnz[i]; k2=M->nnz[i+1];
1058
1059 if (k2-k1 >0){
1060 printf("Row %d, (Variable y%d) : ",i,i+1);
1061 for (j=k1; j<k2; j++)
1062 printf(" %4.2e x%d + ",M->an[j],M->col[j]);
1063 printf("= dobj%d \n",i+1);
1064 }
1065 }
1066
1067 return 0;
1068}
1069
1070#undef __FUNCT__
1071#define __FUNCT__ "LPConeView"
1078int LPConeView(LPCone lpcone){
1079
1080 smatx* A=lpcone->A;
1081 int i,j,jj,info;
1082 const int *row=A->col,*nnz=A->nnz;
1083 int n=A->ncol,m=A->nrow;
1084 const double *an=A->an;
1085 double cc;
1086 DSDPVec C=lpcone->C;
1087 DSDPFunctionBegin;
1088 printf("LPCone Constraint Matrix\n");
1089 printf("Number y variables 1 through %d\n",m);
1090 for (i=0; i<n; i++){
1091 printf("Inequality %d: ",i);
1092 for (j=0;j<m;j++){
1093 for (jj=nnz[j];jj<nnz[j+1];jj++){
1094 if (row[jj]==i){
1095 printf("%4.2e y%d + ",an[jj],j+1);
1096 }
1097 }
1098 }
1099 info=DSDPVecGetElement(C,i,&cc);DSDPCHKERR(info);
1100 printf(" <= %4.2e\n",cc);
1101 }
1102 DSDPFunctionReturn(0);
1103}
1104
1105
The API to DSDP for those applications using DSDP as a subroutine library.
struct LPCone_C * LPCone
The LPCone object points to blocks of data that specify linear scalar inequality constraints.
Definition dsdp5.h:27
DSDPTruth
Boolean variables.
@ DSDP_FALSE
@ DSDP_TRUE
DSDPDualFactorMatrix
DSDP requires two instances of the data structures S.
@ DUAL_FACTOR
struct DSDP_C * DSDP
An implementation of the dual-scaling algorithm for semidefinite programming.
int DSDPConeOpsInitialize(struct DSDPCone_Ops *dops)
Initialize the function pointers to 0.
Definition dsdpcone.c:443
Implementations of a cone (SDP,LP,...) must provide a structure of function pointers.
int DSDPAddCone(DSDP, struct DSDPCone_Ops *, void *)
Apply DSDP to a conic structure.
Definition dsdpcops.c:569
int LPConeGetData(LPCone lpcone, int vari, double vv[], int n)
Get one column (or row) of the LP data.
Definition dsdplp.c:783
int DSDPSchurMatRowColumnScaling(DSDPSchurMat, int, DSDPVec, int *)
Get the scaling and nonzero pattern of each column in this row of the matrix.
int DSDPSchurMatAddRow(DSDPSchurMat, int, double, DSDPVec)
Add elements to a row of the Schur matrix.
struct DSDPSchurMat_C DSDPSchurMat
This object represents the Schur Matrix. Its structure is opaque to the DSDP solver,...
int DSDPSchurMatDiagonalScaling(DSDPSchurMat, DSDPVec)
Get the scaling and nonzero pattern of each diagonal element of the matrix.
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 DSDPGetNumberOfVariables(DSDP dsdp, int *m)
Copy the number of variables y.
int DSDPCreateLPCone(DSDP dsdp, LPCone *dspcone)
Create a new object for linear programs and scalar inequalities.
Definition dsdplp.c:509
int LPConeSetData2(LPCone lpcone, int n, const int ik[], const int cols[], const double vals[])
Set data A and into the LP cone.
Definition dsdplp.c:717
int LPConeGetXArray(LPCone lpcone, double *x[], int *n)
Get the array used to store the x variables.
Definition dsdplp.c:556
int LPConeSetData(LPCone lpcone, int n, const int ik[], const int cols[], const double vals[])
Set data into the LP cone.
Definition dsdplp.c:666
int LPConeView(LPCone lpcone)
Print the data in the LP cone to the screen.
Definition dsdplp.c:1078
int LPConeCopyS(LPCone lpcone, double s[], int n)
Copy the variables s into the spedified array.
Definition dsdplp.c:595
int LPConeGetDimension(LPCone lpcone, int *n)
Get the dimension is the number of variables x, which equals the number of slack variables s.
Definition dsdplp.c:616
int LPConeView2(LPCone lpcone)
Print the data in the LP cone to the screen.
Definition dsdplp.c:744