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 
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
DSDPTruth
Boolean variables.
int LPConeView2(LPCone lpcone)
Print the data in the LP cone to the screen.
Definition: dsdplp.c:744
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.
Definition: dsdpschurmat.h:35
Error handling, printing, and profiling.
Internal structures for the DSDP solver.
Definition: dsdp.h:65
int DSDPCreateLPCone(DSDP dsdp, LPCone *dspcone)
Create a new object for linear programs and scalar inequalities.
Definition: dsdplp.c:509
int DSDPSchurMatRowColumnScaling(DSDPSchurMat, int, DSDPVec, int *)
Get the scaling and nonzero pattern of each column in this row of the matrix.
int LPConeCopyS(LPCone lpcone, double s[], int n)
Copy the variables s into the spedified array.
Definition: dsdplp.c:595
int LPConeGetXArray(LPCone lpcone, double *x[], int *n)
Get the array used to store the x variables.
Definition: dsdplp.c:556
The API to DSDP for those applications using DSDP as a subroutine library.
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 LPConeGetData(LPCone lpcone, int vari, double vv[], int n)
Get one column (or row) of the LP data.
Definition: dsdplp.c:783
int DSDPSchurMatDiagonalScaling(DSDPSchurMat, DSDPVec)
Get the scaling and nonzero pattern of each diagonal element of the matrix.
int LPConeView(LPCone lpcone)
Print the data in the LP cone to the screen.
Definition: dsdplp.c:1078
int DSDPSchurMatAddRow(DSDPSchurMat, int, double, DSDPVec)
Add elements to a row of the Schur matrix.
Implementations of a cone (SDP,LP,...) must provide a structure of function pointers.
DSDPDualFactorMatrix
DSDP requires two instances of the data structures S.
int DSDPGetNumberOfVariables(DSDP dsdp, int *m)
Copy the number of variables y.
Definition: dsdpsetdata.c:707
int DSDPAddCone(DSDP, struct DSDPCone_Ops *, void *)
Apply DSDP to a conic structure.
Definition: dsdpcops.c:569
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 DSDPConeOpsInitialize(struct DSDPCone_Ops *dops)
Initialize the function pointers to 0.
Definition: dsdpcone.c:443
struct LPCone_C * LPCone
The LPCone object points to blocks of data that specify linear scalar inequality constraints.
Definition: dsdp5.h:27