diag.c

Go to the documentation of this file.
00001 
00002 #include "dsdpschurmat_impl.h"
00003 #include "dsdpdualmat_impl.h"
00004 #include "dsdpdsmat_impl.h"
00005 #include "dsdpsys.h"
00012 typedef struct {
00013   int    n;
00014   double *val;
00015   int   owndata;
00016 } diagmat;
00017 
00018 static int DiagMatCreate(int,diagmat**);
00019 static int DiagMatMult(void*,double[],double[],int);
00020 static int DiagMatGetSize(void*, int *);
00021 static int DiagMatAddRow2(void*, int, double, double[], int);
00022 static int DiagMatDestroy(void*);
00023 static int DiagMatView(void*);
00024 static int DiagMatLogDeterminant(void*, double *);
00025 
00026 /* static int DiagMatScale(double *, diagmat *); */
00027 
00028 static int DiagMatCreate(int n,diagmat**M){
00029   int info;
00030   diagmat *M7;
00031 
00032   DSDPCALLOC1(&M7,diagmat,&info);DSDPCHKERR(info);
00033   DSDPCALLOC2(&M7->val,double,n,&info);DSDPCHKERR(info);
00034   if (n>0 && M7->val==NULL) return 1;
00035   M7->owndata=1; M7->n=n;
00036   *M=M7;
00037   return 0;
00038 }
00039 
00040 static int DiagMatDestroy(void* AA){
00041   int info;
00042   diagmat* A=(diagmat*) AA;
00043   if (A->owndata && A->val){ DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
00044   DSDPFREE(&A,&info);DSDPCHKERR(info);
00045   return 0;
00046 }
00047 
00048 
00049 static int DiagMatMult(void* AA, double x[], double y[], int n){
00050 
00051   diagmat* A=(diagmat*) AA;
00052   double *vv=A->val;
00053   int i;
00054 
00055   if (A->n != n) return 1;
00056   if (x==0 && n>0) return 3;
00057   if (y==0 && n>0) return 3;
00058 
00059   for (i=0; i<n; i++){
00060     y[i]=x[i]*vv[i];
00061   }  
00062   return 0;
00063 }
00064 
00065 
00066 static int DiagMatGetSize(void *AA, int *n){
00067   diagmat* A=(diagmat*) AA;
00068   *n=A->n;
00069   return 0;
00070 }
00071 
00072 static int DiagMatView(void* AA){
00073   diagmat* A=(diagmat*) AA;
00074   int i;
00075   for (i=0;i<A->n; i++){
00076     printf(" Row: %d, Column: %d, Value: %8.4e \n",i,i,A->val[i]);
00077   }
00078   return 0;
00079 }
00080 
00081 static int DiagMatAddRow2(void* AA, int nrow, double dd, double row[], int n){
00082   diagmat* A=(diagmat*) AA;
00083   A->val[nrow] += dd*row[nrow];
00084   return 0;
00085 }
00086 
00087 
00088 static int DiagMatAddElement(void*A, int nrow, double dd){
00089   diagmat* AA = (diagmat*)A;
00090   AA->val[nrow] += dd;
00091   return 0;
00092 }
00093 
00094 static int DiagMatZeros(void*A){
00095   diagmat* AA = (diagmat*)A;
00096   int n=AA->n;
00097   memset(AA->val,0,n*sizeof(double));
00098   return 0;
00099 }
00100 
00101 static int DiagMatSolve(void* A, double b[], double x[],int n){
00102   diagmat* AA = (diagmat*)A;
00103   double *v=AA->val;
00104   int i;
00105   for (i=0;i<n;i++){
00106     x[i]=b[i]/v[i];
00107   }
00108   return 0;
00109 }
00110 
00111 static int DiagMatSolve2(void* A, int indx[], int nindx, double b[], double x[],int n){
00112   diagmat* AA = (diagmat*)A;
00113   double *v=AA->val;
00114   int i,j;
00115   memset((void*)x,0,n*sizeof(double));
00116   for (j=0;j<nindx;j++){
00117     i=indx[j];
00118     x[i]=b[i]/v[i];
00119   }
00120   return 0;
00121 }
00122 
00123 static int DiagMatCholeskySolveBackward(void* A, double b[], double x[],int n){
00124   int i;
00125   for (i=0;i<n;i++){
00126     x[i]=b[i];
00127   }
00128   return 0;
00129 }
00130 
00131 static int DiagMatInvert(void *A){
00132   return 0;
00133 }
00134 
00135 static int DiagMatCholeskyFactor(void*A,int *flag){
00136   diagmat* AA = (diagmat*)A;
00137   double *v=AA->val;
00138   int i,n=AA->n;
00139   *flag=0;
00140   for (i=0;i<n;i++){
00141     if (v[i]<=0){ *flag=i+1; break;}
00142   }
00143   return 0;
00144 }
00145 
00146 static int DiagMatLogDeterminant(void*A, double *dd){
00147   diagmat* AA = (diagmat*)A;
00148   double d=0,*val=AA->val;
00149   int i,n=AA->n;
00150   for (i=0;i<n;i++){
00151     if (val[i]<=0) return 1;
00152     d+=log(val[i]);
00153   }
00154   *dd=d;
00155   return 0;
00156 }
00157 
00158 static int DiagMatTakeUREntriesP(void*A, double dd[], int nn, int n){
00159   diagmat* AA = (diagmat*)A;
00160   int i,ii;
00161   double *val=AA->val;
00162   for (i=0;i<n;i++){
00163     ii=(i+1)*(i+2)/2-1;
00164     val[i]=dd[ii];
00165   }
00166   return 0;
00167 }
00168 static int DiagMatTakeUREntriesU(void*A, double dd[], int nn, int n){
00169   diagmat* AA = (diagmat*)A;
00170   int i;
00171   double *val=AA->val;
00172   for (i=0;i<n;i++){
00173     val[i]=dd[i*n+i];
00174   }
00175   return 0;
00176 }
00177 
00178 static int DiagMatInverseAddP(void*A, double alpha, double dd[], int nn, int n){
00179   diagmat* AA = (diagmat*)A;
00180   int i,ii;
00181   double *val=AA->val;
00182   for (i=0;i<n;i++){
00183     ii=(i+1)*(i+2)/2-1;
00184     dd[ii]+=alpha/val[i];
00185   }
00186   return 0;
00187 }
00188 static int DiagMatInverseAddU(void*A, double alpha, double dd[], int nn, int n){
00189   diagmat* AA = (diagmat*)A;
00190   int i;
00191   double *val=AA->val;
00192   for (i=0;i<n;i++){
00193     dd[i*n+i]+=alpha/val[i];
00194   }
00195   return 0;
00196 }
00197 
00198 static int DiagMatFull(void*A, int* dfull){
00199   *dfull=1;
00200   return 0;
00201 }
00202 
00203 static struct  DSDPDualMat_Ops sdmatopsp;
00204 static struct  DSDPDualMat_Ops sdmatopsu;
00205 static const char* diagmatname="DIAGONAL";
00206 
00207 static int DiagDualOpsInitializeP(struct  DSDPDualMat_Ops* sops){
00208   int info;
00209   if (sops==NULL) return 0;
00210   info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
00211   sops->matcholesky=DiagMatCholeskyFactor;
00212   sops->matsolveforward=DiagMatSolve;
00213   sops->matsolvebackward=DiagMatCholeskySolveBackward;
00214   sops->matinvert=DiagMatInvert;
00215   sops->matinverseadd=DiagMatInverseAddP;
00216   sops->matinversemultiply=DiagMatSolve2;
00217   sops->matseturmat=DiagMatTakeUREntriesP;
00218   sops->matfull=DiagMatFull;
00219   sops->matdestroy=DiagMatDestroy;
00220   sops->matgetsize=DiagMatGetSize;
00221   sops->matview=DiagMatView;
00222   sops->matlogdet=DiagMatLogDeterminant;
00223   sops->id=9;
00224   sops->matname=diagmatname;
00225   return 0;
00226 }
00227 static int DiagDualOpsInitializeU(struct  DSDPDualMat_Ops* sops){
00228   int info;
00229   if (sops==NULL) return 0;
00230   info=DSDPDualMatOpsInitialize(sops); DSDPCHKERR(info);
00231   sops->matcholesky=DiagMatCholeskyFactor;
00232   sops->matsolveforward=DiagMatSolve;
00233   sops->matsolvebackward=DiagMatCholeskySolveBackward;
00234   sops->matinvert=DiagMatInvert;
00235   sops->matinversemultiply=DiagMatSolve2;
00236   sops->matseturmat=DiagMatTakeUREntriesU;
00237   sops->matfull=DiagMatFull;
00238   sops->matinverseadd=DiagMatInverseAddU;
00239   sops->matdestroy=DiagMatDestroy;
00240   sops->matgetsize=DiagMatGetSize;
00241   sops->matview=DiagMatView;
00242   sops->matlogdet=DiagMatLogDeterminant;
00243   sops->id=9;
00244   sops->matname=diagmatname;
00245   return 0;
00246 }
00247 
00248 #undef __FUNCT__
00249 #define __FUNCT__ "DSDPDiagDualMatCreateP"
00250 int DSDPDiagDualMatCreateP(int n, 
00251                            struct  DSDPDualMat_Ops **sops1, void**smat1,
00252                            struct  DSDPDualMat_Ops **sops2, void**smat2){
00253   diagmat *AA;
00254   int info;
00255   DSDPFunctionBegin;
00256 
00257   info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00258   info=DiagDualOpsInitializeP(&sdmatopsp); DSDPCHKERR(info);
00259   *sops1=&sdmatopsp;
00260   *smat1=(void*)AA;
00261 
00262   info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00263   info=DiagDualOpsInitializeP(&sdmatopsp); DSDPCHKERR(info);
00264   *sops2=&sdmatopsp;
00265   *smat2=(void*)AA;
00266   DSDPFunctionReturn(0);
00267 }
00268 
00269 #undef __FUNCT__
00270 #define __FUNCT__ "DSDPDiagDualMatCreateU"
00271 int DSDPDiagDualMatCreateU(int n, 
00272                            struct  DSDPDualMat_Ops **sops1, void**smat1,
00273                            struct  DSDPDualMat_Ops **sops2, void**smat2){
00274   diagmat *AA;
00275   int info;
00276   DSDPFunctionBegin;
00277   info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00278   info=DiagDualOpsInitializeU(&sdmatopsu); DSDPCHKERR(info);
00279   *sops1=&sdmatopsu;
00280   *smat1=(void*)AA;
00281   info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00282   info=DiagDualOpsInitializeU(&sdmatopsu); DSDPCHKERR(info);
00283   *sops2=&sdmatopsu;
00284   *smat2=(void*)AA;
00285   DSDPFunctionReturn(0);
00286 }
00287 
00288 
00289 
00290 static int DiagMatVecVec(void*A,double x[], int n, double *vv){
00291   diagmat* AA = (diagmat*)A;
00292   double *v=AA->val,vAv=0;
00293   int i;
00294   for (i=0;i<n;i++){
00295     vAv+=x[i]*x[i]*v[i];
00296   }
00297   *vv=vAv;
00298   return 0;
00299 }
00300 
00301 static int DDiagDSMatOpsInitP(struct  DSDPDSMat_Ops *ddiagops){
00302   int info;
00303   if (ddiagops==NULL) return 0;
00304   info=DSDPDSMatOpsInitialize(ddiagops);DSDPCHKERR(info);
00305   ddiagops->matseturmat=DiagMatTakeUREntriesP;
00306   ddiagops->matview=DiagMatView;
00307   ddiagops->matgetsize=DiagMatGetSize;
00308   ddiagops->matmult=DiagMatMult;
00309   ddiagops->matvecvec=DiagMatVecVec;
00310   ddiagops->matzeroentries=DiagMatZeros;
00311   ddiagops->matdestroy=DiagMatDestroy;
00312   ddiagops->id=9;
00313   ddiagops->matname=diagmatname;
00314   DSDPFunctionReturn(0);
00315 }
00316 static int DDiagDSMatOpsInitU(struct  DSDPDSMat_Ops *ddiagops){
00317   int info;
00318   if (ddiagops==NULL) return 0;
00319   info=DSDPDSMatOpsInitialize(ddiagops);DSDPCHKERR(info);
00320   ddiagops->matseturmat=DiagMatTakeUREntriesU;
00321   ddiagops->matview=DiagMatView;
00322   ddiagops->matgetsize=DiagMatGetSize;
00323   ddiagops->matmult=DiagMatMult;
00324   ddiagops->matvecvec=DiagMatVecVec;
00325   ddiagops->matzeroentries=DiagMatZeros;
00326   ddiagops->matdestroy=DiagMatDestroy;
00327   ddiagops->id=9;
00328   ddiagops->matname=diagmatname;
00329   DSDPFunctionReturn(0);
00330 }
00331 
00332 static struct  DSDPDSMat_Ops dsdiagmatopsp;
00333 static struct  DSDPDSMat_Ops dsdiagmatopsu;
00334 
00335 #undef __FUNCT__
00336 #define __FUNCT__ "DSDPDiagDSMatP"
00337 int DSDPCreateDiagDSMatP(int n,struct  DSDPDSMat_Ops* *dsmatops, void**dsmat){
00338 
00339   int info=0;
00340   diagmat *AA;
00341 
00342   DSDPFunctionBegin;
00343   info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00344   info=DDiagDSMatOpsInitP(&dsdiagmatopsp); DSDPCHKERR(info);
00345   *dsmatops=&dsdiagmatopsp;
00346   *dsmat=(void*)AA;
00347   DSDPFunctionReturn(0);
00348 }
00349 #undef __FUNCT__
00350 #define __FUNCT__ "DSDPDiagDSMatU"
00351 int DSDPCreateDiagDSMatU(int n,struct  DSDPDSMat_Ops* *dsmatops, void**dsmat){
00352 
00353   int info=0;
00354   diagmat *AA;
00355 
00356   DSDPFunctionBegin;
00357   info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00358   info=DDiagDSMatOpsInitU(&dsdiagmatopsu); DSDPCHKERR(info);
00359   *dsmatops=&dsdiagmatopsu;
00360   *dsmat=(void*)AA;
00361   DSDPFunctionReturn(0);
00362 }
00363 
00364 static int DiagRowNonzeros(void*M, int row, double cols[], int *ncols,int nrows){
00365   DSDPFunctionBegin;
00366   *ncols = 1;
00367   cols[row]=1.0;
00368   DSDPFunctionReturn(0);
00369 }
00370 
00371 static int DiagAddDiag(void*M, double diag[], int m){
00372   diagmat* AA = (diagmat*)M;
00373   double *v=AA->val;
00374   int i;
00375   DSDPFunctionBegin;
00376   for (i=0;i<m;i++){
00377     v[i]+=diag[i];
00378   }
00379   DSDPFunctionReturn(0);
00380 }
00381 
00382 static int DiagMultiply(void*M, double xin[], double xout[], int m){
00383   diagmat* AA = (diagmat*)M;
00384   double *v=AA->val;
00385   int i;
00386   DSDPFunctionBegin;
00387   for (i=0;i<m;i++){
00388     xout[i]+=v[i]*xin[i];
00389   }
00390   DSDPFunctionReturn(0);
00391 }
00392 
00393 static int DiagShiftDiag(void*M, double dd){
00394   diagmat* AA = (diagmat*)M;
00395   double *v=AA->val;
00396   int i,m=AA->n;
00397   DSDPFunctionBegin;
00398   for (i=0;i<m;i++){
00399     v[i]+=dd;
00400   }
00401   DSDPFunctionReturn(0);
00402 }
00403 
00404 static int DiagAddElement(void*M, int ii, double dd){
00405   diagmat* AA = (diagmat*)M;
00406   DSDPFunctionBegin;
00407   AA->val[ii]+=dd;
00408   DSDPFunctionReturn(0);
00409 }
00410 
00411 static int DiagMatOnProcessor(void*A,int row,int*flag){
00412   *flag=1;
00413   return 0;
00414 }
00415 
00416 static int DiagAssemble(void*M){
00417   return 0;
00418 }
00419 
00420 static struct  DSDPSchurMat_Ops dsdpdiagschurops;
00421 
00422 #undef __FUNCT__
00423 #define __FUNCT__ "DSDPDiagSchurOps"
00424 static int DiagSchurOps(struct  DSDPSchurMat_Ops *sops){
00425   int info;
00426   DSDPFunctionBegin;
00427   if (!sops) return 0;
00428   info=DSDPSchurMatOpsInitialize(sops); DSDPCHKERR(info);
00429   sops->matzero=DiagMatZeros;
00430   sops->mataddrow=DiagMatAddRow2;
00431   sops->mataddelement=DiagMatAddElement;
00432   sops->matdestroy=DiagMatDestroy;
00433   sops->matfactor=DiagMatCholeskyFactor;
00434   sops->matsolve=DiagMatSolve;
00435   sops->matadddiagonal=DiagAddDiag;
00436   sops->matshiftdiagonal=DiagShiftDiag;
00437   sops->mataddelement=DiagAddElement;
00438   sops->matscaledmultiply=DiagMultiply;
00439   sops->matassemble=DiagAssemble;
00440   sops->pmatonprocessor=DiagMatOnProcessor;
00441   sops->matrownonzeros=DiagRowNonzeros;
00442   sops->id=9;
00443   sops->matname=diagmatname;
00444   DSDPFunctionReturn(0);
00445 }
00446 
00447 #undef __FUNCT__
00448 #define __FUNCT__ "DSDPGetDiagSchurMat"
00449 int DSDPGetDiagSchurMat(int n, struct  DSDPSchurMat_Ops **sops, void **data){
00450   int info=0;
00451   diagmat *AA;
00452   DSDPFunctionBegin;
00453   info=DiagMatCreate(n,&AA); DSDPCHKERR(info);
00454   info=DiagSchurOps(&dsdpdiagschurops); DSDPCHKERR(info);
00455   if (sops){*sops=&dsdpdiagschurops;}
00456   if (data){*data=(void*)AA;}
00457   DSDPFunctionReturn(0);
00458 }

Generated on Wed Nov 5 21:46:05 2008 for DSDP by  doxygen 1.5.6