Main Page | Modules | Alphabetical List | Data Structures | Directories | File List | Globals | Related Pages

dlpack.c

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 

Generated on Sat Oct 15 11:05:41 2005 for DSDP by  doxygen 1.4.2