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