13typedef struct _P_Mat3* Mat3;
17static int ComputeStepFAST(Mat3 A,
SDPConeVec *Q,
int m,
SDPConeVec W,
double *dwork,
int*iwork,
double *maxstep,
double *mineig);
19extern int DSDPGetEigsSTEP(
double[],
int,
double[],
int,
long int[],
int,
20 double[],
int,
double[],
int,
int[],
int);
22int DSDPGetTriDiagonalEigs(
int n,
double *D,
double *E,
double*WORK2N,
int*);
33int DSDPGetTriDiagonalEigs(
int N,
double D[],
double E[],
double WORK[],
int IIWORK[]){
34 ffinteger LDZ=DSDPMax(1,N),INFO,NN=N;
35 ffinteger M,IL=N-1,IU=N,*ISUPPZ=0;
36 ffinteger *IWORK=(ffinteger*)IIWORK,LIWORK,LWORK;
37 double WW[2],VL=-1e10,VU=1e10,*Z=0,ABSTOL=0;
38 char JOBZ=
'N', RANGE=
'I';
40 dstev(&JOBZ,&NN,D,E,Z,&LDZ,WORK,&INFO);
49 dstevr(&JOBZ,&RANGE,&NN,D,E,&VL,&VU,&IL,&IU,&ABSTOL,&M,WW,Z,&LDZ,ISUPPZ,WORK,&LWORK,IWORK,&LIWORK,&INFO);
64#define __FUNCT__ "MatMult3"
68 double minus_one=-1.0;
81 DSDPFunctionReturn(0);
86#define __FUNCT__ "DSDPLanczosInitialize"
108 DSDPFunctionReturn(0);
112#define __FUNCT__ "DSDPSetMaximumLanczosIterations"
121 if (maxlanczos>0) LZ->lanczosm=maxlanczos;
122 DSDPFunctionReturn(0);
126#define __FUNCT__ "DSDPFastLanczosSetup"
137 info=SDPConeVecGetSize(V,&n);DSDPCHKERR(info);
138 LZ->lanczosm=DSDPMin(LZ->maxlanczosm,n);
141 if (LZ->lanczosm<50){
142 DSDPCALLOC2(&LZ->dwork4n,
double,4*(LZ->lanczosm)+2,&info); DSDPCHKERR(info);
143 DSDPCALLOC2(&LZ->iwork10n,
int,1,&info); DSDPCHKERR(info);
145 DSDPCALLOC2(&LZ->dwork4n,
double,23*(LZ->lanczosm)+2,&info); DSDPCHKERR(info);
146 DSDPCALLOC2(&LZ->iwork10n,
int,10*(LZ->lanczosm),&info); DSDPCHKERR(info);
148 DSDPCALLOC2(&LZ->Q,
SDPConeVec,2,&info);DSDPCHKERR(info);
152 DSDPFunctionReturn(0);
156#define __FUNCT__ "DSDPRobustLanczosSetup"
167 info=SDPConeVecGetSize(V,&n);DSDPCHKERR(info);
168 leig=DSDPMin(LZ->maxlanczosm,n);
173 DSDPCALLOC2(&LZ->dwork4n,
double,3*(LZ->lanczosm)+1,&info); DSDPCHKERR(info);
174 DSDPCALLOC2(&LZ->darray,
double,(leig*leig),&info); DSDPCHKERR(info);
175 DSDPCALLOC2(&LZ->Q,
SDPConeVec,(leig+1),&info);DSDPCHKERR(info);
177 for (i=0;i<=leig;i++){
180 info = SDPConeVecCreate(leig,&LZ->Tv);DSDPCHKERR(info);
181 DSDPFunctionReturn(0);
185#define __FUNCT__ "DSDPLanczosDestroy"
195 for (i=0;i<=LZ->lanczosm;i++){
196 info = SDPConeVecDestroy(&LZ->Q[i]);DSDPCHKERR(info);
198 info=SDPConeVecDestroy(&LZ->Tv);DSDPCHKERR(info);
199 DSDPFREE(&LZ->darray,&info);DSDPCHKERR(info);
200 }
else if (LZ->type==1){
201 info = SDPConeVecDestroy(&LZ->Q[1]);DSDPCHKERR(info);
202 info = SDPConeVecDestroy(&LZ->Q[0]);DSDPCHKERR(info);
203 DSDPFREE(&LZ->iwork10n,&info);DSDPCHKERR(info);
205 DSDPFREE(&LZ->Q,&info);DSDPCHKERR(info);
206 DSDPFREE(&LZ->dwork4n,&info);DSDPCHKERR(info);
208 DSDPFunctionReturn(0);
212#define __FUNCT__ "DSDPLanczosMinXEig"
226 info = ComputeStepFAST(A,LZ->Q,m,W1,LZ->dwork4n,LZ->iwork10n,&smaxstep,mineig);DSDPCHKERR(info);
227 }
else if (LZ->type==2){
228 info = ComputeStepROBUST(A,LZ->Q,m,LZ->Q[m],W1,LZ->darray,LZ->Tv,LZ->dwork4n,&smaxstep,mineig);DSDPCHKERR(info);
230 DSDPSETERR1(1,
"Lanczos Step Length Has not been SetUp. Type: %d\n",LZ->type);
232 DSDPFunctionReturn(0);
236#define __FUNCT__ "DSDPLanczosStepSize"
249 double smaxstep,mineig;
260 info = ComputeStepFAST(A,LZ->Q,m,W1,LZ->dwork4n,LZ->iwork10n,&smaxstep,&mineig);DSDPCHKERR(info);
262 }
else if (LZ->type==2){
263 info = ComputeStepROBUST(A,LZ->Q,m,LZ->Q[m],W1,LZ->darray,LZ->Tv,LZ->dwork4n,&smaxstep,&mineig);DSDPCHKERR(info);
266 DSDPSETERR1(1,
"Lanczos Step Length Has not been SetUp. Type: %d\n",LZ->type);
268 DSDPFunctionReturn(0);
274#define __FUNCT__ "ComputeStepROBUST"
278 double tt,wnorm, phi;
279 double lambda1=0,lambda2=0,delta=0;
280 double res1,res2,beta;
287 memset((
void*)darray,0,m*m*
sizeof(
double));
289 for (i=0; i<m; i++){ darray[i*m+i]=-1.0;}
291 for (i=0; i<m; i++){ darray[i*m+i]=1.0;}
298 info = MatMult3(A,Q[i],W);DSDPCHKERR(info);
300 if (phi!=phi){ *maxstep = 0.0;
return 0;}
310 if (wnorm <= 0.8 * phi){
313 if (tt==tt){tt*=-1.0;}
else {tt=0;}
316 if (i!=j){ darray[i*m+j]-=tt; }
322 darray[i*m+i+1]=wnorm;
323 darray[i*m+m+i]=wnorm;
325 if (fabs(wnorm)<=1.0e-14)
break;
339 LWORK=DSDPMax(3*m,1); LDA=DSDPMax(1,m); N=m;
340 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info);
341 info=DSDPGetEigsSTEP(darray,m,0,0,0,0,eigvec,m,dwork,LWORK,0,0);DSDPCHKERR(info);
342 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info);
349 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info);
351 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info);
355 info=SDPConeVecGetArray(QAQTv,&eigvec);DSDPCHKERR(info);
357 lambda1=-eigvec[N-1];
358 lambda2=-eigvec[N-2];
359 info=SDPConeVecRestoreArray(QAQTv,&eigvec);DSDPCHKERR(info);
363 tt=darray[(N-1)*m+i];
366 info = MatMult3(A,W,R);DSDPCHKERR(info);
372 tt=darray[(N-2)*m+i];
375 info = MatMult3(A,W,R);DSDPCHKERR(info);
379 tt = -lambda1 + lambda2 - res2;
382 delta = DSDPMin(res1,sqrt(res1)/beta);
388 *maxstep = 1.0/(delta-lambda1);
392 info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info);
393 DSDPLogInfo(0,19,
"Robust Lanczos StepLength: Iterates %d, Max: %d, BlockSize: %d, Lambda1: %4.2e, Res1: %4.2e, Lambda2: %4.2e, Res2: %4.2e, Delta: %4.2e, MaxStep: %4.2e\n",i,m,n,lambda1,res1*res1,lambda2,res2*res2,delta,*maxstep);
395 DSDPFunctionReturn(0);
399#define __FUNCT__ "ComputeStepFAST"
400static int ComputeStepFAST(Mat3 A,
SDPConeVec *Q,
int m,
SDPConeVec W,
double *dwork,
int *iwork,
double *maxstep ,
double *mineig){
403 double tt,wnorm, phi;
404 double lambda1=0,lambda2=0,delta=0;
405 double res1,res2,beta;
408 double *diag,*subdiag,*ddwork;
416 for (i=0; i<m; i++){ diag[i]=-1; subdiag[i]=0;}
418 for (i=0; i<m; i++){ diag[i]=1.0; subdiag[i]=0;}
424 info = MatMult3(A,Q[0],W);DSDPCHKERR(info);
426 if (phi!=phi){ *maxstep = 0.0;
return 0;}
436 if (wnorm <= 1.0 * phi){
440 if (tt==tt){tt*=-1.0;}
else {tt=0;}
445 if (tt==tt){tt*=-1.0;}
else {tt=0;}
458 if (fabs(wnorm)<=1.0e-10){i++;
break;}
466 info=DSDPGetTriDiagonalEigs(m,diag,subdiag,ddwork,iwork); DSDPCHKERR(info);
483 tt = -lambda1 + lambda2 - res2;
486 delta = DSDPMin(res1,sqrt(res1)/beta);
493 *maxstep = 1.0/(delta-lambda1);
497 info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info);
498 DSDPLogInfo(0,19,
"Step Length: Fast Lanczos Iterates: %2d, Max: %d, Block Size: %d, VNorm: %3.1e, Lambda1: %4.4e, Lambda2: %4.4e, Delta: %4.2e, Maxstep: %4.2e\n",
499 i,m,n,wnorm,lambda1,lambda2,delta,*maxstep);
501 DSDPFunctionReturn(0);
int DSDPDSMatMult(DSDPDSMat A, SDPConeVec X, SDPConeVec Y)
Set values into the matrix.
The interface between the SDPCone and the Delta S matrix.
int DSDPDualMatCholeskySolveBackward(DSDPDualMat S, SDPConeVec B, SDPConeVec X)
Backward triangular solve.
int DSDPDualMatCholeskySolveForward(DSDPDualMat S, SDPConeVec B, SDPConeVec X)
Forward triangular solve.
The interface between the SDPCone and the matrix S.
Lanczos procedure determines the maximum step length.
DSDP uses BLAS and LAPACK for many of its operations.
int DSDPSetMaximumLanczosIterations(DSDPLanczosStepLength *LZ, int maxlanczos)
Set parameter.
int DSDPLanczosInitialize(DSDPLanczosStepLength *LZ)
Initialize Lanczos structure.
int DSDPLanczosDestroy(DSDPLanczosStepLength *LZ)
Free data structure.
int DSDPRobustLanczosSetup(DSDPLanczosStepLength *LZ, SDPConeVec V)
Use slowerer but more robust method.
int DSDPFastLanczosSetup(DSDPLanczosStepLength *LZ, SDPConeVec V)
Use Lanczos procedure. Assume off tridiagonal entries are zero.
int DSDPLanczosStepSize(DSDPLanczosStepLength *LZ, SDPConeVec W1, SDPConeVec W2, DSDPDualMat S, DSDPDSMat DS, double *maxstep)
Compute distance to boundary.
Error handling, printing, and profiling.
int DSDPVMatMult(DSDPVMat X, SDPConeVec Z, SDPConeVec Y)
Multiply X by a vector.
The interface between the SDPCone and the dense matrix array.
int SDPConeVecAXPY(double alpha, SDPConeVec x, SDPConeVec y)
Add a multiple of X to Y.
int SDPConeVecDuplicate(SDPConeVec V1, SDPConeVec *V2)
Allocate another vector with the same structure as the first.
int SDPConeVecNorm2(SDPConeVec VV, double *vnorm)
Compute the Euclidean norm.
int SDPConeVecCopy(SDPConeVec v1, SDPConeVec v2)
Copy v1 to v2.
int SDPConeVecScale(double alpha, SDPConeVec VV)
Compute the Euclidean norm.
int SDPConeVecZero(SDPConeVec V)
Zero the elements of the vector.
int SDPConeVecDot(SDPConeVec V1, SDPConeVec V2, double *ans)
Inner product of two vectors.
int SDPConeVecNormalize(SDPConeVec V)
Scale the vector to norm of 1.
int SDPConeVecSet(double alpha, SDPConeVec V)
Set each element of vector to this number.
Symmetric Delta S matrix for one block in the semidefinite cone.
Represents an S matrix for one block in the semidefinite cone.
Apply Lanczos prodedure to find distance to boundary.
Dense symmetric matrix for one block in the semidefinite cone.
Vector whose length corresponds to dimension of a block in a cone.