cholmat2.c

Go to the documentation of this file.
00001 #include "numchol.h"
00002 #include "dsdpdualmat_impl.h"
00003 #include "dsdpsys.h"
00004 #include "dsdplapack.h"
00005 
00011 typedef struct{
00012   chfac* spsym;
00013   double *sinv;
00014   char   UPLQ;
00015   int n;
00016   int dsinv;
00017 } spmat;
00018 
00019 
00020 static int SMatDestroy(void*S){
00021   spmat* SS=(spmat*)S;
00022   int info;
00023   CfcFree(&SS->spsym);
00024   if (SS->dsinv){
00025     DSDPFREE(&SS->sinv,&info);
00026   }
00027   DSDPFREE(&SS,&info);
00028   return 0;
00029 }
00030 
00031 static int SMatGetSize(void *S, int *n){
00032   spmat* SS=(spmat*)S;
00033   *n=SS->spsym->nrow;
00034   return 0;
00035 }
00036 
00037 static int SMatView(void* S){
00038   spmat* SS=(spmat*)S;
00039   int info;
00040   info=Mat4View(SS->spsym); DSDPCHKERR(info);
00041   return 0;
00042 }
00043 
00044 static int SMatLogDet(void* S, double *dd){
00045   spmat* SS=(spmat*)S;
00046   int info;
00047   info=Mat4LogDet(SS->spsym,dd); DSDPCHKERR(info);
00048   return 0;
00049 }
00050 
00051 
00052 
00053 static int SMatSetURMatP(void*S, double v[], int nn, int n){
00054   spmat* SS=(spmat*)S;
00055   int k,j,row,info;
00056   double *rw1,*rw2,*xr=v;
00057   rw1=SS->spsym->rw;rw2=rw1+n;
00058   info=MatZeroEntries4(SS->spsym); DSDPCHKERR(info);
00059   for (k=0;k<n/2;k++){
00060     row = 2*k;
00061 
00062     xr=v+row*(row+1)/2;
00063     memcpy((void*)rw1,(void*)(xr),(row+1)*sizeof(double));
00064     xr+=row+1;
00065     rw1[row+1]=xr[row];
00066     memcpy((void*)(rw2),(void*)(xr),(row+2)*sizeof(double));
00067     xr+=row+2;
00068 
00069     /*  memset((void*)rw,0,(n-row)*sizeof(double)); */
00070     for (j=row+2;j<n;j++){ 
00071       rw1[j]=xr[row];
00072       rw2[j]=xr[row+1];
00073       xr+=j+1; 
00074     }
00075 
00076     info=MatSetColumn4(SS->spsym,rw1,row);  DSDPCHKERR(info);
00077     info=MatSetColumn4(SS->spsym,rw2,row+1);  DSDPCHKERR(info);
00078   }
00079 
00080   for (row=2*(n/2);row<n;row++){
00081     /*  memset((void*)rw,0,(n-row)*sizeof(double)); */
00082     xr=v+row*(row+1)/2;
00083     memcpy((void*)(rw1),(void*)(xr),(row+1)*sizeof(double));
00084     xr+=row+1;
00085     for (j=row+1;j<n;j++){ rw1[j]=xr[row]; xr+=(j+2); }
00086     info=MatSetColumn4(SS->spsym,rw1,row);  DSDPCHKERR(info);
00087     xr+=(n-row);
00088   }
00089   return 0;
00090 }
00091 
00092 static int SMatSetURMatU(void*S, double v[], int nn, int n){
00093   spmat* SS=(spmat*)S;
00094   int k,j,row,info;
00095   double *rw1,*rw2,*xr=v;
00096   rw1=SS->spsym->rw;rw2=rw1+n;
00097   info=MatZeroEntries4(SS->spsym); DSDPCHKERR(info);
00098   for (k=0;k<n/2;k++){
00099     row = 2*k;
00100 
00101     xr=v+row*n;
00102     memcpy((void*)rw1,(void*)(xr),(row+1)*sizeof(double));
00103     xr+=n;
00104     rw1[row+1]=xr[row];
00105     memcpy((void*)(rw2),(void*)(xr),(row+2)*sizeof(double));
00106     xr+=n;
00107     /*  memset((void*)rw,0,(n-row)*sizeof(double)); */
00108     for (j=row+2;j<n;j++){ 
00109       rw1[j]=xr[row];
00110       rw2[j]=xr[row+1];
00111       xr+=n; 
00112     }
00113 
00114     info=MatSetColumn4(SS->spsym,rw1,row);  DSDPCHKERR(info);
00115     info=MatSetColumn4(SS->spsym,rw2,row+1);  DSDPCHKERR(info);
00116   }
00117 
00118   for (row=2*(n/2);row<n;row++){
00119     /*  memset((void*)rw,0,(n-row)*sizeof(double)); */
00120     xr=v+row*n;
00121     memcpy((void*)(rw1),(void*)(xr),(row+1)*sizeof(double));
00122     xr+=n;
00123     for (j=row+1;j<n;j++){ rw1[j]=xr[row]; xr+=n; }
00124     info=MatSetColumn4(SS->spsym,rw1,row);  DSDPCHKERR(info);
00125   }
00126   return 0;
00127 }
00128 
00129 static int SMatSetURMat(void*S, double v[], int nn, int n){
00130   spmat* SS=(spmat*)S;
00131   int info;
00132   if (SS->UPLQ=='P'){
00133     info=SMatSetURMatP(S,v,nn,n);DSDPCHKERR(info);
00134   } else if (SS->UPLQ=='U'){
00135     info=SMatSetURMatU(S,v,nn,n);DSDPCHKERR(info);
00136   }
00137   return 0;
00138 }
00139 
00140 static int SMatSolve(void *S, int indx[], int nind, double b[], double x[],int n){
00141   spmat* SS=(spmat*)S;
00142   int i,ii;
00143   double alpha,*s1=SS->sinv,*s2;
00144   ffinteger nn,ione;
00145   if (SS->sinv && nind < n/4){
00146     memset((void*)x,0,n*sizeof(double));
00147     for (i=0;i<nind;i++){
00148       ii=indx[i];
00149       ione=1;nn=n;alpha=b[ii];s2=s1+n*ii;
00150       daxpy(&nn,&alpha,s2,&ione,x,&ione);
00151     }
00152   } else {
00153     memcpy(x,b,n*sizeof(double));
00154     ChlSolve(SS->spsym, b, x);
00155   }
00156   return 0;
00157 }
00158 
00159 static int SMatCholeskySolveBackward(void *S, double b[], double x[],int n){
00160   spmat* SS=(spmat*)S;
00161   ChlSolveBackward(SS->spsym, b, x);
00162   return 0;
00163 }
00164 
00165 static int SMatCholeskySolveForward(void *S, double b[], double x[], int n){
00166   spmat* SS=(spmat*)S;
00167   ChlSolveForward(SS->spsym, b, x);
00168   return 0;
00169 }
00170 
00171 static int SMatFull(void *S, int *full){
00172   *full=0;
00173   return 0;
00174 }
00175 
00176 static int SMatCholeskyFactor(void *S, int *flag){
00177   spmat* SS=(spmat*)S;
00178   int *iw;
00179   double *rw;
00180   cfc_sta Cfact;
00181   iw=SS->spsym->iw;
00182   rw=SS->spsym->rw;
00183   Cfact=(cfc_sta)ChlFact(SS->spsym,iw,rw,TRUE);
00184   if (CfcOk!= Cfact){ 
00185     *flag=1;
00186   } else {
00187     *flag=0;
00188   }
00189   return 0;
00190 }
00191 
00192 static int SMatInverseAddP(void *S, double alpha, double v[], int nn, int n){
00193   spmat* SS=(spmat*)S;
00194   ffinteger ii,ione=1;
00195   double *x,*b,*ss=SS->sinv;
00196   int i,j,k=0;
00197 
00198   if (ss){
00199     for (i=0;i<n;i++){
00200       v+=i; ii=i+1;
00201       daxpy(&ii,&alpha,ss,&ione,v,&ione);
00202       ss+=n;
00203     }
00204   } else {
00205     b=SS->spsym->rw;x=b+n;
00206     for (i=0;i<n;i++){
00207       memset((void*)b,0,n*sizeof(double));
00208       b[i]=alpha;
00209       ChlSolve(SS->spsym, b, x);
00210       k=k+i;
00211       for (j=0;j<=i;j++){
00212         v[k+j]+=x[j];
00213       }
00214     }
00215   }
00216   return 0;
00217 }
00218 
00219 static int SMatInverseAddU(void *S, double alpha, double v[], int nn, int n){
00220   spmat* SS=(spmat*)S;
00221   ffinteger n2=n*n,ione=1;
00222   double *x,*b,*ss=SS->sinv;
00223   int i,j,k=0;
00224 
00225   if (ss){
00226     daxpy(&n2,&alpha,ss,&ione,v,&ione);
00227   } else {
00228     b=SS->spsym->rw;x=b+n;
00229     for (i=0;i<n;i++){
00230       memset((void*)b,0,n*sizeof(double));
00231       b[i]=alpha;
00232       ChlSolve(SS->spsym, b, x);
00233       k=i*n;
00234       for (j=0;j<n;j++){
00235         v[k+j]+=x[j];
00236       }
00237     }
00238   }
00239   return 0;
00240 }
00241 
00242 static int SMatInverseAdd(void *S, double alpha, double v[], int nn, int n){
00243   spmat* SS=(spmat*)S;
00244   int info;
00245   if (SS->UPLQ=='P'){
00246     info=SMatInverseAddP(S,alpha,v,nn,n);DSDPCHKERR(info);
00247   } else if (SS->UPLQ=='U'){
00248     info=SMatInverseAddU(S,alpha,v,nn,n);DSDPCHKERR(info);
00249   }
00250   return 0;
00251 }
00252 
00253 static int SMatCholeskyForwardMultiply(void *S, double b[], double x[], int n){
00254   spmat* SS=(spmat*)S;
00255   GetUhat(SS->spsym,b,x);
00256   return 0;
00257 }
00258 
00259 static int SMatInvert(void *S){
00260   spmat* SS=(spmat*)S;
00261   double *x,*b,*v=SS->sinv;
00262   int i,n=SS->n;
00263 
00264   b=SS->spsym->rw;x=b+n;
00265 
00266   if (v){
00267     for (i=0;i<n;i++){
00268       memset((void*)b,0,n*sizeof(double));
00269       b[i]=1.0;
00270       ChlSolve(SS->spsym, b, x);
00271       memcpy((void*)(v+i*n),(void*)x,n*sizeof(double));
00272     }
00273   }
00274   return 0;
00275 }
00276 
00277 static struct  DSDPDualMat_Ops sdmatops;
00278 static const char* tmatname="SPARSE PSD";
00279 static int SDualOpsInitialize(struct  DSDPDualMat_Ops* sops){
00280   int info;
00281   if (sops==NULL) return 0;
00282   info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
00283   sops->matcholesky=SMatCholeskyFactor;
00284   sops->matsolveforward=SMatCholeskySolveForward;
00285   sops->matsolvebackward=SMatCholeskySolveBackward;
00286   sops->matinversemultiply=SMatSolve;
00287   sops->matinvert=SMatInvert;
00288   sops->matinverseadd=SMatInverseAdd;
00289   sops->matforwardmultiply=SMatCholeskyForwardMultiply;
00290   sops->matseturmat=SMatSetURMat;
00291   sops->matfull=SMatFull;
00292   sops->matdestroy=SMatDestroy;
00293   sops->matgetsize=SMatGetSize;
00294   sops->matview=SMatView;
00295   sops->matlogdet=SMatLogDet;
00296   sops->matname=tmatname;
00297   return 0;
00298  }
00299 
00300 static int dcholmatcreate(int n, char UPLQ, chfac *sp, 
00301                            struct  DSDPDualMat_Ops **sops, void**smat){
00302   spmat  *S;
00303   int info;
00304   DSDPCALLOC1(&S,spmat,&info);DSDPCHKERR(info);
00305   S->UPLQ=UPLQ; S->n=n; S->sinv=0; S->dsinv=0; S->spsym=sp;
00306   info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
00307   *sops=&sdmatops;
00308   *smat=(void*)S;
00309   return 0;
00310  }
00311 
00312 static int dcholmatsinverse(int n, spmat *S1, spmat *S2){
00313   int info;
00314   double *ssinv;
00315   DSDPCALLOC2(&ssinv,double,n*n,&info);
00316   S1->sinv=ssinv; S2->sinv=ssinv; S2->dsinv=1;
00317   return 0;
00318 }
00319 
00320 #undef __FUNCT__
00321 #define __FUNCT__ "DSDPDenseDualMatCreate"
00322 int DSDPDenseDualMatCreate(int n, char UPLQ, 
00323                            struct  DSDPDualMat_Ops **sops1, void**smat1,
00324                            struct  DSDPDualMat_Ops **sops2, void**smat2){
00325   int info=0;
00326   chfac *sp;
00327 
00328   DSDPFunctionBegin;
00329   info=MchlSetup2(n,&sp); DSDPCHKERR(info);
00330   info=dcholmatcreate(n,UPLQ,sp,sops1,smat1);DSDPCHKERR(info);
00331   info=MchlSetup2(n,&sp); DSDPCHKERR(info);
00332   info=dcholmatcreate(n,UPLQ,sp,sops1,smat2);DSDPCHKERR(info);
00333   info=dcholmatsinverse(n,(spmat*)(*smat1),(spmat*)(*smat2));DSDPCHKERR(info);
00334 
00335   DSDPFunctionReturn(0);
00336 }
00337 
00338 
00339 #undef __FUNCT__
00340 #define __FUNCT__ "DSDPSparseDualMatCreate"
00341 int DSDPSparseDualMatCreate(int n, int *rnnz, int *snnz, 
00342                             int trank,char UPLQ,int*nnzz, 
00343                             struct  DSDPDualMat_Ops **sops1, void**smat1,
00344                             struct  DSDPDualMat_Ops **sops2, void**smat2){
00345   int nnz,info=0;
00346   chfac *sp;
00347 
00348   DSDPFunctionBegin;
00349   SymbProc(rnnz,snnz,n,&sp); DSDPCHKERR(info);
00350   info=dcholmatcreate(n,UPLQ,sp,sops1,smat1);DSDPCHKERR(info);
00351   SymbProc(rnnz,snnz,n,&sp); DSDPCHKERR(info);
00352   info=dcholmatcreate(n,UPLQ,sp,sops2,smat2);DSDPCHKERR(info);
00353   nnz=sp->unnz;*nnzz=nnz;
00354   if (trank>2*n+2){
00355     info=dcholmatsinverse(n,(spmat*)(*smat1),(spmat*)(*smat2));DSDPCHKERR(info);
00356   }
00357   DSDPFunctionReturn(0);
00358 }

Generated on Thu May 22 09:42:47 2008 for DSDP by  doxygen 1.5.5