cholmat.c

Go to the documentation of this file.
00001 #include "numchol.h"
00002 #include "dsdpschurmat_impl.h"
00003 #include "dsdpsys.h"
00004 #include "dsdpvec.h"
00005 #include "dsdp5.h"
00006 
00011 int DSDPSparsityInSchurMat(DSDP, int, int *, int);
00012 int DSDPSetSchurMatOps(DSDP,struct DSDPSchurMat_Ops*,void*);
00013 
00014 typedef struct {
00015   chfac *M;
00016   int m;
00017   int isdense;
00018   int *rnnz;
00019   int *colnnz;
00020   int nnz;
00021   DSDPVec    D1;
00022   DSDP     dsdp;
00023 } MCholSolverALL;
00024 
00025 
00026 static int dsdpuselapack=1;
00027 #undef __FUNCT__  
00028 #define __FUNCT__ "DSDPUseLAPACKForSchur"
00029 int DSDPUseLAPACKForSchur(DSDP dsdp,int flag){ 
00030   DSDPFunctionBegin;
00031   dsdpuselapack = flag;
00032   DSDPFunctionReturn(0);
00033 }
00034 
00035 #undef __FUNCT__
00036 #define __FUNCT__ "Taddline"
00037 static int Taddline(void* M, int row, double dd, double v[], int m){
00038   int info;
00039   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00040   DSDPFunctionBegin;
00041   info=MatAddColumn4(AMA->M,dd,v,row);  DSDPCHKERR(info);
00042   DSDPFunctionReturn(0);
00043 }
00044 
00045 #undef __FUNCT__
00046 #define __FUNCT__ "Tadddiagonal"
00047 static int Tadddiagonal(void* M, int row, double v){
00048   int info;
00049   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00050   DSDPFunctionBegin;
00051   info=MatAddDiagonalElement(AMA->M,row,v);  DSDPCHKERR(info);
00052   DSDPFunctionReturn(0);
00053 }
00054 
00055 #undef __FUNCT__
00056 #define __FUNCT__ "Tassemble"
00057 static int Tassemble(void*M){
00058   DSDPFunctionBegin;
00059   DSDPFunctionReturn(0);
00060 }
00061 
00062 
00063 #undef __FUNCT__
00064 #define __FUNCT__ "DSDPCheckForSparsity"
00065 static int DSDPCheckForSparsity( DSDP dsdp, int m, int *rnnz, int tnnz[], int *totalnnz){
00066   int info,i,j,tottalnnz=0;
00067   DSDPFunctionBegin;
00068   memset(rnnz,0,(m+1)*sizeof(int));
00069   for (i=0;i<m;i++){
00070     info=DSDPSparsityInSchurMat(dsdp,i,tnnz,m); DSDPCHKERR(info);
00071     for (j=i+1;j<m;j++){ if (tnnz[j]>0) rnnz[i+1]++;} 
00072   }
00073 
00074   for (i=0;i<m;i++){tottalnnz+=rnnz[i+1];}
00075   *totalnnz=tottalnnz;
00076   DSDPFunctionReturn(0);
00077 }
00078 
00079 #undef __FUNCT__
00080 #define __FUNCT__ "DSDPCreateM"
00081 static int DSDPCreateM(MCholSolverALL *ABA, chfac **M, int rrnnz[],int tnnz[], int totalnnz){
00082 
00083   int *snnz,*rnnz;
00084   int i,j,tt,info;
00085   int col,*perm,k;
00086   chfac * sfptr;
00087   int n=ABA->m,m=ABA->m;
00088   DSDP dsdp=ABA->dsdp;
00089  
00090   DSDPFunctionBegin;
00091 
00092   DSDPCALLOC2(&snnz,int,(totalnnz+1),&info); DSDPCHKERR(info);
00093   DSDPCALLOC2(&rnnz,int,(m+1),&info); DSDPCHKERR(info);
00094   for (i=0;i<=m;i++){ rnnz[i]=rrnnz[i];}
00095   tt=0;
00096   for (i=0;i<m;i++){
00097     info=DSDPSparsityInSchurMat(dsdp,i,tnnz,m); DSDPCHKERR(info);
00098     for (j=i+1;j<m;j++){ if (tnnz[j]>0) {snnz[tt]=j; tt++;} } 
00099   }
00100 
00101   printf("Trying Sparse M: Total nonzeros: %d of %d \n",totalnnz,m*(m-1)/2 );
00102   /* Create sparse dual matrix structure */
00103   SymbProc(rnnz+1,snnz,n,&sfptr);
00104   ABA->isdense=0;
00105   ABA->M=sfptr;
00106   ABA->nnz=totalnnz;
00107   ABA->rnnz=rnnz;
00108   ABA->colnnz=snnz;
00109   *M=sfptr;
00110 
00111   for (i=0;i<m;i++){
00112     rnnz[i+1]+=rnnz[i];
00113   }
00114 
00115   perm=sfptr->invp;
00116   for (i=m-1;i>=0;i--){
00117     for (j=rnnz[i+1]-1;j>=rnnz[i];j--){
00118       col=snnz[j];
00119       tt=perm[col];
00120       if (perm[i] > tt ){
00121         for (k=j;k<rnnz[col]-1;k++){ snnz[k]=snnz[k+1];}
00122         for (k=i;k<col;k++){ rnnz[k+1]--;}
00123         snnz[rnnz[col]]=i;
00124       }
00125     }
00126   }
00127   
00128   DSDPFunctionReturn(0);
00129 }
00130 
00131 
00132 #undef __FUNCT__  
00133 #define __FUNCT__ "DSDPLinearSolverPrepare"
00134 static int DSDPLinearSolverPrepare(void* ctx,int*flag){
00135 
00136   cfc_sta Cfact;
00137   chfac *sf;
00138   MCholSolverALL *AMA = (MCholSolverALL *)ctx; 
00139 
00140   DSDPFunctionBegin;
00141   *flag=0;
00142   sf=AMA->M;  
00143   /*
00144   Cfact=(cfc_sta)ChlFact(sf,sf->iw,sf->rw,FALSE);
00145   */
00146   Cfact=(cfc_sta)ChlFact(sf,sf->iw,sf->rw,TRUE);
00147   if (CfcOk!=Cfact ){ *flag=1; /*  printf("Not Good factorization \n"); */ }
00148   DSDPFunctionReturn(0);
00149 }
00150 
00151 
00152 #undef __FUNCT__  
00153 #define __FUNCT__ "DSDPLinearSolve"
00154 static int DSDPLinearSolve2(void* ctx, double dd[], double xx[], int n){
00155 
00156   int i,info;
00157   double *bb;
00158   MCholSolverALL *AMA = (MCholSolverALL *)ctx; 
00159 
00160   DSDPFunctionBegin;
00161   info=DSDPVecGetArray(AMA->D1, &bb);DSDPCHKERR(info);
00162   for (i=0;i<n;i++){ bb[i]=dd[i];}
00163   ChlSolve(AMA->M, bb, xx);
00164   info=DSDPVecRestoreArray(AMA->D1, &bb);
00165   DSDPFunctionReturn(0);
00166 }
00167 
00168 
00169 #undef __FUNCT__
00170 #define __FUNCT__ "DSDPGramMatRowNonzeros"
00171 static int DSDPGramMatRowNonzeros(void*M, int row, double cols[], int *ncols, int nrows){
00172   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00173   int i;
00174   DSDPFunctionBegin;
00175   if (AMA->isdense){
00176     *ncols = nrows-row;
00177     for (i=row;i<nrows;i++){
00178       cols[i]=1.0;
00179     }
00180   } else {
00181     *ncols = AMA->rnnz[row+1] - AMA->rnnz[row]+1;
00182     cols[row]=1.0;
00183     for (i=AMA->rnnz[row]; i<AMA->rnnz[row+1]; i++){
00184       cols[AMA->colnnz[i]]=1.0;
00185     }
00186   }
00187   DSDPFunctionReturn(0);
00188 }
00189 
00190 #undef __FUNCT__
00191 #define __FUNCT__ "Tzero"
00192 static int Tzero(void*M){
00193   int info;
00194   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00195   DSDPFunctionBegin;
00196   info=MatZeroEntries4(AMA->M); DSDPCHKERR(info);
00197   DSDPFunctionReturn(0);
00198 }
00199 
00200 #undef __FUNCT__
00201 #define __FUNCT__ "Tdestroy"
00202 static int Tdestroy(void*M){
00203   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00204   int info;
00205   DSDPFunctionBegin;
00206   CfcFree(&AMA->M);
00207   info = DSDPVecDestroy(&AMA->D1); DSDPCHKERR(info);
00208   if (AMA->isdense==0 && AMA->rnnz ){
00209     DSDPFREE(&AMA->rnnz,&info);DSDPCHKERR(info);
00210     DSDPFREE(&AMA->colnnz,&info);DSDPCHKERR(info);
00211   }
00212   DSDPFREE(&AMA,&info);DSDPCHKERR(info);
00213   DSDPFunctionReturn(0);
00214 }
00215 
00216 static int Tsetup(void*M, int m){
00217   DSDPFunctionBegin;
00218   DSDPFunctionReturn(0);
00219 }
00220 
00221 static int TTTMatMult(void*M,double x[],double y[],int n){
00222   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00223   int info;
00224   DSDPFunctionBegin;
00225   info=MatMult4(AMA->M,x,y,n); DSDPCHKERR(info);
00226   DSDPFunctionReturn(0);
00227 }
00228 
00229 static int TTTMatShiftDiagonal(void *M, double d){
00230   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00231   int info;
00232   DSDPFunctionBegin;
00233   info=Mat4DiagonalShift(AMA->M,d); DSDPCHKERR(info);
00234   DSDPFunctionReturn(0);
00235 }
00236 
00237 static int TTTMatAddDiagonal(void *M, double diag[], int m){
00238   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00239   int info;
00240   DSDPFunctionBegin;
00241   info=Mat4AddDiagonal(AMA->M,diag,m); DSDPCHKERR(info);
00242   DSDPFunctionReturn(0);
00243 }
00244 
00245 static int TTTMatView(void *M){
00246   MCholSolverALL *AMA = (MCholSolverALL *)M; 
00247   int info;
00248   DSDPFunctionBegin;
00249   info=Mat4View(AMA->M); DSDPCHKERR(info);
00250   DSDPFunctionReturn(0);
00251 }
00252 
00253 
00254 /*
00255 static int TTTMatGetDiagonal(void *M, double d[], int n){
00256   chfac*A = (chfac*)M;
00257   int info;
00258   DSDPFunctionBegin;
00259   info=Mat4GetDiagonal(A,d,n); DSDPCHKERR(info);
00260   DSDPFunctionReturn(0);
00261 }
00262 */
00263 static const char* tmatname="SPARSE PSD";
00264 
00265 static int TMatOpsInit(struct  DSDPSchurMat_Ops *mops){  
00266   int info;
00267   DSDPFunctionBegin;
00268   info=DSDPSchurMatOpsInitialize(mops); DSDPCHKERR(info);
00269   mops->matrownonzeros=DSDPGramMatRowNonzeros;
00270   mops->mataddrow=Taddline;
00271   mops->matadddiagonal=TTTMatAddDiagonal;
00272   mops->mataddelement=Tadddiagonal;
00273   mops->matshiftdiagonal=TTTMatShiftDiagonal;
00274   mops->matassemble=Tassemble;
00275   mops->matscaledmultiply=TTTMatMult;
00276   mops->matfactor=DSDPLinearSolverPrepare;
00277   mops->matsolve=DSDPLinearSolve2;
00278   mops->matdestroy=Tdestroy;
00279   mops->matzero=Tzero;
00280   mops->matsetup=Tsetup;
00281   mops->matview=TTTMatView;
00282   mops->id=5;
00283   mops->matname=tmatname;
00284   DSDPFunctionReturn(0);
00285 }
00286 
00287 static struct  DSDPSchurMat_Ops dsdpmatops;
00288 
00289 int DSDPGetDiagSchurMat(int, struct  DSDPSchurMat_Ops **,void **);
00290 int DSDPGetLAPACKSUSchurOps(int,struct DSDPSchurMat_Ops**,void**);
00291 int DSDPGetLAPACKPUSchurOps(int,struct DSDPSchurMat_Ops**,void**);
00292 
00293 static int DSDPCreateM(MCholSolverALL*,chfac**,int[],int[],int);
00294 static int DSDPCreateSchurMatrix(void*,int);
00295 
00296 #undef __FUNCT__
00297 #define __FUNCT__ "DSDPCreateSchurMatrix"
00298 static int DSDPCreateSchurMatrix(void *ctx, int m){
00299 
00300   int info;
00301   int *rnnz,*tnnz,totalnnz;
00302   int gotit=0;
00303   DSDP dsdp=(DSDP)ctx;
00304   chfac *sf;
00305   MCholSolverALL *AMA;
00306   void *tdata;
00307   struct  DSDPSchurMat_Ops *tmatops;
00308 
00309   DSDPFunctionBegin;
00310   if (m<=1){
00311     info=DSDPGetDiagSchurMat(m,&tmatops,&tdata); DSDPCHKERR(info);
00312     info=DSDPSetSchurMatOps(dsdp,tmatops,tdata); DSDPCHKERR(info);
00313     return 0;
00314   }
00315 
00316   DSDPCALLOC2(&rnnz,int,(m+1),&info); DSDPCHKERR(info);
00317   DSDPCALLOC2(&tnnz,int,(m+1),&info); DSDPCHKERR(info);
00318 
00319   info=DSDPCheckForSparsity(dsdp,m,rnnz,tnnz,&totalnnz); DSDPCHKERR(info);
00320 
00321   if (totalnnz*2+m > m*m*0.1 && dsdpuselapack) {
00322     gotit=1;
00323     info=DSDPGetLAPACKSUSchurOps(m,&tmatops,&tdata); 
00324     /* DSDPCHKERR(info); */ if (info) {gotit=0;printf("Try packed format\n"); }
00325     DSDPLogInfo(0,8,"Creating Dense Full LAPACK Schur Matrix\n");
00326     info=DSDPSetSchurMatOps(dsdp,tmatops,tdata); DSDPCHKERR(info);
00327   } 
00328 
00329   if ( 0 && totalnnz*2+m > m*m*0.1 && dsdpuselapack) {
00330 
00331     info=DSDPGetLAPACKPUSchurOps(m,&tmatops,&tdata); DSDPCHKERR(info);
00332     DSDPLogInfo(0,8,"Creating Dense Packed LAPACK Schur Matrix\n");
00333     info=DSDPSetSchurMatOps(dsdp,tmatops,tdata); DSDPCHKERR(info);
00334     gotit=1;
00335 
00336   } 
00337   if (gotit==0){
00338 
00339     DSDPCALLOC1(&AMA,MCholSolverALL,&info);DSDPCHKERR(info);
00340     AMA->dsdp=dsdp; AMA->m=m;
00341     info=DSDPVecCreateSeq(m,&AMA->D1); DSDPCHKERR(info);
00342     if (totalnnz*2+m > m*m * 0.11 ){
00343       info=MchlSetup2(m,&sf); DSDPCHKERR(info);
00344       AMA->M=sf;      AMA->isdense=1;
00345       AMA->rnnz=0;    AMA->colnnz=0;
00346       DSDPLogInfo(0,8,"Creating Dense Full non LAPACK Schur Matrix\n");
00347     } else {
00348       info=DSDPCreateM(AMA,&sf,rnnz,tnnz,totalnnz); DSDPCHKERR(info);
00349       DSDPLogInfo(0,8,"Creating Sparse Schur Matrix\n");
00350     }
00351     AMA->M=sf;  
00352     info=TMatOpsInit(&dsdpmatops); DSDPCHKERR(info);
00353     info=DSDPSetSchurMatOps(dsdp,&dsdpmatops,(void*)AMA); DSDPCHKERR(info);
00354   }
00355   DSDPFREE(&tnnz,&info);DSDPCHKERR(info);
00356   DSDPFREE(&rnnz,&info);DSDPCHKERR(info);
00357   DSDPFunctionReturn(0);
00358 }
00359 
00360 static struct  DSDPSchurMat_Ops dsdpmatops000;
00361 
00362 #undef __FUNCT__
00363 #define __FUNCT__ "DSDPSetDefaultSchurMatrixStructure"
00364 int DSDPSetDefaultSchurMatrixStructure(DSDP dsdp){
00365   int info;
00366   DSDPFunctionBegin;
00367   info=DSDPSchurMatOpsInitialize(&dsdpmatops000); DSDPCHKERR(info);
00368   dsdpmatops000.matsetup=DSDPCreateSchurMatrix;
00369   info=DSDPSetSchurMatOps(dsdp,&dsdpmatops000,(void*)dsdp);DSDPCHKERR(info);
00370   DSDPFunctionReturn(0);
00371 }
00372  

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