DSDP
|
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