00001 #include "dsdp.h"
00002 #include "dsdpsys.h"
00008 #define DSDPCHKCONEERR(a,b); { if (b){ DSDPSETERR1(b,"Cone Number: %d,\n",a);} }
00009
00010 static int ConeSetup=0,ConeComputeS=0,ConeComputeSS=0,ConeComputeH=0,ConeHMultiplyAdd=0,ConeMaxPStep=0,ConeMaxDStep=0,ConePotential=0,ConeComputeX=0,ConeView=0,ConeDestroy=0,ConeXEigs=0,ConeRHS=0,ConeInvertS=0;
00011 static int DSDPRegisterConeEvents(void);
00012 int DSDPSetSchurMatOps(DSDP,struct DSDPSchurMat_Ops*, void*);
00013 int DSDPAddSchurRow(DSDP,int, DSDPVec);
00014
00015
00016
00017
00018 #undef __FUNCT__
00019 #define __FUNCT__ "DSDPZeroConeEvents"
00020 static int DSDPZeroConeEvents(){
00021 DSDPFunctionBegin;
00022 ConeSetup=0;ConeComputeS=0;ConeComputeSS=0;ConeComputeH=0;ConeHMultiplyAdd=0;ConeMaxPStep=0;ConeMaxDStep=0;ConePotential=0;ConeComputeX=0;ConeView=0;ConeDestroy=0;ConeXEigs=0;ConeRHS=0;ConeInvertS=0;
00023 DSDPFunctionReturn(0);
00024 }
00025
00026 #undef __FUNCT__
00027 #define __FUNCT__ "DSDPRegisterConeEvents"
00028 static int DSDPRegisterConeEvents(){
00029
00030 DSDPFunctionBegin;
00031 DSDPEventLogRegister("Cone Setup 1&2",&ConeSetup);
00032 DSDPEventLogRegister("Cone Invert S",&ConeInvertS);
00033 DSDPEventLogRegister("Cone RHS",&ConeRHS);
00034 DSDPEventLogRegister("Cone Compute Newton Eq.",&ConeComputeH);
00035 DSDPEventLogRegister("Cone Newton Multiply-Add",&ConeHMultiplyAdd);
00036 DSDPEventLogRegister("Cone Max P Step Length",&ConeMaxPStep);
00037 DSDPEventLogRegister("Cone Compute and Factor SP",&ConeComputeSS);
00038 DSDPEventLogRegister("Cone Max D Step Length",&ConeMaxDStep);
00039 DSDPEventLogRegister("Cone Compute and Factor S",&ConeComputeS);
00040 DSDPEventLogRegister("Cone Potential",&ConePotential);
00041 DSDPEventLogRegister("Cone View",&ConeView);
00042 DSDPEventLogRegister("Cone Compute X",&ConeComputeX);
00043 DSDPEventLogRegister("Cone X Residuals",&ConeXEigs);
00044 DSDPEventLogRegister("Cone Destroy",&ConeDestroy);
00045
00046 DSDPFunctionReturn(0);
00047 }
00048
00056 #undef __FUNCT__
00057 #define __FUNCT__ "DSDPSetUpCones"
00058 int DSDPSetUpCones(DSDP dsdp){
00059 int info,kk;
00060 DSDPVec yy0=dsdp->y;
00061 DSDPFunctionBegin;
00062 info=DSDPRegisterConeEvents();
00063 DSDPEventLogBegin(ConeSetup);
00064 for (kk=0;kk<dsdp->ncones;kk++){
00065 DSDPEventLogBegin(dsdp->K[kk].coneid);
00066 info=DSDPConeSetUp(dsdp->K[kk].cone,yy0);DSDPCHKCONEERR(kk,info);
00067 DSDPEventLogEnd(dsdp->K[kk].coneid);
00068 }
00069 DSDPEventLogEnd(ConeSetup);
00070 DSDPFunctionReturn(0);
00071 }
00072
00082 #undef __FUNCT__
00083 #define __FUNCT__ "DSDPSetUpCones2"
00084 int DSDPSetUpCones2(DSDP dsdp, DSDPVec yy0, DSDPSchurMat M){
00085 int info,kk;
00086 DSDPFunctionBegin;
00087 DSDPEventLogBegin(ConeSetup);
00088 for (kk=0;kk<dsdp->ncones;kk++){
00089 DSDPEventLogBegin(dsdp->K[kk].coneid);
00090 info=DSDPConeSetUp2(dsdp->K[kk].cone,yy0,M);DSDPCHKCONEERR(kk,info);
00091 DSDPEventLogEnd(dsdp->K[kk].coneid);
00092 }
00093 DSDPEventLogEnd(ConeSetup);
00094 DSDPFunctionReturn(0);
00095 }
00096
00097
00105 #undef __FUNCT__
00106 #define __FUNCT__ "DSDPDestroyCones"
00107 int DSDPDestroyCones(DSDP dsdp){
00108 int info,kk,ncones=dsdp->ncones;
00109 DSDPFunctionBegin;
00110 DSDPEventLogBegin(ConeDestroy);
00111 for (kk=ncones-1;kk>=0; kk--){
00112 DSDPEventLogBegin(dsdp->K[kk].coneid);
00113 info=DSDPConeDestroy(&dsdp->K[kk].cone);DSDPCHKCONEERR(kk,info);
00114 DSDPEventLogEnd(dsdp->K[kk].coneid);
00115 info=DSDPConeInitialize(&dsdp->K[kk].cone);DSDPCHKCONEERR(kk,info);
00116 dsdp->ncones--;
00117 }
00118 if (dsdp->maxcones>0){
00119 DSDPFREE(&dsdp->K,&info);DSDPCHKERR(info);
00120 dsdp->K=0;
00121 dsdp->maxcones=0;
00122 }
00123 DSDPEventLogEnd(ConeDestroy);
00124 info=DSDPZeroConeEvents();DSDPCHKERR(info);
00125 DSDPFunctionReturn(0);
00126 }
00127
00128
00140 #undef __FUNCT__
00141 #define __FUNCT__ "DSDPComputeHessian"
00142 int DSDPComputeHessian( DSDP dsdp , DSDPSchurMat M, DSDPVec vrhs1, DSDPVec vrhs2){
00143 int info,kk; double r;
00144 DSDPFunctionBegin;
00145 DSDPEventLogBegin(ConeComputeH);
00146 dsdp->schurmu=dsdp->mutarget;
00147 info=DSDPVecGetR(dsdp->y,&r);DSDPCHKERR(info);
00148 info=DSDPSchurMatSetR(dsdp->M,r);DSDPCHKERR(info);
00149 info=DSDPSchurMatZeroEntries(M);DSDPCHKERR(info);
00150 info=DSDPVecZero(vrhs1);DSDPCHKERR(info);
00151 info=DSDPVecZero(vrhs2);DSDPCHKERR(info);
00152 info=DSDPVecZero(M.schur->rhs3);DSDPCHKERR(info);
00153 info=DSDPObjectiveGH(dsdp,M,vrhs1); DSDPCHKERR(info);
00154 for (kk=dsdp->ncones-1;kk>=0;kk--){
00155 DSDPEventLogBegin(dsdp->K[kk].coneid);
00156 info=DSDPConeComputeHessian(dsdp->K[kk].cone,dsdp->schurmu,M,vrhs1,vrhs2);DSDPCHKCONEERR(kk,info);
00157 DSDPEventLogEnd(dsdp->K[kk].coneid);
00158 }
00159 info=DSDPSchurMatAssemble(M);DSDPCHKERR(info);
00160
00161 info=DSDPSchurMatReducePVec(M,vrhs1);DSDPCHKERR(info);
00162 info=DSDPSchurMatReducePVec(M,vrhs2);DSDPCHKERR(info);
00163 info=DSDPSchurMatReducePVec(M,M.schur->rhs3);DSDPCHKERR(info);
00164 if (0 && dsdp->UsePenalty==DSDPNever){
00165 info=DSDPVecAXPY(1.0,M.schur->rhs3,vrhs2);DSDPCHKERR(info);
00166 info=DSDPVecZero(M.schur->rhs3);DSDPCHKERR(info);
00167 info=DSDPVecZero(M.schur->dy3);DSDPCHKERR(info);
00168 info=DSDPVecSetR(vrhs1,0);DSDPCHKERR(info);
00169 info=DSDPVecSetR(vrhs2,r);DSDPCHKERR(info);
00170 }
00171 DSDPEventLogEnd(ConeComputeH);
00172 DSDPFunctionReturn(0);
00173 }
00174
00175
00176 #undef __FUNCT__
00177 #define __FUNCT__ "DSDPHessianMultiplyAdd"
00178
00188 int DSDPHessianMultiplyAdd( DSDP dsdp , DSDPVec v, DSDPVec vv){
00189 int info,kk;
00190 DSDPVec vrow=dsdp->sles->BR;
00191 DSDPFunctionBegin;
00192 DSDPEventLogBegin(ConeHMultiplyAdd);
00193
00194 info=DSDPSchurMatRowScaling(dsdp->M,vrow);DSDPCHKERR(info);
00195 for (kk=0;kk<dsdp->ncones;kk++){
00196 DSDPEventLogBegin(dsdp->K[kk].coneid);
00197 info=DSDPConeMultiplyAdd(dsdp->K[kk].cone,dsdp->schurmu,vrow,v,vv);DSDPCHKCONEERR(kk,info);
00198 DSDPEventLogEnd(dsdp->K[kk].coneid);
00199 }
00200 info=DSDPSchurMatReducePVec(dsdp->M,vv);DSDPCHKERR(info);
00201 DSDPEventLogEnd(ConeHMultiplyAdd);
00202 DSDPFunctionReturn(0);
00203 }
00204
00205 #undef __FUNCT__
00206 #define __FUNCT__ "DSDPComputeG"
00207
00215 int DSDPComputeG( DSDP dsdp , DSDPVec vt, DSDPVec vrhs1, DSDPVec vrhs2){
00216 int info,kk; double r;
00217 DSDPFunctionBegin;
00218 DSDPEventLogBegin(ConeRHS);
00219 info=DSDPVecZero(vrhs1);DSDPCHKERR(info);
00220 info=DSDPVecZero(vrhs2);DSDPCHKERR(info);
00221 info=DSDPVecGetR(dsdp->y,&r);DSDPCHKERR(info);
00222 info=DSDPSchurMatSetR(dsdp->M,r);DSDPCHKERR(info);
00223 info=DSDPSchurMatRowScaling(dsdp->M,vt);DSDPCHKERR(info);
00224 info=DSDPObjectiveGH(dsdp,dsdp->M,vrhs1); DSDPCHKERR(info);
00225 if (0 && r==0){info=DSDPVecSetR(vrhs1,0);info=DSDPVecSetR(vrhs2,0);}
00226
00227 for (kk=0;kk<dsdp->ncones;kk++){
00228 DSDPEventLogBegin(dsdp->K[kk].coneid);
00229 info=DSDPConeComputeRHS(dsdp->K[kk].cone,dsdp->schurmu,vt,vrhs1,vrhs2);DSDPCHKCONEERR(kk,info);
00230 DSDPEventLogEnd(dsdp->K[kk].coneid);
00231 }
00232 DSDPEventLogEnd(ConeRHS);
00233 info=DSDPSchurMatReducePVec(dsdp->M,vrhs1);DSDPCHKERR(info);
00234 info=DSDPSchurMatReducePVec(dsdp->M,vrhs2);DSDPCHKERR(info);
00235 DSDPFunctionReturn(0);
00236 }
00237
00238 #undef __FUNCT__
00239 #define __FUNCT__ "DSDPComputeANorm2"
00240
00246 int DSDPComputeANorm2( DSDP dsdp , DSDPVec Anorm2){
00247 int info,kk;
00248 DSDPFunctionBegin;
00249 for (kk=0;kk<dsdp->ncones;kk++){
00250 DSDPEventLogBegin(dsdp->K[kk].coneid);
00251 info=DSDPConeANorm2(dsdp->K[kk].cone,Anorm2);DSDPCHKCONEERR(kk,info);
00252 DSDPEventLogEnd(dsdp->K[kk].coneid);
00253 }
00254 DSDPFunctionReturn(0);
00255 }
00256
00257
00270 #undef __FUNCT__
00271 #define __FUNCT__ "DSDPComputeSS"
00272 int DSDPComputeSS(DSDP dsdp, DSDPVec Y, DSDPDualFactorMatrix flag, DSDPTruth *ispsdefinite){
00273 int info,kk;
00274 DSDPTruth psd=DSDP_TRUE;
00275 DSDPFunctionBegin;
00276 if (flag==DUAL_FACTOR){
00277 DSDPEventLogBegin(ConeComputeS);
00278 } else if (flag==PRIMAL_FACTOR){
00279 DSDPEventLogBegin(ConeComputeSS);
00280 }
00281 for (kk=dsdp->ncones-1; kk>=0 && psd==DSDP_TRUE;kk--){
00282 DSDPEventLogBegin(dsdp->K[kk].coneid);
00283 info=DSDPConeComputeS(dsdp->K[kk].cone,Y,flag,&psd); DSDPCHKCONEERR(kk,info);
00284 DSDPEventLogEnd(dsdp->K[kk].coneid);
00285 }
00286 *ispsdefinite=psd;
00287 if (flag==DUAL_FACTOR){
00288 DSDPEventLogEnd(ConeComputeS);
00289 } else if (flag==PRIMAL_FACTOR){
00290 DSDPEventLogEnd(ConeComputeSS);
00291 }
00292 DSDPFunctionReturn(0);
00293 }
00294
00305 #undef __FUNCT__
00306 #define __FUNCT__ "DSDPInvertS"
00307 int DSDPInvertS(DSDP dsdp){
00308 int info,kk;
00309 DSDPFunctionBegin;
00310 DSDPEventLogBegin(ConeInvertS);
00311 for (kk=0;kk<dsdp->ncones;kk++){
00312 DSDPEventLogBegin(dsdp->K[kk].coneid);
00313 info=DSDPConeInvertS(dsdp->K[kk].cone); DSDPCHKCONEERR(kk,info);
00314 DSDPEventLogEnd(dsdp->K[kk].coneid);
00315 }
00316 DSDPEventLogEnd(ConeInvertS);
00317 DSDPFunctionReturn(0);
00318 }
00319
00334 #undef __FUNCT__
00335 #define __FUNCT__ "DSDPComputeMaxStepLength"
00336 int DSDPComputeMaxStepLength(DSDP dsdp, DSDPVec DY, DSDPDualFactorMatrix flag, double *maxsteplength){
00337 int info,kk;
00338 double msteplength=1.0e30,conesteplength;
00339 DSDPFunctionBegin;
00340 if (flag==DUAL_FACTOR){
00341 DSDPEventLogBegin(ConeMaxDStep);
00342 } else if (flag==PRIMAL_FACTOR){
00343 DSDPEventLogBegin(ConeMaxPStep);
00344 }
00345 for (kk=0;kk<dsdp->ncones;kk++){
00346 DSDPEventLogBegin(dsdp->K[kk].coneid);
00347 conesteplength=1.0e20;
00348 info=DSDPConeComputeMaxStepLength(dsdp->K[kk].cone,DY,flag,&conesteplength);DSDPCHKCONEERR(kk,info);
00349 msteplength=DSDPMin(msteplength,conesteplength);
00350 DSDPEventLogEnd(dsdp->K[kk].coneid);
00351 }
00352 *maxsteplength=msteplength;
00353 if (flag==DUAL_FACTOR){
00354 DSDPEventLogEnd(ConeMaxDStep);
00355 } else if (flag==PRIMAL_FACTOR){
00356 DSDPEventLogEnd(ConeMaxPStep);
00357 }
00358 DSDPFunctionReturn(0);
00359 }
00360
00361
00376 #undef __FUNCT__
00377 #define __FUNCT__ "DSDPPassXVectors"
00378 int DSDPPassXVectors(DSDP dsdp, double mu, DSDPVec Y, DSDPVec DY){
00379 int info,kk;
00380 DSDPFunctionBegin;
00381 for (kk=0;kk<dsdp->ncones;kk++){
00382 DSDPEventLogBegin(dsdp->K[kk].coneid);
00383 info=DSDPConeSetXMaker(dsdp->K[kk].cone,mu,Y,DY);DSDPCHKCONEERR(kk,info);
00384 DSDPEventLogEnd(dsdp->K[kk].coneid);
00385 }
00386 DSDPFunctionReturn(0);
00387 }
00388
00389
00399 #undef __FUNCT__
00400 #define __FUNCT__ "DSDPGetConicDimension"
00401 int DSDPGetConicDimension(DSDP dsdp, double *n){
00402 int info,kk;
00403 double nn,nnn=0;
00404 DSDPFunctionBegin;
00405 for (kk=0;kk<dsdp->ncones;kk++){
00406 nn=0;
00407 info=DSDPConeGetDimension(dsdp->K[kk].cone,&nn);DSDPCHKCONEERR(kk,info);
00408 nnn+=nn;
00409 }
00410 *n=nnn;
00411 DSDPFunctionReturn(0);
00412 }
00413
00414
00422 #undef __FUNCT__
00423 #define __FUNCT__ "DSDPViewCones"
00424 int DSDPViewCones(DSDP dsdp){
00425 int info,kk;
00426 DSDPFunctionBegin;
00427 DSDPEventLogBegin(ConeView);
00428 for (kk=0;kk<dsdp->ncones;kk++){
00429 DSDPEventLogBegin(dsdp->K[kk].coneid);
00430 info=DSDPConeView(dsdp->K[kk].cone);DSDPCHKCONEERR(kk,info);
00431 DSDPEventLogEnd(dsdp->K[kk].coneid);
00432 }
00433 DSDPEventLogEnd(ConeView);
00434 DSDPFunctionReturn(0);
00435 }
00436
00437
00448 #undef __FUNCT__
00449 #define __FUNCT__ "DSDPMonitorCones"
00450 int DSDPMonitorCones(DSDP dsdp,int tag){
00451 int info,kk;
00452 DSDPFunctionBegin;
00453 DSDPEventLogBegin(ConeView);
00454 for (kk=0;kk<dsdp->ncones;kk++){
00455 DSDPEventLogBegin(dsdp->K[kk].coneid);
00456 info=DSDPConeMonitor(dsdp->K[kk].cone,tag);DSDPCHKCONEERR(kk,info);
00457 DSDPEventLogEnd(dsdp->K[kk].coneid);
00458 }
00459 DSDPEventLogEnd(ConeView);
00460 DSDPFunctionReturn(0);
00461 }
00462
00463
00472 #undef __FUNCT__
00473 #define __FUNCT__ "DSDPSparsityInSchurMat"
00474 int DSDPSchurSparsity(DSDP dsdp, int row, int rnnz[], int m){
00475 int info,kk;
00476 DSDPFunctionBegin;
00477 for (kk=0;kk<dsdp->ncones;kk++){
00478
00479 info=DSDPConeSparsityInSchurMat(dsdp->K[kk].cone,row,rnnz,m+2);DSDPCHKCONEERR(kk,info);
00480
00481 }
00482 DSDPFunctionReturn(0);
00483 }
00484
00493 #undef __FUNCT__
00494 #define __FUNCT__ "DSDPComputeLogSDeterminant"
00495 int DSDPComputeLogSDeterminant(DSDP dsdp, double *logdet){
00496 int info,kk;
00497 double coneobjective,conepotential,llogdet=0;
00498 DSDPFunctionBegin;
00499 DSDPEventLogBegin(ConePotential);
00500 for (kk=0;kk<dsdp->ncones;kk++){
00501 DSDPEventLogBegin(dsdp->K[kk].coneid);
00502 coneobjective=0;conepotential=0;
00503 info=DSDPConeComputeLogSDeterminant(dsdp->K[kk].cone,&coneobjective,&conepotential);DSDPCHKCONEERR(kk,info);
00504 llogdet+=conepotential;
00505 DSDPEventLogEnd(dsdp->K[kk].coneid);
00506 }
00507 *logdet=llogdet;
00508 DSDPEventLogEnd(ConePotential);
00509 DSDPFunctionReturn(0);
00510 }
00511
00512
00520 #undef __FUNCT__
00521 #define __FUNCT__ "DSDPSetCone"
00522 int DSDPSetCone(DSDP dsdp,DSDPCone tcone){
00523 int info,i,tc;
00524 char conename[100];
00525 DCone *ccones;
00526 DSDPFunctionBegin;
00527 if (dsdp->ncones>=dsdp->maxcones){
00528 tc=2*dsdp->maxcones+4;
00529
00530 DSDPCALLOC2(&ccones,DCone,tc,&info);DSDPCHKERR(info);
00531 for (i=0;i<dsdp->ncones;i++){ccones[i].cone=dsdp->K[i].cone;}
00532 for (i=0;i<dsdp->ncones;i++){ccones[i].coneid=dsdp->K[i].coneid;}
00533 DSDPFREE(&dsdp->K,&info);DSDPCHKERR(info);
00534 dsdp->K=ccones;
00535 dsdp->maxcones=tc;
00536 }
00537 info=DSDPGetConeName(tcone,conename,100);DSDPCHKERR(info);
00538 DSDPEventLogRegister(conename,&tc);
00539 dsdp->K[dsdp->ncones].cone=tcone;
00540 dsdp->K[dsdp->ncones].coneid=tc;
00541 dsdp->ncones++;
00542 DSDPFunctionReturn(0);
00543 }
00544
00567 #undef __FUNCT__
00568 #define __FUNCT__ "DSDPAddCone"
00569 int DSDPAddCone(DSDP dsdp,struct DSDPCone_Ops* dsdpops, void* dsdpcone){
00570 int info;
00571 DSDPCone K;
00572 DSDPFunctionBegin;
00573 info=DSDPConeInitialize(&K); DSDPCHKERR(info);
00574 info=DSDPConeSetData(&K,dsdpops,dsdpcone); DSDPCHKERR(info);
00575 info=DSDPSetCone(dsdp,K); DSDPCHKERR(info);
00576 DSDPFunctionReturn(0);
00577 }
00578
00579
00600 #undef __FUNCT__
00601 #define __FUNCT__ "DSDPSetSchurMatOps"
00602 int DSDPSetSchurMatOps(DSDP dsdp,struct DSDPSchurMat_Ops* sops, void* mdata){
00603 int info;
00604 DSDPFunctionBegin;
00605 info=DSDPSchurMatSetData(&dsdp->M,sops,mdata);DSDPCHKERR(info);
00606 DSDPFunctionReturn(0);
00607 }
00608
00609
00620 #undef __FUNCT__
00621 #define __FUNCT__ "DSDPSetSchurRow"
00622 int DSDPAddSchurRow(DSDP dsdp,int row, DSDPVec R){
00623 int info;
00624 DSDPFunctionBegin;
00625 info=DSDPSchurMatAddRow(dsdp->M,row,1.0,R);DSDPCHKERR(info);
00626 DSDPFunctionReturn(0);
00627 }
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00652 #undef __FUNCT__
00653 #define __FUNCT__ "DSDPComputeXVariables"
00654 int DSDPComputeXVariables(DSDP dsdp, double xmakermu, DSDPVec xmakery, DSDPVec xmakerdy, DSDPVec AX, double *tracexs){
00655 int kk,info;
00656 double ttracexs=0,tttracexs=0,tracex;
00657
00658 DSDPFunctionBegin;
00659 DSDPEventLogBegin(ConeComputeX);
00660 info=DSDPVecZero(AX);DSDPCHKERR(info);
00661 for (kk=0;kk<dsdp->ncones;kk++){
00662 DSDPEventLogBegin(dsdp->K[kk].coneid);
00663 tttracexs=0;
00664 info=DSDPConeComputeX(dsdp->K[kk].cone,xmakermu,xmakery,xmakerdy,AX,&tttracexs);DSDPCHKCONEERR(kk,info);
00665 ttracexs+=tttracexs;
00666 DSDPEventLogEnd(dsdp->K[kk].coneid);
00667 }
00668 info=DSDPVecGetR(AX,&tracex); DSDPCHKERR(info);
00669 DSDPLogInfo(0,2,"Trace(X): %4.2e\n",dsdp->tracex);
00670 info=DSDPVecAXPY(-1.0,dsdp->b,AX); DSDPCHKERR(info);
00671 info=DSDPComputeFixedYX(dsdp->M,AX); DSDPCHKERR(info);
00672 *tracexs=ttracexs;
00673 info=DSDPVecSetR(AX,tracex); DSDPCHKERR(info);
00674 DSDPEventLogEnd(ConeComputeX);
00675 DSDPFunctionReturn(0);
00676 }
00677
00678