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