DSDP
|
00001 #include "dsdpsdp.h" 00002 #include "dsdpsys.h" 00003 00004 static int sdpvecvecevent=0,sdpdotevent=0; 00010 #undef __FUNCT__ 00011 #define __FUNCT__ "DSDPBlockASum" 00012 00020 int DSDPBlockASum(DSDPBlockData *ADATA, double aa, DSDPVec Yk, DSDPVMat XX){ 00021 00022 double *xx,ytmp,scl=ADATA->scl; 00023 int ii,vari,n,nn,info; 00024 00025 DSDPFunctionBegin; 00026 info=DSDPVMatGetSize(XX, &n); DSDPCHKERR(info); 00027 info=DSDPVMatGetArray(XX, &xx, &nn); DSDPCHKERR(info); 00028 for (ii=0;ii<ADATA->nnzmats;ii++){ 00029 vari=ADATA->nzmat[ii]; 00030 info=DSDPVecGetElement(Yk,vari,&ytmp);DSDPCHKVARERR(vari,info); 00031 if (ytmp==0) continue; 00032 info = DSDPDataMatAddMultiple(ADATA->A[ii], -aa*scl*ytmp, xx,nn,n); DSDPCHKVARERR(vari,info); 00033 } 00034 info=DSDPVMatRestoreArray(XX, &xx, &nn); DSDPCHKERR(info); 00035 DSDPFunctionReturn(0); 00036 } 00037 00038 #undef __FUNCT__ 00039 #define __FUNCT__ "DSDPBlockADot" 00040 00049 int DSDPBlockADot(DSDPBlockData *ADATA, double aa, DSDPVec Alpha, DSDPVMat X, DSDPVec AX){ 00050 00051 int ii,vari,n,nn,info; 00052 double *x,sum=0,aalpha=0,scl=ADATA->scl; 00053 00054 DSDPFunctionBegin; 00055 DSDPEventLogBegin(sdpdotevent); 00056 info=DSDPVMatScaleDiagonal(X,0.5); DSDPCHKERR(info); 00057 info=DSDPVMatGetSize(X, &n); DSDPCHKERR(info); 00058 info=DSDPVMatGetArray(X, &x, &nn); DSDPCHKERR(info); 00059 for (ii=0;ii<ADATA->nnzmats; ii++){ /* Matrix Entries */ 00060 vari=ADATA->nzmat[ii]; 00061 info=DSDPVecGetElement(Alpha,vari,&aalpha);DSDPCHKVARERR(vari,info); 00062 if (aalpha==0.0) continue; 00063 info=DSDPDataMatDot(ADATA->A[ii],x,nn,n,&sum);DSDPCHKVARERR(vari,info); 00064 info=DSDPVecAddElement(AX,vari,aa*aalpha*sum*scl);DSDPCHKVARERR(vari,info); 00065 } 00066 info=DSDPVMatRestoreArray(X, &x, &nn); DSDPCHKERR(info); 00067 info=DSDPVMatScaleDiagonal(X,2.0); DSDPCHKERR(info); 00068 DSDPEventLogEnd(sdpdotevent); 00069 DSDPFunctionReturn(0); 00070 } 00071 00072 #undef __FUNCT__ 00073 #define __FUNCT__ "DSDPBlockvAv" 00074 00084 int DSDPBlockvAv(DSDPBlockData *ADATA, double aa, DSDPVec Alpha, SDPConeVec V, DSDPVec VAV){ 00085 00086 int ii,vari,info; 00087 double sum=0,aalpha=0,scl=ADATA->scl; 00088 00089 DSDPFunctionBegin; 00090 DSDPEventLogBegin(sdpvecvecevent); 00091 if (aa==0){DSDPFunctionReturn(0);} 00092 for (ii=0;ii<ADATA->nnzmats; ii++){ /* Matrix Entries */ 00093 vari=ADATA->nzmat[ii]; 00094 info=DSDPVecGetElement(Alpha,vari,&aalpha);DSDPCHKVARERR(vari,info); 00095 if (aalpha==0.0) continue; 00096 info=DSDPDataMatVecVec(ADATA->A[ii],V,&sum);DSDPCHKVARERR(vari,info); 00097 info=DSDPVecAddElement(VAV,vari,aa*aalpha*sum*scl);DSDPCHKVARERR(vari,info); 00098 } 00099 DSDPEventLogEnd(sdpvecvecevent); 00100 DSDPFunctionReturn(0); 00101 } 00102 00103 #undef __FUNCT__ 00104 #define __FUNCT__ "DSDPBlockFactorData" 00105 00113 int DSDPBlockFactorData(DSDPBlockData *ADATA, 00114 DSDPVMat X, SDPConeVec W){ 00115 00116 int ii,vari,n,nn,info,*iwork3n,i13,n26; 00117 double *x,*dwork3n; 00118 DSDPDataMat AA; 00119 00120 DSDPFunctionBegin; 00121 info=DSDPVMatGetSize(X, &n); DSDPCHKERR(info); 00122 i13=13*n+1;n26=26*n+1; 00123 DSDPCALLOC2(&dwork3n,double,n26,&info);DSDPCHKERR(info); 00124 DSDPCALLOC2(&iwork3n,int,i13,&info);DSDPCHKERR(info); 00125 info=DSDPVMatGetArray(X, &x, &nn); DSDPCHKERR(info); 00126 for (ii=0;ii<ADATA->nnzmats; ii++){ /* Matrix Entries */ 00127 info=DSDPBlockGetMatrix(ADATA,ii,&vari,0,&AA);DSDPCHKVARERR(vari,info); 00128 DSDPLogInfo(0,39,"SDP Data Mat Setup: %d\n",vari); 00129 if (vari==0) continue; 00130 info=DSDPDataMatFactor(AA,W,x,nn,dwork3n,n26,iwork3n,i13); DSDPCHKVARERR(vari,info); 00131 } 00132 info=DSDPVMatRestoreArray(X, &x, &nn); DSDPCHKERR(info); 00133 DSDPFREE(&dwork3n,&info);DSDPCHKERR(info); 00134 DSDPFREE(&iwork3n,&info);DSDPCHKERR(info); 00135 DSDPFunctionReturn(0); 00136 } 00137 00138 00139 #undef __FUNCT__ 00140 #define __FUNCT__ "DSDPBlockEventZero" 00141 int DSDPBlockEventZero(void){ 00142 DSDPFunctionBegin; 00143 sdpvecvecevent=0;sdpdotevent=0; 00144 DSDPFunctionReturn(0); 00145 } 00146 00147 #undef __FUNCT__ 00148 #define __FUNCT__ "DSDPBlockEventInitialize" 00149 int DSDPBlockEventInitialize(void){ 00150 DSDPFunctionBegin; 00151 if (sdpvecvecevent==0){DSDPEventLogRegister("SDP VecMatVec",&sdpvecvecevent);} 00152 if (sdpdotevent==0){DSDPEventLogRegister("SDP Dot",&sdpdotevent);} 00153 DSDPFunctionReturn(0); 00154 } 00155 00156 #undef __FUNCT__ 00157 #define __FUNCT__ "DSDPBlockDataInitialize" 00158 00163 int DSDPBlockDataInitialize(DSDPBlockData *ADATA){ 00164 DSDPFunctionBegin; 00165 ADATA->nnzmats=0; 00166 ADATA->maxnnzmats=0; 00167 ADATA->nzmat=0; 00168 ADATA->A=0; 00169 ADATA->r=1.0; 00170 ADATA->scl=1.0; 00171 /* ADATA->n=0; */ 00172 DSDPFunctionReturn(0); 00173 } 00174 00175 #undef __FUNCT__ 00176 #define __FUNCT__ "DSDPBlockTakeDownData" 00177 00182 int DSDPBlockTakeDownData(DSDPBlockData *ADATA){ 00183 DSDPFunctionBegin; 00184 DSDPFunctionReturn(0); 00185 } 00186 00187 00188 #undef __FUNCT__ 00189 #define __FUNCT__ "DSDPBlockDataDestroy" 00190 00195 int DSDPBlockDataDestroy(DSDPBlockData *ADATA){ 00196 int ii,vari,info; 00197 DSDPFunctionBegin; 00198 if (!ADATA){DSDPFunctionReturn(0);} 00199 DSDPLogInfo(0,18,"Destroying All Existing Data Matrices \n"); 00200 for (ii=0; ii<ADATA->nnzmats; ii++){ 00201 vari=ADATA->nzmat[ii]; 00202 info = DSDPDataMatDestroy(&ADATA->A[ii]);DSDPCHKVARERR(vari,info); 00203 ADATA->nzmat[ii]=0; 00204 } 00205 ADATA->nnzmats=0; 00206 info=DSDPBlockTakeDownData(ADATA);DSDPCHKERR(info); 00207 DSDPFREE(&ADATA->nzmat,&info);DSDPCHKERR(info); 00208 DSDPFREE(&ADATA->A,&info);DSDPCHKERR(info); 00209 info=DSDPBlockDataInitialize(ADATA);DSDPCHKERR(info); 00210 DSDPFunctionReturn(0); 00211 } 00212 00213 #undef __FUNCT__ 00214 #define __FUNCT__ "DSDPBlockDataAllocate" 00215 00221 int DSDPBlockDataAllocate(DSDPBlockData *ADATA, int nnz){ 00222 int j,info,*nzmat; 00223 DSDPDataMat *A; 00224 DSDPFunctionBegin; 00225 if (!ADATA){DSDPFunctionReturn(0);} 00226 if (nnz<=ADATA->maxnnzmats){DSDPFunctionReturn(0);} 00227 DSDPLogInfo(0,18,"REALLOCATING SPACE FOR %d SDP BLOCK MATRICES! Previously allocated: %d \n",nnz,ADATA->maxnnzmats); 00228 DSDPCALLOC2(&A,struct DSDPDataMat_C,nnz,&info);DSDPCHKERR(info); 00229 DSDPCALLOC2(&nzmat,int,nnz,&info);DSDPCHKERR(info); 00230 for (j=0;j<nnz;j++){nzmat[j]=0;} 00231 for (j=0;j<nnz;j++){info = DSDPDataMatInitialize(&A[j]);DSDPCHKERR(info);} 00232 if (ADATA->maxnnzmats>0){ 00233 for (j=0;j<ADATA->nnzmats;j++){nzmat[j]=ADATA->nzmat[j];} 00234 for (j=0;j<ADATA->nnzmats;j++){A[j]=ADATA->A[j];} 00235 DSDPFREE(&ADATA->A,&info);DSDPCHKERR(info); 00236 DSDPFREE(&ADATA->nzmat,&info);DSDPCHKERR(info); 00237 } else { 00238 ADATA->nnzmats=0; 00239 } 00240 ADATA->maxnnzmats=nnz; 00241 ADATA->nzmat=nzmat; 00242 ADATA->A=A; 00243 DSDPFunctionReturn(0); 00244 } 00245 00246 #undef __FUNCT__ 00247 #define __FUNCT__ "DSDPBlockDataMarkNonzeroMatrices" 00248 00254 int DSDPBlockDataMarkNonzeroMatrices(DSDPBlockData *ADATA,int*annz){ 00255 int i; 00256 DSDPFunctionBegin; 00257 for (i=0; i<ADATA->nnzmats; i++){ 00258 annz[ADATA->nzmat[i]]++; 00259 } 00260 DSDPFunctionReturn(0); 00261 } 00262 00263 #undef __FUNCT__ 00264 #define __FUNCT__ "DSDPBlockCountNonzerosMatrices" 00265 00272 int DSDPBlockCountNonzeroMatrices(DSDPBlockData *ADATA,int*nzmats){ 00273 DSDPFunctionBegin; 00274 *nzmats=ADATA->nnzmats; 00275 DSDPFunctionReturn(0); 00276 } 00277 00278 #undef __FUNCT__ 00279 #define __FUNCT__ "DSDPBlockDataRank" 00280 int DSDPBlockDataRank(DSDPBlockData *ADATA, int *trank, int n){ 00281 int ii,vari,info,ri,r2=0; 00282 DSDPDataMat AA; 00283 00284 DSDPFunctionBegin; 00285 for (ii=0;ii<ADATA->nnzmats;ii++){ 00286 info=DSDPBlockGetMatrix(ADATA,ii,&vari,0,&AA);DSDPCHKVARERR(vari,info); 00287 if (vari==0) continue; 00288 info=DSDPDataMatGetRank(AA,&ri,n); DSDPCHKVARERR(vari,info); 00289 r2+=ri; 00290 } 00291 *trank=r2; 00292 DSDPFunctionReturn(0); 00293 } 00294 00295 #undef __FUNCT__ 00296 #define __FUNCT__ "DSDPBlockGetMatrix" 00297 00307 int DSDPBlockGetMatrix(DSDPBlockData *ADATA,int id, int *vari, double *scl, DSDPDataMat *A){ 00308 DSDPFunctionBegin; 00309 if (id>=0 && id < ADATA->nnzmats){ 00310 if (vari) *vari=ADATA->nzmat[id]; 00311 if (scl) *scl=ADATA->scl; 00312 if (A) *A=ADATA->A[id]; 00313 } else { 00314 DSDPSETERR2(2,"Invalid Matrix request. 0 <= %d < %d\n",id,ADATA->nnzmats); 00315 } 00316 DSDPFunctionReturn(0); 00317 } 00318 00319 #undef __FUNCT__ 00320 #define __FUNCT__ "DSDPBlockDataRowSparsity" 00321 00330 int DSDPBlockDataRowSparsity(DSDPBlockData *ADATA,int row, int ai[], int rnnz[],int n){ 00331 int info,i,vari,rn; 00332 DSDPFunctionBegin; 00333 if (ai){ 00334 for (i=0; i<ADATA->nnzmats; i++){ 00335 vari=ADATA->nzmat[i]; 00336 if (ai[vari]==0){continue;} 00337 info=DSDPDataMatGetRowNonzeros(ADATA->A[i],row, n, rnnz, &rn); DSDPCHKVARERR(vari,info); 00338 } 00339 } 00340 DSDPFunctionReturn(0); 00341 } 00342 00343 #undef __FUNCT__ 00344 #define __FUNCT__ "DSDPBlockRemoveDataMatrix" 00345 00351 int DSDPBlockRemoveDataMatrix(DSDPBlockData *ADATA,int vari){ 00352 int info,ii,k; 00353 DSDPFunctionBegin; 00354 for (ii=0;ii<ADATA->nnzmats;ii++){ 00355 if (ADATA->nzmat[ii]==vari){ 00356 info=DSDPDataMatDestroy(&ADATA->A[ii]);DSDPCHKVARERR(vari,info); 00357 info=DSDPSetDataMatZero(&ADATA->A[ii]);DSDPCHKVARERR(vari,info); 00358 for (k=ii;k<ADATA->nnzmats;k++){ 00359 ADATA->A[k]=ADATA->A[k+1]; 00360 ADATA->nzmat[k]=ADATA->nzmat[k+1]; 00361 } 00362 ADATA->nnzmats--; 00363 info=DSDPSetDataMatZero(&ADATA->A[ADATA->nnzmats]);DSDPCHKERR(info); 00364 DSDPFunctionReturn(0); 00365 } 00366 } 00367 DSDPFunctionReturn(0); 00368 } 00369 00370 #undef __FUNCT__ 00371 #define __FUNCT__ "DSDPBlockAddDataMatrix" 00372 00381 int DSDPBlockAddDataMatrix(DSDPBlockData *ADATA,int vari, struct DSDPDataMat_Ops* dsdpdataops, void* data){ 00382 int info,ii; 00383 DSDPFunctionBegin; 00384 if (ADATA->nnzmats>=ADATA->maxnnzmats){ 00385 info=DSDPBlockDataAllocate(ADATA,2*ADATA->maxnnzmats+7);DSDPCHKERR(info); 00386 } 00387 ii=ADATA->nnzmats; 00388 info=DSDPDataMatDestroy(&ADATA->A[ii]);DSDPCHKERR(info); 00389 info=DSDPDataMatSetData(&ADATA->A[ii], dsdpdataops, data);DSDPCHKVARERR(vari,info); 00390 ADATA->nzmat[ii]=vari; 00391 ADATA->nnzmats++; 00392 DSDPFunctionReturn(0); 00393 } 00394 00395 #undef __FUNCT__ 00396 #define __FUNCT__ "DSDPBlockSetDataMatrix" 00397 00406 int DSDPBlockSetDataMatrix(DSDPBlockData *ADATA,int vari, struct DSDPDataMat_Ops* dsdpdataops, void* data){ 00407 int info; 00408 DSDPFunctionBegin; 00409 info=DSDPBlockRemoveDataMatrix(ADATA,vari);DSDPCHKERR(info); 00410 info=DSDPBlockAddDataMatrix(ADATA,vari,dsdpdataops,data);DSDPCHKERR(info); 00411 DSDPFunctionReturn(0); 00412 } 00413 00414 #undef __FUNCT__ 00415 #define __FUNCT__ "DSDPBlockNorm2" 00416 int DSDPBlockNorm2(DSDPBlockData *ADATA, int n){ 00417 double fn2,tt=0; 00418 int ii,info; 00419 DSDPFunctionBegin; 00420 for (ii=0;ii<ADATA->nnzmats;ii++){ 00421 info=DSDPDataMatFNorm2(ADATA->A[ii],n,&fn2); DSDPCHKERR(info); 00422 tt+=fn2; 00423 } 00424 DSDPFunctionReturn(0); 00425 } 00426 00427 #undef __FUNCT__ 00428 #define __FUNCT__ "DSDPBlockANorm2" 00429 int DSDPBlockANorm2(DSDPBlockData *ADATA, DSDPVec ANorm2, int n){ 00430 00431 double fn2,scl=ADATA->scl; 00432 int ii,vari,info; 00433 00434 DSDPFunctionBegin; 00435 info=DSDPBlockNorm2(ADATA,n);DSDPCHKERR(info); 00436 scl=ADATA->scl; 00437 for (ii=0;ii<ADATA->nnzmats;ii++){ 00438 vari=ADATA->nzmat[ii]; 00439 info=DSDPDataMatFNorm2(ADATA->A[ii],n,&fn2); DSDPCHKVARERR(vari,info); 00440 info=DSDPVecAddElement(ANorm2,vari,fn2*scl);DSDPCHKVARERR(vari,info); 00441 } 00442 DSDPFunctionReturn(0); 00443 } 00444 00445 00446 #undef __FUNCT__ 00447 #define __FUNCT__ "DSDPBlockView" 00448 00454 int DSDPBlockView(DSDPBlockData *ADATA){ 00455 int ii,kk; 00456 00457 DSDPFunctionBegin; 00458 for (ii=0;ii<ADATA->nnzmats;ii++){ 00459 kk=ADATA->nzmat[ii]; 00460 if (kk==0){ printf("+ C\n");} 00461 else { printf(" - A[%d] y%d\n",kk,kk);} 00462 } 00463 printf(" = S >= 0\n"); 00464 DSDPFunctionReturn(0); 00465 } 00466 #undef __FUNCT__ 00467 #define __FUNCT__ "DSDPBlockView2" 00468 00474 int DSDPBlockView2(DSDPBlockData *ADATA){ 00475 int ii,vari,info; 00476 00477 DSDPFunctionBegin; 00478 for (ii=0;ii<ADATA->nnzmats;ii++){ 00479 vari=ADATA->nzmat[ii]; 00480 printf("A[%d] y%d \n",vari,vari); 00481 info=DSDPDataMatView(ADATA->A[ii]); DSDPCHKERR(info); 00482 } 00483 DSDPFunctionReturn(0); 00484 } 00485 00486 00487 #undef __FUNCT__ 00488 #define __FUNCT__ "DSDPDataMatCheck" 00489 00498 int DSDPDataMatCheck(DSDPDataMat AA, SDPConeVec W, DSDPIndex IS, DSDPVMat XX){ 00499 00500 double *xx,ack,vAv=0,esum=0,sum,eignorm,fnorm22,dnorm,scl=1; 00501 int k,n,nn,rank,info; 00502 00503 DSDPFunctionBegin; 00504 info=SDPConeVecGetSize(W,&n);DSDPCHKERR(info); 00505 00506 info=DSDPVMatZeroEntries(XX);DSDPCHKERR(info); 00507 info=DSDPDataMatGetRank(AA,&rank,n);DSDPCHKERR(info); 00508 for (k=0; k<rank; k++){ 00509 info=DSDPDataMatGetEig(AA,k,W,IS,&ack); DSDPCHKERR(info); 00510 info=SDPConeVecDot(W,W,&eignorm);DSDPCHKERR(info); 00511 info=DSDPVMatAddOuterProduct(XX,scl*ack,W);DSDPCHKERR(info); 00512 info=DSDPDataMatVecVec(AA,W,&sum);DSDPCHKERR(info); 00513 vAv+=ack*ack*eignorm*eignorm*scl; 00514 } 00515 info=DSDPDataMatFNorm2(AA,n,&fnorm22); DSDPCHKERR(info); 00516 00517 info=DSDPVMatScaleDiagonal(XX,0.5); DSDPCHKERR(info); 00518 info=DSDPVMatGetArray(XX, &xx, &nn); DSDPCHKERR(info); 00519 info=DSDPDataMatDot(AA,xx,nn,n,&esum);DSDPCHKERR(info); 00520 info=DSDPVMatRestoreArray(XX, &xx, &nn); DSDPCHKERR(info); 00521 info=DSDPVMatScaleDiagonal(XX,2.0); DSDPCHKERR(info); 00522 00523 info=DSDPVMatGetArray(XX, &xx, &nn); DSDPCHKERR(info); 00524 info=DSDPDataMatAddMultiple(AA,-1.0,xx,nn,n); DSDPCHKERR(info); 00525 info=DSDPVMatRestoreArray(XX, &xx, &nn); DSDPCHKERR(info); 00526 if (0==1){info=DSDPVMatView(XX);DSDPCHKERR(info);} 00527 info=DSDPVMatNormF2(XX,&dnorm); DSDPCHKERR(info); 00528 printf(" %4.4e, %4.4e %4.4e\n",esum,vAv,fnorm22); 00529 printf(" error1: %4.4e, error2: %4.4e, error3: %4.4e\n",sqrt(dnorm),fabs(esum-vAv),fabs(fnorm22-vAv)); 00530 if (dnorm>1) printf("Check Add or eigs\n"); 00531 if (fabs(esum-vAv) > 1.0) printf("Check vAv \n"); 00532 if (fabs(fnorm22-vAv) > 1.0) printf("Check fnorm22\n"); 00533 00534 DSDPFunctionReturn(0); 00535 } 00536