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
00038
00039
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
00434
00435
00436
00437
00438
00439
00440
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
00791
00792
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
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
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
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
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