11 static int hfactorevent=0,hsolveevent=0;
13 #define DSDPNoOperationError(a); { DSDPSETERR1(10,"Schur matrix type: %s, Operation not defined\n",(a).dsdpops->matname); }
14 #define DSDPChkMatError(a,b); { if (b){ DSDPSETERR1(b,"Schur matrix type: %s,\n",(a).dsdpops->matname);} }
20 #define __FUNCT__ "DSDPSchurMatSetData"
32 DSDPFunctionReturn(0);
35 static const char* schurmatname=
"NOT NAMED YET";
38 #define __FUNCT__ "DSDPSchurMatOpsInitialize"
46 if (dops==NULL)
return 0;
48 dops->matrownonzeros=0;
50 dops->mataddelement=0;
51 dops->matadddiagonal=0;
52 dops->matshiftdiagonal=0;
54 dops->matscaledmultiply=0;
58 dops->pmatonprocessor=0;
59 dops->pmatwhichdiag=0;
60 dops->pmatdistributed=0;
65 dops->matname=schurmatname;
66 DSDPFunctionReturn(0);
69 static struct DSDPSchurMat_Ops dsdpmops;
73 #define __FUNCT__ "DSDPSchurMatOpsInitialize"
84 DSDPCALLOC1(&M->schur,DSDPSchurInfo,&info);DSDPCHKERR(info);
85 M->schur->m=0; M->schur->r=0; M->schur->dd=0;
86 info=DSDPInitializeFixedVariable(&M->schur->fv);DSDPCHKERR(info);
87 DSDPFunctionReturn(0);
91 #define __FUNCT__ "DSDPSchurMatZeroEntries"
100 if (M.dsdpops->matzero){
101 info=(M.dsdpops->matzero)(M.data); DSDPChkMatError(M,info);
103 DSDPNoOperationError(M);
105 DSDPFunctionReturn(0);
110 #define __FUNCT__ "DSDPSchurMatShiftDiagonal"
123 if (dd==0){DSDPFunctionReturn(0);}
125 if (M.dsdpops->matshiftdiagonal){
127 info=(M.dsdpops->matshiftdiagonal)(M.data,dd); DSDPChkMatError(M,info);
128 DSDPLogInfo(0,2,
"Add %4.4e to the Diagonal of Schur Matrix\n",dd);
130 DSDPNoOperationError(M);
132 DSDPFunctionReturn(0);
137 #define __FUNCT__ "DSDPSchurMatInParallel"
152 if (M.dsdpops->pmatdistributed){
153 info=(M.dsdpops->pmatdistributed)(M.data,&flg); DSDPChkMatError(M,info);
158 DSDPFunctionReturn(0);
164 #define __FUNCT__ "DSDPSchurMatAssemble"
177 if (M.dsdpops->matassemble){
178 info=(M.dsdpops->matassemble)(M.data); DSDPChkMatError(M,info);
180 DSDPNoOperationError(M);
182 DSDPFunctionReturn(0);
186 #define __FUNCT__ "DSDPSchurMatFactor"
198 DSDPVec rhs3=M.schur->rhs3,dy3=M.schur->dy3;
201 DSDPEventLogBegin(hfactorevent);
202 if (M.dsdpops->matfactor){
203 info=(M.dsdpops->matfactor)(M.data,&flag); DSDPChkMatError(M,info);
206 DSDPLogInfo(0,2,
"Indefinite Schur Matrix -- Bad Factorization\n");
209 DSDPNoOperationError(M);
211 DSDPEventLogEnd(hfactorevent);
213 info=DSDPSchurMatSolveM(M,rhs3,dy3);DSDPCHKERR(info);}
214 else {info=DSDPVecZero(dy3);DSDPCHKERR(info);}
215 DSDPFunctionReturn(0);
220 #define __FUNCT__ "DSDPSchurMatMultiply"
233 double *xx,*yy,r=M.schur->r;
238 if (M.dsdpops->matscaledmultiply){
239 info=DSDPVecGetSize(x,&n); DSDPCHKERR(info);
240 info=DSDPVecGetArray(x,&xx); DSDPCHKERR(info);
241 info=DSDPVecGetArray(y,&yy); DSDPCHKERR(info);
242 info=(M.dsdpops->matscaledmultiply)(M.data,xx+1,yy+1,n-2); DSDPChkMatError(M,info);
245 info=DSDPVecRestoreArray(y,&yy); DSDPCHKERR(info);
246 info=DSDPVecRestoreArray(x,&xx); DSDPCHKERR(info);
248 DSDPNoOperationError(M);
252 info=DSDPVecGetR(rhs3,&r2);DSDPCHKERR(info);
253 info=DSDPVecGetR(x,&r1);DSDPCHKERR(info);
254 info=DSDPVecAXPY(r1,rhs3,y);DSDPCHKERR(info);
255 info=DSDPVecDot(rhs3,x,&dd);DSDPCHKERR(info);
256 info=DSDPVecAddR(y,dd-r1*r2);DSDPCHKERR(info);
258 DSDPFunctionReturn(0);
262 #define __FUNCT__ "DSDPSchurMatMultR"
265 double *xx,*yy,r=M.schur->r;
270 if (M.dsdpops->matmultr){
271 info=DSDPVecGetSize(x,&n); DSDPCHKERR(info);
272 info=DSDPVecGetArray(x,&xx); DSDPCHKERR(info);
273 info=DSDPVecGetArray(y,&yy); DSDPCHKERR(info);
274 info=(M.dsdpops->matmultr)(M.data,xx+1,yy+1,n-2); DSDPChkMatError(M,info);
277 info=DSDPVecRestoreArray(y,&yy); DSDPCHKERR(info);
278 info=DSDPVecRestoreArray(x,&xx); DSDPCHKERR(info);
281 info=DSDPVecGetR(rhs3,&r2);DSDPCHKERR(info);
282 info=DSDPVecGetR(x,&r1);DSDPCHKERR(info);
283 info=DSDPVecAXPY(r1,rhs3,y);DSDPCHKERR(info);
284 info=DSDPVecDot(rhs3,x,&dd);DSDPCHKERR(info);
285 info=DSDPVecAddR(y,dd-r1*r2);DSDPCHKERR(info);
288 info=DSDPVecZero(y);DSDPCHKERR(info);
291 DSDPFunctionReturn(0);
295 #define __FUNCT__ "DSDPSchurMatReducePVec"
313 if (M.dsdpops->pmatreduction){
314 info=DSDPVecGetSize(x,&n); DSDPCHKERR(info);
315 info=DSDPVecGetArray(x,&xx); DSDPCHKERR(info);
316 info=(M.dsdpops->pmatreduction)(M.data,xx+1,n-2); DSDPChkMatError(M,info);
317 info=DSDPVecRestoreArray(x,&xx); DSDPCHKERR(info);
321 DSDPNoOperationError(M);
324 info=DSDPZeroFixedVariables(M,x);DSDPCHKERR(info);
325 DSDPFunctionReturn(0);
331 #define __FUNCT__ "DSDPSchurMatSetR"
341 DSDPFunctionReturn(0);
345 #define __FUNCT__ "DSDPSchurMatSetup"
355 info=DSDPVecDuplicate(Y,&M.schur->rhs3);
356 info=DSDPVecDuplicate(Y,&M.schur->dy3);
357 info=DSDPVecGetSize(Y,&m);DSDPCHKERR(info);
358 if (M.dsdpops->matsetup){
359 info=(M.dsdpops->matsetup)(M.data,m-2); DSDPChkMatError(M,info);
361 DSDPNoOperationError(M);
363 DSDPEventLogRegister(
"Factor Newton Eq.",&hfactorevent);
364 DSDPEventLogRegister(
"Solve Newton Eq.",&hsolveevent);
365 DSDPFunctionReturn(0);
370 #define __FUNCT__ "DSDPSchurMatView"
379 if (M.dsdpops->matview){
380 info=(M.dsdpops->matview)(M.data); DSDPChkMatError(M,info);
382 DSDPNoOperationError(M);
384 info=DSDPVecView(M.schur->rhs3);DSDPCHKERR(info);
385 DSDPFunctionReturn(0);
390 #define __FUNCT__ "DSDPSchurMatRowScaling"
403 info=DSDPZeroFixedVariables(M,D);DSDPCHKERR(info);
404 DSDPFunctionReturn(0);
408 #define __FUNCT__ "DSDPSchurMatDestroy"
417 if ((*M).dsdpops->matdestroy){
418 info=((*M).dsdpops->matdestroy)((*M).data); DSDPChkMatError(*M,info);
424 info=DSDPVecDestroy(&M->schur->rhs3);DSDPCHKERR(info);
425 info=DSDPVecDestroy(&M->schur->dy3);DSDPCHKERR(info);
429 DSDPFREE(&M->schur,&info);DSDPCHKERR(info);
430 DSDPFunctionReturn(0);
434 #define __FUNCT__ "DSDPSchurMatSolveM"
439 info=DSDPEventLogBegin(hsolveevent);
440 if (M.dsdpops->matsolve){
441 info=DSDPVecGetArray(b,&bb); DSDPCHKERR(info);
442 info=DSDPVecGetSize(x,&n); DSDPCHKERR(info);
443 info=DSDPVecZero(x);DSDPCHKERR(info);
444 info=DSDPVecGetArray(x,&xx); DSDPCHKERR(info);
445 info=(M.dsdpops->matsolve)(M.data,bb+1,xx+1,n-2); DSDPChkMatError(M,info);
446 info=DSDPVecRestoreArray(b,&bb); DSDPCHKERR(info);
447 info=DSDPVecRestoreArray(x,&xx); DSDPCHKERR(info);
449 DSDPNoOperationError(M);
451 info=DSDPVecSetR(x,0.0);DSDPCHKERR(info);
452 info=DSDPVecSetC(x,0.0);DSDPCHKERR(info);
453 info=DSDPEventLogEnd(hsolveevent);
454 DSDPFunctionReturn(0);
458 #define __FUNCT__ "DSDPSchurMatSolve"
469 info=DSDPSchurMatSolveM(M,b,x);DSDPCHKERR(info);
470 info=DSDPApplySMW(M,b,x);DSDPCHKERR(info);
471 info=DSDPZeroFixedVariables(M,x);DSDPCHKERR(info);
472 DSDPFunctionReturn(0);
476 #define __FUNCT__ "DSDPApplySMW"
479 double r=M.schur->r,rr,dr,rhsr,rssr;
480 double rhsnorm,rhsnorm3,rhs1mrhs3=0,rhs3mrhs3=0;
481 DSDPVec rhs3=M.schur->rhs3,dy3=M.schur->dy3;
484 info=DSDPVecNormInfinity(rhs,&rhsnorm);DSDPCHKERR(info);
485 info=DSDPVecNormInfinity(rhs3,&rhsnorm3);DSDPCHKERR(info);
486 if (r==0 || rhsnorm==0){
487 info=DSDPVecSetR(dy,0); DSDPCHKERR(info);
488 info=DSDPVecSetR(rhs,0); DSDPCHKERR(info);
489 }
else if (0 && rhsnorm3==0){
490 info=DSDPVecGetR(rhs,&rr); DSDPCHKERR(info);
491 info=DSDPVecSetR(dy,rr); DSDPCHKERR(info);
494 info=DSDPVecGetR(rhs,&rhsr); DSDPCHKERR(info);
495 info=DSDPVecGetR(rhs3,&rssr); DSDPCHKERR(info);
496 info=DSDPVecDot(rhs3,dy,&rhs1mrhs3); DSDPCHKERR(info);
497 info=DSDPVecDot(rhs3,dy3,&rhs3mrhs3); DSDPCHKERR(info);
498 if (rssr-rhs3mrhs3==0) rssr*=(1.00001);
499 dr=-(rhs1mrhs3-rhsr )/(rssr-rhs3mrhs3);
500 info=DSDPVecAXPY(-dr,dy3,dy);DSDPCHKERR(info);
501 info=DSDPVecSetR(dy,dr); DSDPCHKERR(info);
502 info=DSDPVecSetR(rhs,rhsr); DSDPCHKERR(info);
503 info=DSDPVecDot(rhs,dy,&rhs3mrhs3); DSDPCHKERR(info);
505 DSDPLogInfo(0,3,
"DSDP Step Direction Not Descent, Adjusting. \n");
506 info=DSDPVecAddR(rhs3,rssr*0.1);DSDPCHKERR(info);
507 info=DSDPVecAXPY(dr,dy3,dy);DSDPCHKERR(info);
508 info=DSDPVecSetR(dy,0); DSDPCHKERR(info);
509 info=DSDPApplySMW(M,rhs,dy);DSDPCHKERR(info);
512 DSDPFunctionReturn(0);
516 #define __FUNCT__ "DSDPZeroFixedVariables"
519 FixedVariables *fv=&M.schur->fv;
521 for (i=0;i<fv->nvars;i++){
522 info=DSDPVecSetElement(dy,fv->var[i],0.0);DSDPCHKERR(info);
524 DSDPFunctionReturn(0);
528 #define __FUNCT__ "DSDPApplyFixedVariables"
532 FixedVariables *fv=&M.schur->fv;
533 info=DSDPVecGetC(y,&scl);DSDPCHKERR(info);
535 for (i=0;i<fv->nvars;i++){
536 vv=fv->fval[i]*fabs(scl);
538 info=DSDPVecSetElement(y,jj,vv);DSDPCHKERR(info);
540 DSDPFunctionReturn(0);
544 #define __FUNCT__ "DSDPFixedVariableNorm"
548 FixedVariables *fv=&M.schur->fv;
550 for (i=0;i<fv->nvars;i++){
551 jj=fv->var[i]; vv=fv->fval[i];
552 info=DSDPVecAddC(y,1.0);DSDPCHKERR(info);
553 info=DSDPVecAddElement(y,jj,vv*vv);DSDPCHKERR(info);
555 DSDPFunctionReturn(0);
559 #define __FUNCT__ "DSDPComputeFixedYX"
563 FixedVariables *fv=&M.schur->fv;
565 for (i=0;i<fv->nvars;i++){
567 info=DSDPVecGetElement(berr,jj,&vv);DSDPCHKERR(info);
568 info=DSDPVecSetElement(berr,jj,0);DSDPCHKERR(info);
569 info=DSDPVecAddC(berr,-vv*fv->fval[i]);DSDPCHKERR(info);
570 info=DSDPVecAddR(berr,fabs(vv));DSDPCHKERR(info);
572 if (fv->xout) fv->xout[i]=-vv;
573 DSDPLogInfo(0,2,
"FIXED VAR DUAL: %d %4.4f, ADD %4.4f to objective.\n",jj,vv,-vv*fv->fval[i]);
575 DSDPFunctionReturn(0);
579 #define __FUNCT__ "DSDPIsFixed"
582 FixedVariables *fv=&M.schur->fv;
585 for (i=0;i<fv->nvars;i++){
586 if (fv->var[i]==vari){
591 DSDPFunctionReturn(0);
594 #define __FUNCT__ "DSDPInitializeFixedVariables"
595 int DSDPInitializeFixedVariable( FixedVariables *fv){
602 DSDPFunctionReturn(0);
606 #define __FUNCT__ "DSDPAddFixedVariables"
607 int DSDPAddFixedVariable(
DSDPSchurMat M,
int vari,
double val){
608 int i,t,*iinew,info,nvars;
609 double *ddnew,*vvnew;
610 FixedVariables *fv=&M.schur->fv;
613 if (nvars>=fv->nmaxvars){
615 DSDPCALLOC2(&iinew,
int,t,&info);
616 DSDPCALLOC2(&ddnew,
double,t,&info);
617 DSDPCALLOC2(&vvnew,
double,t,&info);
618 for (i=0;i<nvars;i++){
620 ddnew[i]=fv->fval[i];
621 vvnew[i]=fv->fdual[i];
623 DSDPFREE(&fv->var,&info);DSDPCHKERR(info);
624 DSDPFREE(&fv->fval,&info);DSDPCHKERR(info);
625 DSDPFREE(&fv->fdual,&info);DSDPCHKERR(info);
631 fv->var[fv->nvars]=vari;
632 fv->fval[fv->nvars]=val;
634 DSDPFunctionReturn(0);
640 #define __FUNCT__ "DSDPSparsityInSchurMat"
650 int info,*iptr,m=mm+2;
654 info=DSDPVecZero(R);DSDPCHKERR(info);
655 info=DSDPVecGetArray(R,&dd);DSDPCHKERR(info);
658 memcpy((
void*)rnnz,(
void*)(iptr+1),(mm)*
sizeof(
int));
659 info=DSDPVecRestoreArray(R,&dd);DSDPCHKERR(info);
660 DSDPFunctionReturn(0);
665 #define __FUNCT__ "DSDPSetFixedVariable"
678 DSDPLogInfo(0,2,
"Set Fixed Variable: %d, %12.8f\n",vari,val);
679 info= DSDPAddFixedVariable(dsdp->M,vari,val);DSDPCHKERR(info);
680 DSDPFunctionReturn(0);
684 #define __FUNCT__ "DSDPSetFixedVariables"
698 for (i=0;i<nvars;i++){
700 dsdp->M.schur->fv.xout=xout;
702 DSDPFunctionReturn(0);
706 #define __FUNCT__ "DSDPGetFixedYX"
707 int DSDPGetFixedYX(
DSDP dsdp,
int vari,
double *dd){
709 FixedVariables *fv=&dsdp->M.schur->fv;
711 for (i=0;i<fv->nvars;i++){
712 if (vari==fv->var[i]){
713 *dd=fv->fdual[i];
break;
716 DSDPFunctionReturn(0);