DSDP
|
00001 #include "dsdpsys.h" 00002 #include "dsdpvec.h" 00003 #include "dsdplapack.h" 00004 00010 typedef struct{ 00011 char UPLO; 00012 double *val; 00013 double *v2; 00014 double *sscale; 00015 int scaleit; 00016 int n; 00017 int owndata; 00018 } dtpumat; 00019 00020 static const char* lapackname="DENSE,SYMMETRIC,PACKED STORAGE"; 00021 00022 int DTPUMatEigs(void*AA,double W[],double IIWORK[], int nn1, double *mineig){ 00023 dtpumat* AAA=(dtpumat*) AA; 00024 ffinteger info,INFO=0,M,N=AAA->n; 00025 ffinteger IL=1,IU=1,LDZ=1,IFAIL; 00026 ffinteger *IWORK=(ffinteger*)IIWORK; 00027 double *AP=AAA->val,ABSTOL=1e-13; 00028 double Z=0,VL=-1e10,VU=1; 00029 double *WORK; 00030 char UPLO=AAA->UPLO,JOBZ='N',RANGE='I'; 00031 00032 DSDPCALLOC2(&WORK,double,7*N,&info);DSDPCHKERR(info); 00033 DSDPCALLOC2(&IWORK,ffinteger,5*N,&info);DSDPCHKERR(info); 00034 dspevx(&JOBZ,&RANGE,&UPLO,&N,AP,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,&Z,&LDZ,WORK,IWORK,&IFAIL,&INFO); 00035 00036 /* 00037 DSDPCALLOC2(&WORK,double,2*N,&info); 00038 LWORK=2*N; 00039 dspevd(&JOBZ,&UPLO,&N,AP,W,&Z,&LDZ,WORK,&LWORK,IWORK,&LIWORK,&INFO); 00040 */ 00041 *mineig=W[0]; 00042 DSDPFREE(&WORK,&info);DSDPCHKERR(info); 00043 DSDPFREE(&IWORK,&info);DSDPCHKERR(info); 00044 return INFO; 00045 } 00046 00047 00048 #undef __FUNCT__ 00049 #define __FUNCT__ "DSDPLAPACKROUTINE" 00050 static void dtpuscalevec(double alpha, double v1[], double v2[], double v3[], int n){ 00051 int i; 00052 for (i=0;i<n;i++){ 00053 v3[i] = (alpha*v1[i]*v2[i]); 00054 } 00055 } 00056 00057 static void dtpuscalemat(double vv[], double ss[], int n){ 00058 int i; 00059 for (i=0;i<n;i++,vv+=i){ 00060 dtpuscalevec(ss[i],vv,ss,vv,i+1); 00061 } 00062 } 00063 00064 static int DTPUMatCreateWData(int n,double nz[], int nnz, dtpumat**M){ 00065 int i,nn=(n*n+n)/2,info; 00066 double dtmp; 00067 dtpumat*M23; 00068 if (nnz<nn){DSDPSETERR1(2,"Array must have length of : %d \n",nn);} 00069 for (i=0;i<nnz;i++) dtmp=nz[i]; 00070 DSDPCALLOC1(&M23,dtpumat,&info);DSDPCHKERR(info); 00071 DSDPCALLOC2(&M23->sscale,double,n,&info);DSDPCHKERR(info); 00072 M23->owndata=0; M23->val=nz; M23->n=n; M23->UPLO='U'; 00073 for (i=0;i<n;i++)M23->sscale[i]=1.0; 00074 M23->scaleit=0; 00075 *M=M23; 00076 return 0; 00077 } 00078 00079 00080 00081 #undef __FUNCT__ 00082 #define __FUNCT__ "DSDPLAPACK ROUTINE" 00083 00084 00085 static int DTPUMatMult(void* AA, double x[], double y[], int n){ 00086 dtpumat* A=(dtpumat*) AA; 00087 ffinteger ione=1,N=n; 00088 double BETA=0.0,ALPHA=1.0; 00089 double *AP=A->val,*Y=y,*X=x; 00090 char UPLO=A->UPLO; 00091 00092 if (A->n != n) return 1; 00093 if (x==0 && n>0) return 3; 00094 dspmv(&UPLO,&N,&ALPHA,AP,X,&ione,&BETA,Y,&ione); 00095 return 0; 00096 } 00097 00098 static int DTPUMatSolve(void* AA, double b[], double x[], int n){ 00099 dtpumat* A=(dtpumat*) AA; 00100 ffinteger INFO,NRHS=1,LDB=A->n,N=A->n; 00101 double *AP=A->val,*ss=A->sscale; 00102 char UPLO=A->UPLO; 00103 00104 if (N>0) LDB=N; 00105 dtpuscalevec(1.0,ss,b,x,n); 00106 dpptrs(&UPLO, &N, &NRHS, AP, x, &LDB, &INFO ); 00107 dtpuscalevec(1.0,ss,x,x,n); 00108 return INFO; 00109 } 00110 00111 static int DTPUMatCholeskyFactor(void* AA, int *flag){ 00112 dtpumat* A=(dtpumat*) AA; 00113 int i; 00114 ffinteger INFO,LDA=1,N=A->n; 00115 double *AP=A->val,*ss=A->sscale,*v; 00116 char UPLO=A->UPLO; 00117 00118 if (N<=0) LDA=1; 00119 else LDA=N; 00120 if (A->scaleit){ 00121 for (v=AP,i=0;i<N;i++){ ss[i]=*v;v+=(i+2);} 00122 for (i=0;i<N;i++){ ss[i]=1.0/sqrt(fabs(ss[i])+1.0e-8); } 00123 dtpuscalemat(AP,ss,N); 00124 } 00125 dpptrf(&UPLO, &N, AP, &INFO ); 00126 *flag=INFO; 00127 return 0; 00128 } 00129 00130 static int DTPUMatShiftDiagonal(void* AA, double shift){ 00131 dtpumat* A=(dtpumat*) AA; 00132 int i,n=A->n; 00133 double *v=A->val; 00134 for (i=0; i<n; i++){ 00135 *v+=shift; 00136 v+=i+2; 00137 } 00138 return 0; 00139 } 00140 00141 00142 #undef __FUNCT__ 00143 #define __FUNCT__ "DTPUMatAssemble" 00144 static int DTPUMatAssemble(void*M){ 00145 int info; 00146 double shift=1.0e-15; 00147 DSDPFunctionBegin; 00148 info= DTPUMatShiftDiagonal(M, shift); DSDPCHKERR(info); 00149 DSDPFunctionReturn(0); 00150 } 00151 00152 static int DTPUMatRowNonzeros(void*M, int row, double cols[], int *ncols,int nrows){ 00153 int i; 00154 DSDPFunctionBegin; 00155 *ncols = row+1; 00156 for (i=0;i<=row;i++){ 00157 cols[i]=1.0; 00158 } 00159 for (i=row+1;i<nrows;i++){ 00160 cols[i]=0.0; 00161 } 00162 DSDPFunctionReturn(0); 00163 } 00164 00165 00166 #undef __FUNCT__ 00167 #define __FUNCT__ "DTPUMatDiag" 00168 static int DTPUMatDiag(void*M, int row, double dd){ 00169 int ii; 00170 dtpumat* ABA=(dtpumat*)M; 00171 DSDPFunctionBegin; 00172 ii=row*(row+1)/2+row; 00173 ABA->val[ii] +=dd; 00174 DSDPFunctionReturn(0); 00175 } 00176 #undef __FUNCT__ 00177 #define __FUNCT__ "DTPUMatDiag2" 00178 static int DTPUMatDiag2(void*M, double diag[], int m){ 00179 int row,ii; 00180 dtpumat* ABA=(dtpumat*)M; 00181 DSDPFunctionBegin; 00182 for (row=0;row<m;row++){ 00183 ii=row*(row+1)/2+row; 00184 ABA->val[ii] +=diag[row]; 00185 } 00186 DSDPFunctionReturn(0); 00187 } 00188 00189 static int DTPUMatAddRow(void* AA, int nrow, double dd, double row[], int n){ 00190 dtpumat* A=(dtpumat*) AA; 00191 ffinteger ione=1,nn,nnn; 00192 double *vv=A->val; 00193 00194 nnn=nrow*(nrow+1)/2; 00195 nn=nrow+1; 00196 daxpy(&nn,&dd,row,&ione,vv+nnn,&ione); 00197 00198 return 0; 00199 } 00200 00201 static int DTPUMatZero(void* AA){ 00202 dtpumat* A=(dtpumat*) AA; 00203 int mn=A->n*(A->n+1)/2; 00204 double *vv=A->val; 00205 memset((void*)vv,0,mn*sizeof(double)); 00206 return 0; 00207 } 00208 00209 static int DTPUMatGetSize(void *AA, int *n){ 00210 dtpumat* A=(dtpumat*) AA; 00211 *n=A->n; 00212 return 0; 00213 } 00214 00215 static int DTPUMatDestroy(void* AA){ 00216 int info; 00217 dtpumat* A=(dtpumat*) AA; 00218 if (A && A->owndata){DSDPFREE(&A->val,&info);DSDPCHKERR(info);} 00219 if (A && A->sscale) {DSDPFREE(&A->sscale,&info);DSDPCHKERR(info);} 00220 if (A) {DSDPFREE(&A,&info);DSDPCHKERR(info);} 00221 return 0; 00222 } 00223 00224 static int DTPUMatView(void* AA){ 00225 dtpumat* M=(dtpumat*) AA; 00226 int i,j,kk=0; 00227 double *val=M->val; 00228 for (i=0; i<M->n; i++){ 00229 for (j=0; j<=i; j++){ 00230 printf(" %9.2e",val[kk]); 00231 kk++; 00232 } 00233 printf("\n"); 00234 } 00235 return 0; 00236 } 00237 00238 00239 #include "dsdpschurmat_impl.h" 00240 #include "dsdpsys.h" 00241 static struct DSDPSchurMat_Ops dsdpmmatops; 00242 00243 static int DSDPInitSchurOps(struct DSDPSchurMat_Ops* mops){ 00244 int info; 00245 DSDPFunctionBegin; 00246 info=DSDPSchurMatOpsInitialize(mops);DSDPCHKERR(info); 00247 mops->matrownonzeros=DTPUMatRowNonzeros; 00248 mops->matscaledmultiply=DTPUMatMult; 00249 mops->mataddrow=DTPUMatAddRow; 00250 mops->mataddelement=DTPUMatDiag; 00251 mops->matadddiagonal=DTPUMatDiag2; 00252 mops->matshiftdiagonal=DTPUMatShiftDiagonal; 00253 mops->matassemble=DTPUMatAssemble; 00254 mops->matfactor=DTPUMatCholeskyFactor; 00255 mops->matsolve=DTPUMatSolve; 00256 mops->matdestroy=DTPUMatDestroy; 00257 mops->matzero=DTPUMatZero; 00258 mops->matview=DTPUMatView; 00259 mops->id=1; 00260 mops->matname=lapackname; 00261 DSDPFunctionReturn(0); 00262 } 00263 00264 #undef __FUNCT__ 00265 #define __FUNCT__ "DSDPGetLAPACKPUSchurOps" 00266 int DSDPGetLAPACKPUSchurOps(int n,struct DSDPSchurMat_Ops** sops, void** mdata){ 00267 int info,nn=n*(n+1)/2; 00268 double *vv; 00269 dtpumat*AA; 00270 DSDPFunctionBegin; 00271 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info); 00272 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info); 00273 AA->owndata=1; 00274 AA->scaleit=1; 00275 info=DSDPInitSchurOps(&dsdpmmatops); DSDPCHKERR(info); 00276 *sops=&dsdpmmatops; 00277 *mdata=(void*)AA; 00278 DSDPFunctionReturn(0); 00279 } 00280 00281 00282 static void daddrow(double *v, double alpha, int i, double row[], int n){ 00283 double *s1; 00284 ffinteger j,nn=n,ione=1; 00285 nn=i+1; s1=v+i*(i+1)/2; 00286 daxpy(&nn,&alpha,s1,&ione,row,&ione); 00287 for (j=i+1;j<n;j++){ 00288 s1+=j; 00289 row[j]+=alpha*s1[i]; 00290 } 00291 return; 00292 } 00293 00294 static int DTPUMatInverseMult(void* AA, int indx[], int nind, double x[], double y[], int n){ 00295 dtpumat* A=(dtpumat*) AA; 00296 ffinteger ione=1,N=n; 00297 double BETA=0.0,ALPHA=1.0; 00298 double *AP=A->v2,*Y=y,*X=x; 00299 int i,ii; 00300 char UPLO=A->UPLO; 00301 00302 if (A->n != n) return 1; 00303 if (x==0 && n>0) return 3; 00304 00305 if (nind<n/4 ){ 00306 memset((void*)y,0,n*sizeof(double)); 00307 for (ii=0;ii<nind;ii++){ 00308 i=indx[ii]; ALPHA=x[i]; 00309 daddrow(AP,ALPHA,i,y,n); 00310 } 00311 } else { 00312 ALPHA=1.0; 00313 dspmv(&UPLO,&N,&ALPHA,AP,X,&ione,&BETA,Y,&ione); 00314 } 00315 return 0; 00316 } 00317 00318 00319 static int DTPUMatCholeskyBackward(void* AA, double b[], double x[], int n){ 00320 dtpumat* M=(dtpumat*) AA; 00321 ffinteger N=M->n,INCX=1; 00322 double *AP=M->val,*ss=M->sscale; 00323 char UPLO=M->UPLO,TRANS='N',DIAG='N'; 00324 memcpy(x,b,N*sizeof(double)); 00325 dtpsv(&UPLO,&TRANS,&DIAG, &N, AP, x, &INCX ); 00326 dtpuscalevec(1.0,ss,x,x,n); 00327 return 0; 00328 } 00329 00330 00331 static int DTPUMatCholeskyForward(void* AA, double b[], double x[], int n){ 00332 dtpumat* M=(dtpumat*) AA; 00333 ffinteger N=M->n,INCX=1; 00334 double *AP=M->val,*ss=M->sscale; 00335 char UPLO=M->UPLO,TRANS='T',DIAG='N'; 00336 dtpuscalevec(1.0,ss,b,x,n); 00337 dtpsv(&UPLO,&TRANS,&DIAG, &N, AP, x, &INCX); 00338 return 0; 00339 } 00340 00341 static int DenseSymPSDCholeskyForwardMultiply(void* AA, double x[], double y[], int n){ 00342 dtpumat* A=(dtpumat*) AA; 00343 int i,j,k=0; 00344 ffinteger N=A->n; 00345 double *AP=A->val,*ss=A->sscale; 00346 00347 if (x==0 && N>0) return 3; 00348 for (i=0;i<N;i++){ 00349 for (j=0;j<=i;j++){ 00350 y[i]+=AP[k]*x[j]; 00351 k++; 00352 } 00353 } 00354 for (i=0;i<N;i++){ y[i]=y[i]/ss[i];} 00355 return 0; 00356 } 00357 00358 static int DTPUMatLogDet(void* AA, double *dd){ 00359 dtpumat* A=(dtpumat*) AA; 00360 int i,n=A->n; 00361 double d=0,*v=A->val,*ss=A->sscale; 00362 for (i=0; i<n; i++){ 00363 if (*v<=0) return 1; 00364 d+=2*log(*v/ss[i]); 00365 v+=i+2; 00366 } 00367 *dd=d; 00368 return 0; 00369 } 00370 00371 static int DTPUMatInvert(void* AA){ 00372 dtpumat* A=(dtpumat*) AA; 00373 ffinteger INFO,N=A->n,nn=N*(N+1)/2; 00374 double *v=A->val,*AP=A->v2,*ss=A->sscale; 00375 char UPLO=A->UPLO; 00376 memcpy((void*)AP,(void*)v,nn*sizeof(double)); 00377 dpptri(&UPLO, &N, AP, &INFO ); 00378 if (INFO){ 00379 INFO=DTPUMatShiftDiagonal(AA,1e-11); 00380 memcpy((void*)AP,(void*)v,nn*sizeof(double)); 00381 dpptri(&UPLO, &N, AP, &INFO ); 00382 } 00383 if (A->scaleit){ 00384 dtpuscalemat(AP,ss,N); 00385 } 00386 return INFO; 00387 } 00388 00389 static int DTPUMatInverseAdd(void* AA, double alpha, double y[], int nn, int n){ 00390 dtpumat* A=(dtpumat*) AA; 00391 ffinteger N=nn,ione=1; 00392 double *v2=A->v2; 00393 daxpy(&N,&alpha,v2,&ione,y,&ione); 00394 return 0; 00395 } 00396 00397 00398 static int DTPUMatScaleDiagonal(void* AA, double dd){ 00399 dtpumat* A=(dtpumat*) AA; 00400 int i,n=A->n; 00401 double *v=A->val; 00402 for (i=0; i<n; i++){ 00403 *v*=dd; 00404 v+=i+2; 00405 } 00406 return 0; 00407 } 00408 00409 static int DTPUMatOuterProduct(void* AA, double alpha, double x[], int n){ 00410 dtpumat* A=(dtpumat*) AA; 00411 ffinteger ione=1,N=n; 00412 double *v=A->val; 00413 char UPLO=A->UPLO; 00414 dspr(&UPLO,&N,&alpha,x,&ione,v); 00415 return 0; 00416 } 00417 00418 00419 static int DenseSymPSDNormF2(void* AA, int n, double *dddot){ 00420 dtpumat* A=(dtpumat*) AA; 00421 ffinteger ione=1,nn=A->n*(A->n+1)/2; 00422 double dd,tt=sqrt(0.5),*val=A->val; 00423 int info; 00424 info=DTPUMatScaleDiagonal(AA,tt); 00425 dd=dnrm2(&nn,val,&ione); 00426 info=DTPUMatScaleDiagonal(AA,1.0/tt); 00427 *dddot=dd*dd*2; 00428 return 0; 00429 } 00430 00431 00432 /* 00433 static int DTPUMatFNorm2(void* AA, double *mnorm){ 00434 dtpumat* A=(dtpumat*) AA; 00435 ffinteger ione=1,n; 00436 double vv=0,*AP=A->val; 00437 n=A->n*(A->n+1)/2; 00438 vv=dnrm2(&n,AP,&ione); 00439 *mnorm=2.0*vv; 00440 return 1; 00441 } 00442 */ 00443 00444 #include "dsdpdualmat_impl.h" 00445 #include "dsdpdatamat_impl.h" 00446 #include "dsdpxmat_impl.h" 00447 #include "dsdpdsmat_impl.h" 00448 00449 00450 00451 static int DTPUMatFull(void*A, int*full){ 00452 *full=1; 00453 return 0; 00454 } 00455 00456 00457 static int DTPUMatGetDenseArray(void* A, double *v[], int*n){ 00458 dtpumat* ABA=(dtpumat*)A; 00459 *v=ABA->val; 00460 *n=(ABA->n)*(ABA->n+1)/2; 00461 return 0; 00462 } 00463 00464 static int DTPUMatRestoreDenseArray(void* A, double *v[], int *n){ 00465 *v=0;*n=0; 00466 return 0; 00467 } 00468 00469 static int DDenseSetXMat(void*A, double v[], int nn, int n){ 00470 double *vv; 00471 dtpumat* ABA=(dtpumat*)A; 00472 vv=ABA->val; 00473 if (v!=vv){ 00474 memcpy((void*)vv,(void*)v,nn*sizeof(double)); 00475 } 00476 return 0; 00477 } 00478 00479 static int DDenseVecVec(void* A, double x[], int n, double *v){ 00480 dtpumat* ABA=(dtpumat*)A; 00481 int i,j,k=0; 00482 double dd=0,*val=ABA->val; 00483 *v=0.0; 00484 for (i=0; i<n; i++){ 00485 for (j=0;j<i;j++){ 00486 dd+=2*x[i]*x[j]*val[k]; 00487 k++; 00488 } 00489 dd+=x[i]*x[i]*val[k]; 00490 k++; 00491 } 00492 *v=dd; 00493 return 0; 00494 } 00495 00496 static struct DSDPDSMat_Ops tdsdensematops; 00497 static int DSDPDSDenseInitializeOps(struct DSDPDSMat_Ops* densematops){ 00498 int info; 00499 if (!densematops) return 0; 00500 info=DSDPDSMatOpsInitialize(densematops); DSDPCHKERR(info); 00501 densematops->matseturmat=DDenseSetXMat; 00502 densematops->matview=DTPUMatView; 00503 densematops->matdestroy=DTPUMatDestroy; 00504 densematops->matgetsize=DTPUMatGetSize; 00505 densematops->matzeroentries=DTPUMatZero; 00506 densematops->matmult=DTPUMatMult; 00507 densematops->matvecvec=DDenseVecVec; 00508 densematops->id=1; 00509 densematops->matname=lapackname; 00510 return 0; 00511 } 00512 00513 #undef __FUNCT__ 00514 #define __FUNCT__ "DSDPCreateDSMatWithArray" 00515 int DSDPCreateDSMatWithArray(int n,double vv[], int nnz, struct DSDPDSMat_Ops* *dsmatops, void**dsmat){ 00516 int info; 00517 dtpumat*AA; 00518 DSDPFunctionBegin; 00519 info=DTPUMatCreateWData(n,vv,nnz,&AA); DSDPCHKERR(info); 00520 AA->owndata=0; 00521 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info); 00522 *dsmatops=&tdsdensematops; 00523 *dsmat=(void*)AA; 00524 DSDPFunctionReturn(0); 00525 } 00526 00527 00528 #undef __FUNCT__ 00529 #define __FUNCT__ "DSDPCreateDSMat" 00530 int DSDPCreateDSMat(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){ 00531 int info,nn=n*(n+1)/2; 00532 double *vv; 00533 dtpumat*AA; 00534 DSDPFunctionBegin; 00535 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info); 00536 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info); 00537 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info); 00538 *dsmatops=&tdsdensematops; 00539 *dsmat=(void*)AA; 00540 AA->owndata=1; 00541 DSDPFunctionReturn(0); 00542 } 00543 00544 static struct DSDPVMat_Ops turdensematops; 00545 00546 static int DSDPDenseXInitializeOps(struct DSDPVMat_Ops* densematops){ 00547 int info; 00548 if (!densematops) return 0; 00549 info=DSDPVMatOpsInitialize(densematops); DSDPCHKERR(info); 00550 densematops->matview=DTPUMatView; 00551 densematops->matscalediagonal=DTPUMatScaleDiagonal; 00552 densematops->matshiftdiagonal=DTPUMatShiftDiagonal; 00553 densematops->mataddouterproduct=DTPUMatOuterProduct; 00554 densematops->matdestroy=DTPUMatDestroy; 00555 densematops->matfnorm2=DenseSymPSDNormF2; 00556 densematops->matgetsize=DTPUMatGetSize; 00557 densematops->matzeroentries=DTPUMatZero; 00558 densematops->matgeturarray=DTPUMatGetDenseArray; 00559 densematops->matrestoreurarray=DTPUMatRestoreDenseArray; 00560 densematops->matmineig=DTPUMatEigs; 00561 densematops->matmult=DTPUMatMult; 00562 densematops->id=1; 00563 densematops->matname=lapackname; 00564 return 0; 00565 } 00566 00567 #undef __FUNCT__ 00568 #define __FUNCT__ "DSDPXMatCreate" 00569 int DSDPXMatPCreate(int n,struct DSDPVMat_Ops* *xops, void * *xmat){ 00570 int info,nn=n*(n+1)/2; 00571 double *vv; 00572 dtpumat*AA; 00573 DSDPFunctionBegin; 00574 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info); 00575 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info); 00576 AA->owndata=1; 00577 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info); 00578 *xops=&turdensematops; 00579 *xmat=(void*)AA; 00580 DSDPFunctionReturn(0); 00581 } 00582 00583 #undef __FUNCT__ 00584 #define __FUNCT__ "DSDPXMatCreate" 00585 int DSDPXMatPCreateWithData(int n,double nz[],int nnz,struct DSDPVMat_Ops* *xops, void * *xmat){ 00586 int i,info; 00587 double dtmp; 00588 dtpumat*AA; 00589 DSDPFunctionBegin; 00590 for (i=0;i<n*(n+1)/2;i++) dtmp=nz[i]; 00591 info=DTPUMatCreateWData(n,nz,nnz,&AA); DSDPCHKERR(info); 00592 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info); 00593 *xops=&turdensematops; 00594 *xmat=(void*)AA; 00595 DSDPFunctionReturn(0); 00596 } 00597 00598 00599 static struct DSDPDualMat_Ops sdmatops; 00600 static int SDualOpsInitialize(struct DSDPDualMat_Ops* sops){ 00601 int info; 00602 if (sops==NULL) return 0; 00603 info=DSDPDualMatOpsInitialize(sops);DSDPCHKERR(info); 00604 sops->matseturmat=DDenseSetXMat; 00605 sops->matcholesky=DTPUMatCholeskyFactor; 00606 sops->matsolveforward=DTPUMatCholeskyForward; 00607 sops->matsolvebackward=DTPUMatCholeskyBackward; 00608 sops->matinvert=DTPUMatInvert; 00609 sops->matinverseadd=DTPUMatInverseAdd; 00610 sops->matinversemultiply=DTPUMatInverseMult; 00611 sops->matforwardmultiply=DenseSymPSDCholeskyForwardMultiply; 00612 sops->matfull=DTPUMatFull; 00613 sops->matdestroy=DTPUMatDestroy; 00614 sops->matgetsize=DTPUMatGetSize; 00615 sops->matview=DTPUMatView; 00616 sops->matlogdet=DTPUMatLogDet; 00617 sops->matname=lapackname; 00618 sops->id=1; 00619 return 0; 00620 } 00621 00622 00623 #undef __FUNCT__ 00624 #define __FUNCT__ "DSDPLAPACKDualMatCreate" 00625 int DSDPLAPACKPUDualMatCreate(int n,struct DSDPDualMat_Ops **sops, void**smat){ 00626 int info,nn=n*(n+1)/2; 00627 double *vv; 00628 dtpumat*AA; 00629 DSDPFunctionBegin; 00630 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info); 00631 info=DTPUMatCreateWData(n,vv,nn,&AA); DSDPCHKERR(info); 00632 AA->owndata=1;; 00633 AA->scaleit=1; 00634 info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info); 00635 *sops=&sdmatops; 00636 *smat=(void*)AA; 00637 DSDPFunctionReturn(0); 00638 } 00639 00640 static int switchptr(void* SD,void *SP){ 00641 dtpumat *s1,*s2; 00642 s1=(dtpumat*)(SD); 00643 s2=(dtpumat*)(SP); 00644 s1->v2=s2->val; 00645 s2->v2=s1->val; 00646 return 0; 00647 } 00648 00649 00650 #undef __FUNCT__ 00651 #define __FUNCT__ "DSDPLAPACKDualMatCreate2" 00652 int DSDPLAPACKPUDualMatCreate2(int n, 00653 struct DSDPDualMat_Ops **sops1, void**smat1, 00654 struct DSDPDualMat_Ops **sops2, void**smat2){ 00655 int info; 00656 DSDPFunctionBegin; 00657 info=DSDPLAPACKPUDualMatCreate(n,sops1,smat1);DSDPCHKERR(info); 00658 info=DSDPLAPACKPUDualMatCreate(n,sops2,smat2);DSDPCHKERR(info); 00659 info=switchptr(*smat1,*smat2); 00660 DSDPFunctionReturn(0); 00661 } 00662 00663 00664 typedef struct { 00665 int neigs; 00666 double *eigval; 00667 double *an; 00668 } Eigen; 00669 00670 typedef struct { 00671 dtpumat* AA; 00672 double alpha; 00673 Eigen Eig; 00674 } dvechmat; 00675 00676 #undef __FUNCT__ 00677 #define __FUNCT__ "CreateDvechmatWData" 00678 static int CreateDvechmatWdata(int n, double alpha, double vv[], dvechmat **A){ 00679 int info,nn=(n*n+n)/2; 00680 dvechmat* V; 00681 DSDPCALLOC1(&V,dvechmat,&info);DSDPCHKERR(info); 00682 info=DTPUMatCreateWData(n,vv,nn,&V->AA); DSDPCHKERR(info); 00683 V->Eig.neigs=-1; 00684 V->Eig.eigval=0; 00685 V->Eig.an=0; 00686 V->alpha=alpha; 00687 *A=V; 00688 return 0; 00689 } 00690 00691 00692 static int DvechmatGetRowNnz(void* AA, int trow, int nz[], int *nnzz,int n){ 00693 int k; 00694 *nnzz=n; 00695 for (k=0;k<n;k++) nz[k]++; 00696 return 0; 00697 } 00698 00699 static int DTPUMatGetRowAdd(void* AA, int nrow, double ytmp, double row[], int n){ 00700 dtpumat* A=(dtpumat*) AA; 00701 ffinteger i,k,nnn=n; 00702 double *v=A->val; 00703 00704 nnn=nrow*(nrow+1)/2; 00705 for (i=0;i<nrow;i++){ 00706 row[i]+=ytmp*v[nnn+i]; 00707 } 00708 row[nrow]+=ytmp*v[nnn+nrow]; 00709 for (i=nrow+1;i<n;i++){ 00710 k=i*(i+1)/2+nrow; 00711 row[i]+=ytmp*v[k]; 00712 } 00713 return 0; 00714 } 00715 00716 static int DvechmatGetRowAdd(void* AA, int trow, double scl, double r[], int m){ 00717 int info; 00718 dvechmat* A=(dvechmat*)AA; 00719 info=DTPUMatGetRowAdd((void*)A->AA ,trow,scl*A->alpha,r,m); 00720 return 0; 00721 } 00722 00723 static int DvechmatAddMultiple(void* AA, double alpha, double r[], int nnn, int n){ 00724 dvechmat* A=(dvechmat*)AA; 00725 ffinteger nn=nnn, ione=1; 00726 double *val=A->AA->val; 00727 alpha*=A->alpha; 00728 daxpy(&nn,&alpha,val,&ione,r,&ione); 00729 return 0; 00730 } 00731 00732 00733 static int DvechEigVecVec(void*, double[], int, double*); 00734 static int DvechmatVecVec(void* AA, double x[], int n, double *v){ 00735 dvechmat* A=(dvechmat*)AA; 00736 int i,j,k=0; 00737 double dd=0,*val=A->AA->val; 00738 *v=0.0; 00739 if (A->Eig.neigs<n/5){ 00740 i=DvechEigVecVec(AA,x,n,&dd); 00741 *v=dd*A->alpha; 00742 } else { 00743 for (i=0; i<n; i++){ 00744 for (j=0;j<i;j++){ 00745 dd+=2*x[i]*x[j]*val[k]; 00746 k++; 00747 } 00748 dd+=x[i]*x[i]*val[k]; 00749 k++; 00750 } 00751 *v=dd*A->alpha; 00752 } 00753 return 0; 00754 } 00755 00756 00757 static int DvechmatFNorm2(void* AA, int n, double *v){ 00758 dvechmat* A=(dvechmat*)AA; 00759 long int i,j,k=0; 00760 double dd=0,*x=A->AA->val; 00761 for (i=0; i<n; i++){ 00762 for (j=0;j<i;j++){ 00763 dd+=2*x[k]*x[k]; 00764 k++; 00765 } 00766 dd+=x[k]*x[k]; 00767 k++; 00768 } 00769 *v=dd*A->alpha*A->alpha; 00770 return 0; 00771 } 00772 00773 00774 static int DvechmatCountNonzeros(void* AA, int *nnz, int n){ 00775 *nnz=n*(n+1)/2; 00776 return 0; 00777 } 00778 00779 00780 static int DvechmatDot(void* AA, double x[], int nn, int n, double *v){ 00781 dvechmat* A=(dvechmat*)AA; 00782 ffinteger ione=1,nnn=nn; 00783 double dd,*val=A->AA->val; 00784 dd=ddot(&nnn,val,&ione,x,&ione); 00785 *v=2*dd*A->alpha; 00786 return 0; 00787 } 00788 00789 /* 00790 static int DvechmatNormF2(void* AA, int n, double *v){ 00791 dvechmat* A=(dvechmat*)AA; 00792 return(DTPUMatNormF2((void*)(A->AA), n,v)); 00793 } 00794 */ 00795 #undef __FUNCT__ 00796 #define __FUNCT__ "DvechmatDestroy" 00797 static int DvechmatDestroy(void* AA){ 00798 dvechmat* A=(dvechmat*)AA; 00799 int info; 00800 info=DTPUMatDestroy((void*)(A->AA)); 00801 if (A->Eig.an){DSDPFREE(&A->Eig.an,&info);DSDPCHKERR(info);} 00802 if (A->Eig.eigval){DSDPFREE(&A->Eig.eigval,&info);DSDPCHKERR(info);} 00803 A->Eig.neigs=-1; 00804 DSDPFREE(&A,&info);DSDPCHKERR(info); 00805 return 0; 00806 } 00807 00808 00809 static int DvechmatView(void* AA){ 00810 dvechmat* A=(dvechmat*)AA; 00811 dtpumat* M=A->AA; 00812 int i,j,kk=0; 00813 double *val=M->val; 00814 for (i=0; i<M->n; i++){ 00815 for (j=0; j<=i; j++){ 00816 printf(" %4.2e",A->alpha*val[kk]); 00817 kk++; 00818 } 00819 printf(" \n"); 00820 } 00821 return 0; 00822 } 00823 00824 00825 #undef __FUNCT__ 00826 #define __FUNCT__ "DSDPCreateDvechmatEigs" 00827 static int CreateEigenLocker(Eigen *E,int neigs, int n){ 00828 int info; 00829 DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info); 00830 DSDPCALLOC2(&E->an,double,n*neigs,&info);DSDPCHKERR(info); 00831 E->neigs=neigs; 00832 return 0; 00833 } 00834 00835 00836 static int EigMatSetEig(Eigen* A,int row, double eigv, double v[], int n){ 00837 double *an=A->an; 00838 A->eigval[row]=eigv; 00839 memcpy((void*)(an+n*row),(void*)v,n*sizeof(double)); 00840 return 0; 00841 } 00842 00843 00844 static int EigMatGetEig(Eigen* A,int row, double *eigenvalue, double eigenvector[], int n){ 00845 double* an=A->an; 00846 *eigenvalue=A->eigval[row]; 00847 memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double)); 00848 return 0; 00849 } 00850 00851 00852 static int DvechmatComputeEigs(dvechmat*,double[],int,double[],int,double[],int,int[],int); 00853 00854 static int DvechmatFactor(void*AA, double dmatp[], int nn0, double dwork[], int n, double ddwork[], int n1, int iptr[], int n2){ 00855 00856 int info; 00857 dvechmat* A=(dvechmat*)AA; 00858 if (A->Eig.neigs>=0) return 0; 00859 info=DvechmatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2);DSDPCHKERR(info); 00860 return 0; 00861 } 00862 00863 static int DvechmatGetRank(void *AA,int *rank, int n){ 00864 dvechmat* A=(dvechmat*)AA; 00865 if (A->Eig.neigs>=0){ 00866 *rank=A->Eig.neigs; 00867 } else { 00868 DSDPSETERR(1,"Vech Matrix not factored yet\n"); 00869 } 00870 return 0; 00871 } 00872 00873 static int DvechmatGetEig(void* AA, int rank, double *eigenvalue, double vv[], int n, int indz[], int *nind){ 00874 dvechmat* A=(dvechmat*)AA; 00875 int i,info; 00876 double dd; 00877 if (A->Eig.neigs>0){ 00878 info=EigMatGetEig(&A->Eig,rank,&dd,vv,n);DSDPCHKERR(info); 00879 *nind=n; 00880 *eigenvalue=dd*A->alpha; 00881 for (i=0;i<n;i++){ indz[i]=i;} 00882 } else { 00883 DSDPSETERR(1,"Vech Matrix not factored yet\n"); 00884 } 00885 return 0; 00886 } 00887 00888 static int DvechEigVecVec(void* AA, double v[], int n, double *vv){ 00889 dvechmat* A=(dvechmat*)AA; 00890 int i,rank,neigs; 00891 double *an,dd,ddd=0,*eigval; 00892 if (A->Eig.neigs>=0){ 00893 an=A->Eig.an; 00894 neigs=A->Eig.neigs; 00895 eigval=A->Eig.eigval; 00896 for (rank=0;rank<neigs;rank++){ 00897 for (dd=0,i=0;i<n;i++){ 00898 dd+=v[i]*an[i]; 00899 } 00900 an+=n; 00901 ddd+=dd*dd*eigval[rank]; 00902 } 00903 *vv=ddd*A->alpha; 00904 } else { 00905 DSDPSETERR(1,"Vech Matrix not factored yet\n"); 00906 } 00907 return 0; 00908 } 00909 00910 00911 static struct DSDPDataMat_Ops dvechmatops; 00912 static const char *datamatname="DENSE VECH MATRIX"; 00913 00914 static int DvechmatOpsInitialize(struct DSDPDataMat_Ops *sops){ 00915 int info; 00916 if (sops==NULL) return 0; 00917 info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info); 00918 sops->matvecvec=DvechmatVecVec; 00919 sops->matdot=DvechmatDot; 00920 sops->mataddrowmultiple=DvechmatGetRowAdd; 00921 sops->mataddallmultiple=DvechmatAddMultiple; 00922 sops->matview=DvechmatView; 00923 sops->matdestroy=DvechmatDestroy; 00924 sops->matfactor2=DvechmatFactor; 00925 sops->matgetrank=DvechmatGetRank; 00926 sops->matgeteig=DvechmatGetEig; 00927 sops->matrownz=DvechmatGetRowNnz; 00928 sops->matfnorm2=DvechmatFNorm2; 00929 sops->matnnz=DvechmatCountNonzeros; 00930 sops->id=1; 00931 sops->matname=datamatname; 00932 return 0; 00933 } 00934 00935 #undef __FUNCT__ 00936 #define __FUNCT__ "DSDPGetDmat" 00937 int DSDPGetDMat(int n,double alpha, double *val, struct DSDPDataMat_Ops**sops, void**smat){ 00938 int info,k; 00939 double dtmp; 00940 dvechmat* A; 00941 DSDPFunctionBegin; 00942 00943 for (k=0;k<n*(n+1)/2;++k) dtmp=val[k]; 00944 info=CreateDvechmatWdata(n,alpha,val,&A); DSDPCHKERR(info); 00945 A->Eig.neigs=-1; 00946 info=DvechmatOpsInitialize(&dvechmatops); DSDPCHKERR(info); 00947 if (sops){*sops=&dvechmatops;} 00948 if (smat){*smat=(void*)A;} 00949 DSDPFunctionReturn(0); 00950 } 00951 00952 00953 #undef __FUNCT__ 00954 #define __FUNCT__ "DvechmatComputeEigs" 00955 static int DvechmatComputeEigs(dvechmat* AA,double DD[], int nn0, double W[], int n, double WORK[], int n1, int iiptr[], int n2){ 00956 00957 int i,j,k,neigs,info; 00958 long int *i2darray=(long int*)DD; 00959 int ownarray1=0,ownarray2=0,ownarray3=0; 00960 double *val=AA->AA->val; 00961 double *dmatarray=0,*dworkarray=0,eps=1.0e-12; 00962 int nn1=0,nn2=0; 00963 00964 /* create a dense array in which to put numbers */ 00965 if (n*n>nn1){ 00966 DSDPCALLOC2(&dmatarray,double,(n*n),&info); DSDPCHKERR(info); 00967 ownarray1=1; 00968 } 00969 memset((void*)dmatarray,0,n*n*sizeof(double)); 00970 00971 if (n*n>nn2){ 00972 DSDPCALLOC2(&dworkarray,double,(n*n),&info); DSDPCHKERR(info); 00973 ownarray2=1; 00974 } 00975 00976 if (n*n*sizeof(long int)>nn0*sizeof(double)){ 00977 DSDPCALLOC2(&i2darray,long int,(n*n),&info); DSDPCHKERR(info); 00978 ownarray3=1; 00979 } 00980 00981 00982 for (k=0,i=0; i<n; i++){ 00983 for (j=0; j<=i; j++){ 00984 dmatarray[i*n+j] += val[k]; 00985 if (i!=j){ 00986 dmatarray[j*n+i] += val[k]; 00987 } 00988 k++; 00989 } 00990 } 00991 /* Call LAPACK to compute the eigenvalues */ 00992 info=DSDPGetEigs(dmatarray,n,dworkarray,n*n,i2darray,n*n, 00993 W,n,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info); 00994 00995 /* Count the nonzero eigenvalues */ 00996 for (neigs=0,i=0;i<n;i++){ 00997 if (fabs(W[i])> eps ){ neigs++;} 00998 } 00999 01000 info=CreateEigenLocker(&AA->Eig,neigs,n);DSDPCHKERR(info); 01001 01002 /* Copy into structure */ 01003 for (neigs=0,i=0; i<n; i++){ 01004 if (fabs(W[i]) > eps){ 01005 info=EigMatSetEig(&AA->Eig,neigs,W[i],dmatarray+n*i,n);DSDPCHKERR(info); 01006 neigs++; 01007 } 01008 } 01009 01010 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);} 01011 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);} 01012 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);} 01013 return 0; 01014 } 01015