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>
14 typedef 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 
31 struct 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 
44 static int CreateSpRowMatWdata(int,int,const double[], const int[], const int[],smatx **);
45 static int SpRowMatMult(smatx*,const double[], int , double[], int);
46 static int SpRowMatMultTrans(smatx *, const double[],int, double[],int);
47 static int SpRowMatGetRowVector(smatx*, int, double*,int);
48 static int SpRowMatGetScaledRowVector(smatx*, int, const double[], double*, int);
49 static int SpRowMatDestroy(smatx*);
50 static int SpRowMatView(smatx*);
51 /*
52 static int SpRowMatGetSize(smatx *, int *, int *);
53 static int SpRowMatZero(smatx*);
54 static int SpRowMatAddRowMultiple(smatx*, int, double, double[], int);
55 */
56 static int SpRowIMultAdd(smatx*,int*,int,int *,int);
57 static int SpRowMatRowNnz(smatx*, int, int*,int);
58 static int SpRowMatNorm2(smatx*, int, double*);
59 
60 
61 #undef __FUNCT__
62 #define __FUNCT__ "LPConeSetUp"
63 static 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"
82 static 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"
92 static 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"
113 static 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"
124 static 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"
145 static 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"
167 static 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"
221 static 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"
244 static 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"
263 static int LPConeSetX(void* dcone,double mu, DSDPVec Y,DSDPVec DY){
264  DSDPFunctionBegin;
265  DSDPFunctionReturn(0);
266 }
267 
268 #undef __FUNCT__
269 #define __FUNCT__ "LPConeX"
270 static 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"
314 static 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"
343 static int LPConeInvertS(void* dcone){
344  DSDPFunctionBegin;
345  DSDPFunctionReturn(0);
346 }
347 
348 #undef __FUNCT__
349 #define __FUNCT__ "LPConeComputeMaxStepLength"
350 static 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"
380 static 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"
400 static 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"
423 static int LPConeMonitor( void *dcone, int tag){
424  DSDPFunctionBegin;
425  DSDPFunctionReturn(0);
426 }
427 
428 #undef __FUNCT__
429 #define __FUNCT__ "LPANorm2"
430 static 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 
447 static struct DSDPCone_Ops kops;
448 static const char *lpconename="LP Cone";
449 
450 #undef __FUNCT__
451 #define __FUNCT__ "LPConeOperationsInitialize"
452 static 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"
479 int 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"
489 
509 int 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"
539 
556 int 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  */
575 int 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"
585 
595 int 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"
608 
616 int 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"
625 int 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"
635 
666 int 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"
687 
717 int 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"
738 
744 int 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"
758 int 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"
775 
783 int 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"
794 int LPConeSetXVec(LPCone lpcone,double *xout, int n){
795  DSDPFunctionBegin;
796  if (n==lpcone->n) lpcone->xout=xout;
797  DSDPFunctionReturn(0);
798 }
799 
800 
801 static int vsdot(const int*,const double *,int,const double *, int, double *);
802 static int checkvsparse(smatx *);
803 
804 
805 static 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 
822 static 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 
835 static 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"
862 static 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"
885 static 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 
907 static 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 
920 static 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 /*
937 static 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 */
949 static 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"
965 static 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"
987 static 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"
1012 static 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"
1024 static 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"
1033 static 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"
1052 static 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"
1072 
1078 int 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 
DSDP_FALSE
@ DSDP_FALSE
Definition: dsdpbasictypes.h:19
DSDPCreateLPCone
int DSDPCreateLPCone(DSDP dsdp, LPCone *dspcone)
Create a new object for linear programs and scalar inequalities.
Definition: dsdplp.c:509
LPConeSetData2
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
DSDPGetNumberOfVariables
int DSDPGetNumberOfVariables(DSDP dsdp, int *m)
Copy the number of variables y.
Definition: dsdpsetdata.c:707
DUAL_FACTOR
@ DUAL_FACTOR
Definition: dsdpbasictypes.h:26
DSDPSchurMatRowColumnScaling
int DSDPSchurMatRowColumnScaling(DSDPSchurMat, int, DSDPVec, int *)
Get the scaling and nonzero pattern of each column in this row of the matrix.
Definition: dsdpschurmatadd.c:33
DSDPVec
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
LPConeGetData
int LPConeGetData(LPCone lpcone, int vari, double vv[], int n)
Get one column (or row) of the LP data.
Definition: dsdplp.c:783
DSDPTruth
DSDPTruth
Boolean variables.
Definition: dsdpbasictypes.h:19
DSDPSchurMatDiagonalScaling
int DSDPSchurMatDiagonalScaling(DSDPSchurMat, DSDPVec)
Get the scaling and nonzero pattern of each diagonal element of the matrix.
Definition: dsdpschurmatadd.c:235
DSDP_C
Internal structures for the DSDP solver.
Definition: dsdp.h:65
LPCone
struct LPCone_C * LPCone
The LPCone object points to blocks of data that specify linear scalar inequality constraints.
Definition: dsdp5.h:27
LPConeGetDimension
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
LPConeView
int LPConeView(LPCone lpcone)
Print the data in the LP cone to the screen.
Definition: dsdplp.c:1078
dsdpsys.h
Error handling, printing, and profiling.
DSDPConeOpsInitialize
int DSDPConeOpsInitialize(struct DSDPCone_Ops *dops)
Initialize the function pointers to 0.
Definition: dsdpcone.c:443
DSDPSchurMatAddRow
int DSDPSchurMatAddRow(DSDPSchurMat, int, double, DSDPVec)
Add elements to a row of the Schur matrix.
Definition: dsdpschurmatadd.c:76
LPConeGetXArray
int LPConeGetXArray(LPCone lpcone, double *x[], int *n)
Get the array used to store the x variables.
Definition: dsdplp.c:556
dsdp5.h
The API to DSDP for those applications using DSDP as a subroutine library.
LPConeView2
int LPConeView2(LPCone lpcone)
Print the data in the LP cone to the screen.
Definition: dsdplp.c:744
LPConeSetData
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
dsdpcone_impl.h
Implementations of a cone (SDP,LP,...) must provide a structure of function pointers.
DSDP_TRUE
@ DSDP_TRUE
Definition: dsdpbasictypes.h:19
DSDPSchurMat_C
Schur complement matrix whose solution is the Newton direction.
Definition: dsdpschurmat.h:35
LPConeCopyS
int LPConeCopyS(LPCone lpcone, double s[], int n)
Copy the variables s into the spedified array.
Definition: dsdplp.c:595
DSDPAddCone
int DSDPAddCone(DSDP, struct DSDPCone_Ops *, void *)
Apply DSDP to a conic structure.
Definition: dsdpcops.c:569
DSDPDualFactorMatrix
DSDPDualFactorMatrix
DSDP requires two instances of the data structures S.
Definition: dsdpbasictypes.h:25