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
DSDPDualFactorMatrix
DSDP requires two instances of the data structures S.
@ DUAL_FACTOR
DSDPTruth
Boolean variables.
@ DSDP_FALSE
@ DSDP_TRUE
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.
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
Schur complement matrix whose solution is the Newton direction.
Internal structures for the DSDP solver.
Definition dsdp.h:65