drowcol.c

Go to the documentation of this file.
00001 
00002 typedef struct {
00003   int    rc;
00004   const double *val;
00005   int    n;
00006   double x,y;
00007 } rcmat;
00008 
00009 #include "dsdpdatamat_impl.h"
00010 #include "dsdpsys.h"
00017 static int RCMatDestroy(void*);
00018 static int RCMatView(void*);
00019 static int RCMatVecVec(void*, double[], int, double *);
00020 static int RCMatDot(void*, double[], int, int, double *);
00021 static int RCMatGetRank(void*, int*, int);
00022 static int RCMatFactor(void*);
00023 static int RCMatGetEig(void*, int, double*, double[], int, int[], int*);
00024 static int RCMatAddRowMultiple(void*, int, double, double[], int);
00025 static int RCMatAddMultiple(void*, double, double[], int, int);
00026 static int RCMatGetRowNnz(void*, int, int[], int*, int);
00027 
00028 static int RCMatCreate(int n,int rowcol, const double v[], rcmat**M){
00029   rcmat *M16;
00030   M16=(rcmat*) malloc(1*sizeof(rcmat));
00031   M16->val=v;
00032   M16->rc=rowcol;
00033   M16->n=n;
00034   *M=M16;
00035   return 0;
00036 }
00037 
00038 static int RCMatDestroy(void* AA){
00039   free(AA);
00040   return 0;
00041 }
00042 static int RCMatVecVec(void* A, double x[], int n, double *v){
00043   rcmat*RC=(rcmat*)A;
00044   int i;
00045   double vv=0;
00046   const double *val=RC->val;
00047   for (i=0;i<n;i++){ vv+=val[i]*x[i];}
00048   *v=2*vv*x[RC->rc];
00049   return 0;
00050 }
00051 static int RCMatDot(void* A, double x[], int nn, int n1, double *v){
00052   rcmat*RC=(rcmat*)A;
00053   int i,k=0,rc=RC->rc,n=RC->n;
00054   double vv=0;
00055   const double *val=RC->val;
00056   k=rc*(rc+1)/2;
00057   for (i=0;i<=rc;i++,k++){ vv+=v[k]*val[i]; }
00058   for (i=rc+1;i<n;i++,k+=i){ vv+=v[k+rc]*val[i];}
00059   *v=2*vv;
00060   return 0;
00061 }
00062 static int RCMatView(void* A){
00063   rcmat*RC=(rcmat*)A;
00064   int i;
00065   printf("Row Col %d\n",RC->rc);
00066   for (i=0;i<RC->n;i++){
00067     printf("%4.4e ",RC->val[i]);
00068   }
00069   printf("\n");
00070   return 0;
00071 }
00072 
00073 static int RCMatFactor(void* A){
00074   /* Must compute x,y such that eigs are orthogonal */
00075   rcmat*RC=(rcmat*)A;
00076   int i,rc=RC->rc;
00077   double vnorm2=0;
00078   const double *val=RC->val;
00079   for (i=0;i<RC->n;i++){ vnorm2+=val[i]*val[i];} 
00080   vnorm2=sqrt(vnorm2);
00081   if (val[rc]>=0){
00082     RC->y=-sqrt(vnorm2);
00083     RC->x=sqrt(2*val[rc]+RC->y*RC->y);
00084   } else {
00085     RC->x=sqrt(vnorm2);
00086     RC->y=-sqrt(-2*val[rc]+RC->x*RC->x);
00087   }
00088   return 0;
00089 }
00090 static int RCMatGetRank(void *A, int*rank, int n){
00091   *rank=2;
00092   return 0;
00093 }
00094 
00095 static int RCMatGetEig(void*A, int neig, double *eig, double v[], int n,int spind[], int *nind){
00096   rcmat*RC=(rcmat*)A;
00097   int i,rc=RC->rc;
00098   double x=RC->x,y=RC->y,xmy=x-y;
00099   const double *val=RC->val;
00100   if (neig==0){
00101     for (i=0;i<n;i++){spind[i]=i;}
00102     for (i=0;i<n;i++){v[i]=val[i]/xmy;}
00103     v[rc]=RC->x;
00104     *nind=n;
00105     *eig=1.0;
00106   } else if (neig==1){
00107     for (i=0;i<n;i++){spind[i]=i;}
00108     for (i=0;i<n;i++){v[i]=val[i]/xmy;}
00109     v[rc]=RC->y;
00110     *nind=n;
00111     *eig=-1.0;
00112   } else { *eig=0;*nind=0;}
00113   return 0;
00114 }
00115 
00116 static int RCMatGetRowNnz(void*A, int nrow, int nz[], int *nnzz, int n){
00117   rcmat*RC=(rcmat*)A;
00118   int i;
00119   *nnzz=1;
00120   if (nrow==RC->rc){for (i=0;i<n;i++){nz[i]++;}*nnzz=n;}
00121   nz[nrow]++;
00122   return 0;
00123 }
00124 
00125 static int RCMatAddRowMultiple(void*A, int nrow, double dd, double rrow[], int n){
00126   rcmat*RC=(rcmat*)A;
00127   int i;
00128   if (nrow==RC->rc){ 
00129     for (i=0;i<n;i++){rrow[i]+=dd*RC->val[i];}
00130   }
00131   rrow[nrow]+=dd*RC->val[nrow];
00132   return 0;
00133 }
00134 static int RCMatAddMultiple(void*A, double dd, double vv[], int nn, int n1){
00135   rcmat*RC=(rcmat*)A;
00136   int i,rc=RC->rc,k=rc*(rc+1)/2,n=RC->n;
00137   const double *val=RC->val;
00138   for (i=0;i<=rc;i++,k++){
00139     vv[k]+=dd*val[i];
00140   }
00141   for (i=rc+1;i<n;i++,k+=i){
00142     vv[k+rc]+=dd*val[i];
00143   }
00144   return 0;
00145 }
00146 static int RCMatFNorm(void*A, int n, double *fnorm){
00147   rcmat*RC=(rcmat*)A;
00148   int i,rc=RC->rc;
00149   double ff=0;
00150   const double *val=RC->val;
00151   for (i=0;i<n;i++){
00152     ff+=val[i]*val[i];
00153   }
00154   ff*=2;
00155   ff+=2*val[rc]*val[rc];
00156   *fnorm=ff;
00157   return 0;
00158 }
00159 static int RCMatCountNonzeros(void*A, int *nnz, int n){
00160   *nnz=2*n-1;
00161   return 0;
00162 }
00163 static struct  DSDPDataMat_Ops rcmatops;
00164 static const char *datamatname="One Row and Column matrix";
00165 static int RCMatOperationsInitialize(struct  DSDPDataMat_Ops* rcmatoperator){
00166   int info;
00167   if (rcmatoperator==NULL) return 0;
00168   info=DSDPDataMatOpsInitialize(rcmatoperator); if (info){ return info;}
00169   rcmatoperator->matfactor1=RCMatFactor;
00170   rcmatoperator->matgetrank=RCMatGetRank;
00171   rcmatoperator->matgeteig=RCMatGetEig;
00172   rcmatoperator->matvecvec=RCMatVecVec;
00173   rcmatoperator->matrownz=RCMatGetRowNnz;
00174   rcmatoperator->matdot=RCMatDot;
00175   rcmatoperator->matfnorm2=RCMatFNorm;
00176   rcmatoperator->matnnz=RCMatCountNonzeros;
00177   rcmatoperator->mataddrowmultiple=RCMatAddRowMultiple;
00178   rcmatoperator->mataddallmultiple=RCMatAddMultiple;
00179   rcmatoperator->matdestroy=RCMatDestroy;
00180   rcmatoperator->matview=RCMatView;
00181   rcmatoperator->matname=datamatname;
00182   rcmatoperator->id=27;
00183   return 0;
00184 }
00185 
00186 
00187 #undef __FUNCT__
00188 #define __FUNCT__ "DSDPGetRCMat"
00189 int DSDPGetRCMat(int n,const double val[],int rowcol, struct  DSDPDataMat_Ops**sops, void**smat){
00190   rcmat *AA;
00191   int info;
00192   DSDPFunctionBegin;
00193   info=RCMatCreate(n,rowcol,val,&AA);DSDPCHKERR(info);
00194   info=RCMatOperationsInitialize(&rcmatops); DSDPCHKERR(info);
00195   if (sops){*sops=&rcmatops;}
00196   if (smat){*smat=(void*)AA;}
00197   DSDPFunctionReturn(0);
00198 }
00199 

Generated on Wed Nov 11 20:41:02 2009 for DSDP by  doxygen 1.6.1