DSDP
sdpcompute.c
Go to the documentation of this file.
1 #include "dsdpsdp.h"
2 #include "dsdpsys.h"
8 /*
9 static int tt1=0,tt2=0;
10  DSDPEventLogBegin(tt1);
11  DSDPEventLogBegin(tt2);
12  DSDPEventLogEnd(tt2);
13  DSDPEventLogEnd(tt1);
14  if (tt1==0){DSDPEventLogRegister("SDP T1 Event",&tt1);}
15  if (tt2==0){DSDPEventLogRegister("SDP T2 Event",&tt2);}
16 */
17 #undef __FUNCT__
18 #define __FUNCT__ "SDPConeComputeHessian"
19 
30 int SDPConeComputeHessian( SDPCone sdpcone, double mu, DSDPSchurMat M, DSDPVec vrhs1, DSDPVec vrhs2){
31 
32  int i,k,kt,kk,m,n,rank,info;
33  int ncols,ii,idA;
34  double rtemp,ack,ggamma,bmu,scl;
35  double rhs1i,rhs2i;
36  DSDPTruth method1;
37  SDPConeVec W,W2;
38  DSDPVec MRowI=sdpcone->Work, Select=sdpcone->Work2;
39  DSDPDataMat AA;
40  DSDPDualMat S;
41  DSDPVMat T;
42  DSDPDataTranspose ATranspose=sdpcone->ATR;
43  SDPblk *blk=sdpcone->blk;
44  DSDPIndex IS;
45 
46  /* Evaluate M */
47  DSDPFunctionBegin;
48  SDPConeValid(sdpcone);
49  info=DSDPVecGetSize(vrhs1,&m);DSDPCHKERR(info);
50 
51  for (i=0; i<m; i++){ /* One row at a time */
52 
53  /* Which Coluns */
54  rhs1i=0;rhs2i=0;
55  info=DSDPVecZero(MRowI);DSDPCHKERR(info);
56  info=DSDPSchurMatRowColumnScaling(M,i,Select,&ncols); DSDPCHKERR(info);
57  if (ncols==0){continue;}
58 
59  for (kt=0; kt<ATranspose.nnzblocks[i]; kt++){ /* Loop over all blocks */
60  kk=ATranspose.nzblocks[i][kt];
61  idA=ATranspose.idA[i][kt];
62  info=DSDPBlockGetMatrix(&blk[kk].ADATA,idA,&ii,&scl,&AA);DSDPCHKBLOCKERR(kk,info);
63  if (ii!=i){DSDPSETERR2(8,"Data Transpose Error: var %d does not equal %d.\n",i,ii);}
64  info = DSDPDataMatGetRank(AA,&rank,blk[kk].n);DSDPCHKBLOCKERR(kk,info);
65  if (rank==0) continue;
66 
67  T=blk[kk].T; S=blk[kk].S; W=blk[kk].W; W2=blk[kk].W2;
68  n=blk[kk].n; ggamma=blk[kk].gammamu; bmu=blk[kk].bmu; IS=blk[kk].IS;
69 
70  method1=DSDP_TRUE; /* Simple heuristic */
71  if (rank==1) method1=DSDP_FALSE;
72  if (rank==2 && ncols<=n) method1=DSDP_FALSE;
73  if (rank*rank*ncols<=n+1)method1=DSDP_FALSE;
74  if (ncols*blk[kk].nnz < n*n/10) method1=DSDP_FALSE;
75  if (ncols==1 && i==m-1)method1=DSDP_FALSE;
76  if (n<5) method1=DSDP_TRUE;
77  if (0==1) method1=DSDP_FALSE;
78  if (method1==DSDP_TRUE){info=DSDPVMatZeroEntries(T);DSDPCHKBLOCKERR(kk,info);}
79  for (k=0; k<rank; k++){
80 
81  info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKBLOCKERR(kk,info);
82  if (ack==0.0) continue;
83  ack*=scl;
84  info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKBLOCKERR(kk,info);
85 
86  /* RHS terms */
87  info = SDPConeVecDot(W,W2,&rtemp); DSDPCHKBLOCKERR(kk,info);
88  if (rtemp==0.0) continue;
89  rhs1i+=rtemp*ack*bmu; rhs2i+=rtemp*ack*ggamma*mu;
90  ack*=(ggamma+bmu);
91 
92  if (method1==DSDP_TRUE){
93  info=DSDPVMatAddOuterProduct(T,ack*mu,W2);DSDPCHKBLOCKERR(kk,info);
94  } else {
95  info=DSDPBlockvAv(&blk[kk].ADATA,ack*mu,Select,W2,MRowI);DSDPCHKBLOCKERR(kk,info);
96  } /* End row computations for rank kk of block kk */
97 
98  } /* End row computations for all of block kk */
99 
100  if (method1==DSDP_TRUE){
101  info=DSDPBlockADot(&blk[kk].ADATA,1.0,Select,T,MRowI);DSDPCHKBLOCKERR(kk,info);
102  } /* End row computations for all of block ll */
103  } /* End row computations for all blocks */
104  info=DSDPVecAddElement(vrhs1,i,rhs1i);DSDPCHKERR(info);
105  info=DSDPVecAddElement(vrhs2,i,rhs2i);DSDPCHKERR(info);
106  info=DSDPSchurMatAddRow(M,i,1.0,MRowI);DSDPCHKERR(info);
107  }
108 
109  DSDPFunctionReturn(0);
110 }
111 
112 #undef __FUNCT__
113 #define __FUNCT__ "SDPConeComputeRHS"
114 
125 int SDPConeComputeRHS( SDPCone sdpcone, int blockj, double mu, DSDPVec vrow, DSDPVec vrhs1, DSDPVec vrhs2){
126 
127  int info,i,ii,k,rank,nnzmats;
128  double dtmp,dyiscale=1,ack,scl,rtemp;
129  SDPblk *sdp=&sdpcone->blk[blockj];
130  SDPConeVec W=sdp->W,W2=sdp->W2;
131  DSDPDataMat AA;
132  DSDPVMat T=sdp->T;
133  DSDPDualMat S=sdp->S;
134  DSDPTruth method1;
135  DSDPIndex IS=sdp->IS;
136 
137  DSDPFunctionBegin;
138 
139  info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKERR(info);
140  method1=DSDP_TRUE;
141  method1=DSDP_FALSE;
142 
143  if (method1==DSDP_TRUE){
144  info=DSDPBlockCountNonzeroMatrices(&sdp->ADATA,&nnzmats);DSDPCHKERR(info);
145  for (i=0;i<nnzmats;i++){
146  info=DSDPBlockGetMatrix(&sdp->ADATA,i,&ii,&scl,&AA);DSDPCHKERR(info);
147  info=DSDPVecGetElement(vrow,ii,&dyiscale);DSDPCHKVARERR(ii,info);
148  if (dyiscale==0) continue;
149  info=DSDPDataMatGetRank(AA,&rank,sdp->n);DSDPCHKVARERR(ii,info);
150  for (k=0; k<rank; k++){
151  info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKVARERR(ii,info);
152  if (ack==0) continue;
153  info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKVARERR(ii,info);
154  info=SDPConeVecDot(W,W2,&rtemp); DSDPCHKVARERR(ii,info);
155  dtmp=rtemp*ack*mu*dyiscale*scl;
156  info=DSDPVecAddElement(vrhs2,ii,dtmp);DSDPCHKVARERR(ii,info);
157  }
158  }
159 
160  } else {
161  info=DSDPVMatZeroEntries(T);DSDPCHKERR(info);
162  info=DSDPDualMatInverseAdd(S,mu,T);DSDPCHKERR(info);
163  info=DSDPBlockADot(&sdp->ADATA,1.0,vrow,T,vrhs2);DSDPCHKERR(info);
164  }
165 
166  DSDPFunctionReturn(0);
167 }
168 
169 #undef __FUNCT__
170 #define __FUNCT__ "SDPConeMultiply"
171 
182 int SDPConeMultiply(SDPCone sdpcone, int blockj, double mu, DSDPVec vrow, DSDPVec vin, DSDPVec vout){
183 
184  int info,i,ii,k,rank,nnzmats;
185  double dd2,dtmp,dyiscale,ack,scl,vv;
186  SDPblk *sdp=&sdpcone->blk[blockj];
187  SDPConeVec W=sdp->W,W2=sdp->W2;
188  DSDPDataMat AA;
189  DSDPVMat T=sdp->T;
190  DSDPDSMat DS=sdp->DS;
191  DSDPIndex IS=sdp->IS;
192  DSDPDualMat S=sdp->S;
193 
194  DSDPFunctionBegin;
195  info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKERR(info);
196  info=DSDPVMatZeroEntries(T); DSDPCHKERR(info);
197  info=DSDPBlockASum(&sdp->ADATA,-1.0,vin,T); DSDPCHKERR(info);
198  info=DSDPDSMatSetArray(DS,T); DSDPCHKERR(info);
199  info=DSDPBlockCountNonzeroMatrices(&sdp->ADATA,&nnzmats);DSDPCHKERR(info);
200  for (i=0;i<nnzmats;i++){
201  info=DSDPBlockGetMatrix(&sdp->ADATA,i,&ii,&scl,&AA);DSDPCHKERR(info);
202  info=DSDPVecGetElement(vrow,ii,&dyiscale);DSDPCHKVARERR(ii,info);
203  if (dyiscale==0) continue;
204 
205  info=DSDPDataMatGetRank(AA,&rank,sdp->n);DSDPCHKVARERR(ii,info);
206  for (dd2=0,k=0; k<rank; k++){
207  info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKVARERR(ii,info);
208  if (ack==0) continue;
209  info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKVARERR(ii,info);
210  info=DSDPDSMatVecVec(DS,W2,&vv);DSDPCHKVARERR(ii,info);
211  dd2+=vv*ack;
212  }
213  dtmp = dd2 * dyiscale *mu *scl;
214  info=DSDPVecAddElement(vout,ii,dtmp);DSDPCHKVARERR(ii,info);
215  }
216 
217  DSDPFunctionReturn(0);
218 }
219 
220 
221 
222 #undef __FUNCT__
223 #define __FUNCT__ "SDPConeComputeXX"
224 
235 int SDPConeComputeXX(SDPCone sdpcone, int blockj, DSDPVec DY, double mu,DSDPDualMat S, DSDPVMat X){
236 
237  int info, i, ii,k, rank, n, nnzmats;
238  double dtmp,dyiscale,ack,scl;
239  SDPblk *sdp=&sdpcone->blk[blockj];
240  SDPConeVec W=sdp->W,W2=sdp->W2;
241  DSDPDataMat AA;
242  DSDPIndex IS=sdp->IS;
243 
244  DSDPFunctionBegin;
245  info=SDPConeCheckJ(sdpcone,blockj);DSDPCHKERR(info);
246  mu=mu*sdp->gammamu; n=sdp->n;
247  info=DSDPVMatZeroEntries(X);DSDPCHKERR(info);
248  info=DSDPBlockCountNonzeroMatrices(&sdp->ADATA,&nnzmats);DSDPCHKERR(info);
249  for (i=0;i<nnzmats;i++){
250  info=DSDPBlockGetMatrix(&sdp->ADATA,i,&ii,&scl,&AA);DSDPCHKVARERR(ii,info);
251  info=DSDPVecGetElement(DY,ii,&dyiscale);DSDPCHKVARERR(ii,info);
252  if (dyiscale==0) continue;
253  info=DSDPDataMatGetRank(AA,&rank,n);DSDPCHKVARERR(ii,info);
254  for (k=0; k<rank; k++){
255  info = DSDPDataMatGetEig(AA,k,W,IS,&ack);DSDPCHKVARERR(ii,info);
256  if (ack==0) continue;
257  info=DSDPDualMatInverseMultiply(S,IS,W,W2);DSDPCHKVARERR(ii,info);
258  dtmp = ack * dyiscale * mu * scl;
259  info=DSDPVMatAddOuterProduct(X,dtmp,W2);DSDPCHKVARERR(ii,info);
260  }
261  }
262 
263  info=DSDPDualMatInverseAdd(S,mu,X);DSDPCHKERR(info);
264  DSDPFunctionReturn(0);
265 }
266 
Internal structure for transpose of data.
Definition: dsdpsdp.h:25
DSDPTruth
Boolean variables.
struct DSDPVec_C DSDPVec
This object hold m+2 variables: a scaling of C, the y variables, and r.
Definition: dsdpvec.h:25
int DSDPBlockASum(DSDPBlockData *ADATA, double aa, DSDPVec Yk, DSDPVMat XX)
Sum the data matrices.
Definition: dsdpblock.c:20
int DSDPDualMatInverseMultiply(DSDPDualMat S, DSDPIndex IS, SDPConeVec B, SDPConeVec X)
Multiply the inverse by a vector or solve the system of equations.
Definition: dsdpdualmat.c:236
Schur complement matrix whose solution is the Newton direction.
Definition: dsdpschurmat.h:35
Error handling, printing, and profiling.
int DSDPDSMatSetArray(DSDPDSMat A, DSDPVMat T)
Set values into the matrix.
Definition: dsdpdsmat.c:130
int DSDPDualMatInverseAdd(DSDPDualMat S, double alpha, DSDPVMat T)
Add a multiple of the inverse to T.
Definition: dsdpdualmat.c:209
int SDPConeComputeRHS(SDPCone sdpcone, int blockj, double mu, DSDPVec vrow, DSDPVec vrhs1, DSDPVec vrhs2)
Compute the gradient to the barrier term.
Definition: sdpcompute.c:125
int DSDPDataMatGetRank(DSDPDataMat A, int *rank, int n)
Get the number of nonzero eigenvalues/eigenvectors for the matrix.
Definition: dsdpdatamat.c:129
int DSDPSchurMatRowColumnScaling(DSDPSchurMat, int, DSDPVec, int *)
Get the scaling and nonzero pattern of each column in this row of the matrix.
int DSDPDataMatGetEig(DSDPDataMat A, int rr, SDPConeVec V, DSDPIndex S, double *eigenvalue)
Get an eigenvalue/vector pair.
Definition: dsdpdatamat.c:204
int SDPConeMultiply(SDPCone sdpcone, int blockj, double mu, DSDPVec vrow, DSDPVec vin, DSDPVec vout)
Compute the gradient to the barrier term.
Definition: sdpcompute.c:182
int DSDPBlockvAv(DSDPBlockData *ADATA, double aa, DSDPVec Alpha, SDPConeVec V, DSDPVec VAV)
Set VAV[i] to aa * Alpha[i] * V&#39; A[i] V.
Definition: dsdpblock.c:84
int SDPConeComputeHessian(SDPCone sdpcone, double mu, DSDPSchurMat M, DSDPVec vrhs1, DSDPVec vrhs2)
Compute the Hessian to the barrier term.
Definition: sdpcompute.c:30
int DSDPBlockADot(DSDPBlockData *ADATA, double aa, DSDPVec Alpha, DSDPVMat X, DSDPVec AX)
Compute inner product of XX with data matrices.
Definition: dsdpblock.c:49
int DSDPBlockCountNonzeroMatrices(DSDPBlockData *ADATA, int *nzmats)
Count how many data matrices are in a block of data.
Definition: dsdpblock.c:272
Vector whose length corresponds to dimension of a block in a cone.
Definition: sdpconevec.h:13
Symmetric data matrix for one block in the semidefinite cone.
Definition: dsdpdatamat.h:15
int SDPConeComputeXX(SDPCone sdpcone, int blockj, DSDPVec DY, double mu, DSDPDualMat S, DSDPVMat X)
Compute X.
Definition: sdpcompute.c:235
Internal structure for block of semidefinite cone.
Definition: dsdpsdp.h:52
int SDPConeVecDot(SDPConeVec V1, SDPConeVec V2, double *ans)
Inner product of two vectors.
Definition: sdpconevec.c:125
int DSDPDSMatVecVec(DSDPDSMat A, SDPConeVec X, double *vAv)
Compute the product x&#39; A x.
Definition: dsdpdsmat.c:181
int SDPConeCheckJ(SDPCone sdpcone, int blockj)
Check validity of parameter.
Definition: dsdpadddata.c:31
Represents an S matrix for one block in the semidefinite cone.
Definition: dsdpdualmat.h:18
int DSDPVMatZeroEntries(DSDPVMat X)
Zero matrix.
Definition: dsdpxmat.c:125
int DSDPSchurMatAddRow(DSDPSchurMat, int, double, DSDPVec)
Add elements to a row of the Schur matrix.
int DSDPVMatAddOuterProduct(DSDPVMat X, double alpha, SDPConeVec V)
Add outer product of a vector to the matrix.
Definition: dsdpxmat.c:275
int DSDPBlockGetMatrix(DSDPBlockData *ADATA, int id, int *vari, double *scl, DSDPDataMat *A)
Get a data matrix from a block of data.
Definition: dsdpblock.c:307
Internal structure for semidefinite cone.
Definition: dsdpsdp.h:80
Symmetric Delta S matrix for one block in the semidefinite cone.
Definition: dsdpdsmat.h:23
Dense symmetric matrix for one block in the semidefinite cone.
Definition: dsdpxmat.h:17
Internal SDPCone data structures and routines.