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
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
00145
00146 Cfact=(cfc_sta)ChlFact(sf,sf->iw,sf->rw,TRUE);
00147 if (CfcOk!=Cfact ){ *flag=1; }
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
00256
00257
00258
00259
00260
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 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