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