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