DSDP
drowcol.c
Go to the documentation of this file.
1 
2 typedef struct {
3  int rc;
4  const double *val;
5  int n;
6  double x,y;
7 } rcmat;
8 
9 #include "dsdpdatamat_impl.h"
10 #include "dsdpsys.h"
17 static int RCMatDestroy(void*);
18 static int RCMatView(void*);
19 static int RCMatVecVec(void*, double[], int, double *);
20 static int RCMatDot(void*, double[], int, int, double *);
21 static int RCMatGetRank(void*, int*, int);
22 static int RCMatFactor(void*);
23 static int RCMatGetEig(void*, int, double*, double[], int, int[], int*);
24 static int RCMatAddRowMultiple(void*, int, double, double[], int);
25 static int RCMatAddMultiple(void*, double, double[], int, int);
26 static int RCMatGetRowNnz(void*, int, int[], int*, int);
27 
28 static int RCMatCreate(int n,int rowcol, const double v[], rcmat**M){
29  rcmat *M16;
30  M16=(rcmat*) malloc(1*sizeof(rcmat));
31  M16->val=v;
32  M16->rc=rowcol;
33  M16->n=n;
34  *M=M16;
35  return 0;
36 }
37 
38 static int RCMatDestroy(void* AA){
39  free(AA);
40  return 0;
41 }
42 static int RCMatVecVec(void* A, double x[], int n, double *v){
43  rcmat*RC=(rcmat*)A;
44  int i;
45  double vv=0;
46  const double *val=RC->val;
47  for (i=0;i<n;i++){ vv+=val[i]*x[i];}
48  *v=2*vv*x[RC->rc];
49  return 0;
50 }
51 static int RCMatDot(void* A, double x[], int nn, int n1, double *v){
52  rcmat*RC=(rcmat*)A;
53  int i,k=0,rc=RC->rc,n=RC->n;
54  double vv=0;
55  const double *val=RC->val;
56  k=rc*(rc+1)/2;
57  for (i=0;i<=rc;i++,k++){ vv+=v[k]*val[i]; }
58  for (i=rc+1;i<n;i++,k+=i){ vv+=v[k+rc]*val[i];}
59  *v=2*vv;
60  return 0;
61 }
62 static int RCMatView(void* A){
63  rcmat*RC=(rcmat*)A;
64  int i;
65  printf("Row Col %d\n",RC->rc);
66  for (i=0;i<RC->n;i++){
67  printf("%4.4e ",RC->val[i]);
68  }
69  printf("\n");
70  return 0;
71 }
72 
73 static int RCMatFactor(void* A){
74  /* Must compute x,y such that eigs are orthogonal */
75  rcmat*RC=(rcmat*)A;
76  int i,rc=RC->rc;
77  double vnorm2=0;
78  const double *val=RC->val;
79  for (i=0;i<RC->n;i++){ vnorm2+=val[i]*val[i];}
80  vnorm2=sqrt(vnorm2);
81  if (val[rc]>=0){
82  RC->y=-sqrt(vnorm2);
83  RC->x=sqrt(2*val[rc]+RC->y*RC->y);
84  } else {
85  RC->x=sqrt(vnorm2);
86  RC->y=-sqrt(-2*val[rc]+RC->x*RC->x);
87  }
88  return 0;
89 }
90 static int RCMatGetRank(void *A, int*rank, int n){
91  *rank=2;
92  return 0;
93 }
94 
95 static int RCMatGetEig(void*A, int neig, double *eig, double v[], int n,int spind[], int *nind){
96  rcmat*RC=(rcmat*)A;
97  int i,rc=RC->rc;
98  double x=RC->x,y=RC->y,xmy=x-y;
99  const double *val=RC->val;
100  if (neig==0){
101  for (i=0;i<n;i++){spind[i]=i;}
102  for (i=0;i<n;i++){v[i]=val[i]/xmy;}
103  v[rc]=RC->x;
104  *nind=n;
105  *eig=1.0;
106  } else if (neig==1){
107  for (i=0;i<n;i++){spind[i]=i;}
108  for (i=0;i<n;i++){v[i]=val[i]/xmy;}
109  v[rc]=RC->y;
110  *nind=n;
111  *eig=-1.0;
112  } else { *eig=0;*nind=0;}
113  return 0;
114 }
115 
116 static int RCMatGetRowNnz(void*A, int nrow, int nz[], int *nnzz, int n){
117  rcmat*RC=(rcmat*)A;
118  int i;
119  *nnzz=1;
120  if (nrow==RC->rc){for (i=0;i<n;i++){nz[i]++;}*nnzz=n;}
121  nz[nrow]++;
122  return 0;
123 }
124 
125 static int RCMatAddRowMultiple(void*A, int nrow, double dd, double rrow[], int n){
126  rcmat*RC=(rcmat*)A;
127  int i;
128  if (nrow==RC->rc){
129  for (i=0;i<n;i++){rrow[i]+=dd*RC->val[i];}
130  }
131  rrow[nrow]+=dd*RC->val[nrow];
132  return 0;
133 }
134 static int RCMatAddMultiple(void*A, double dd, double vv[], int nn, int n1){
135  rcmat*RC=(rcmat*)A;
136  int i,rc=RC->rc,k=rc*(rc+1)/2,n=RC->n;
137  const double *val=RC->val;
138  for (i=0;i<=rc;i++,k++){
139  vv[k]+=dd*val[i];
140  }
141  for (i=rc+1;i<n;i++,k+=i){
142  vv[k+rc]+=dd*val[i];
143  }
144  return 0;
145 }
146 static int RCMatFNorm(void*A, int n, double *fnorm){
147  rcmat*RC=(rcmat*)A;
148  int i,rc=RC->rc;
149  double ff=0;
150  const double *val=RC->val;
151  for (i=0;i<n;i++){
152  ff+=val[i]*val[i];
153  }
154  ff*=2;
155  ff+=2*val[rc]*val[rc];
156  *fnorm=ff;
157  return 0;
158 }
159 static int RCMatCountNonzeros(void*A, int *nnz, int n){
160  *nnz=2*n-1;
161  return 0;
162 }
163 static struct DSDPDataMat_Ops rcmatops;
164 static const char *datamatname="One Row and Column matrix";
165 static int RCMatOperationsInitialize(struct DSDPDataMat_Ops* rcmatoperator){
166  int info;
167  if (rcmatoperator==NULL) return 0;
168  info=DSDPDataMatOpsInitialize(rcmatoperator); if (info){ return info;}
169  rcmatoperator->matfactor1=RCMatFactor;
170  rcmatoperator->matgetrank=RCMatGetRank;
171  rcmatoperator->matgeteig=RCMatGetEig;
172  rcmatoperator->matvecvec=RCMatVecVec;
173  rcmatoperator->matrownz=RCMatGetRowNnz;
174  rcmatoperator->matdot=RCMatDot;
175  rcmatoperator->matfnorm2=RCMatFNorm;
176  rcmatoperator->matnnz=RCMatCountNonzeros;
177  rcmatoperator->mataddrowmultiple=RCMatAddRowMultiple;
178  rcmatoperator->mataddallmultiple=RCMatAddMultiple;
179  rcmatoperator->matdestroy=RCMatDestroy;
180  rcmatoperator->matview=RCMatView;
181  rcmatoperator->matname=datamatname;
182  rcmatoperator->id=27;
183  return 0;
184 }
185 
186 
187 #undef __FUNCT__
188 #define __FUNCT__ "DSDPGetRCMat"
189 int DSDPGetRCMat(int n,const double val[],int rowcol, struct DSDPDataMat_Ops**sops, void**smat){
190  rcmat *AA;
191  int info;
192  DSDPFunctionBegin;
193  info=RCMatCreate(n,rowcol,val,&AA);DSDPCHKERR(info);
194  info=RCMatOperationsInitialize(&rcmatops); DSDPCHKERR(info);
195  if (sops){*sops=&rcmatops;}
196  if (smat){*smat=(void*)AA;}
197  DSDPFunctionReturn(0);
198 }
199 
Error handling, printing, and profiling.
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition: dsdpdatamat.c:47
Structure of function pointers that each SDP data matrix type (sparse, dense, constant,...
Table of function pointers that operate on the data matrix.