DSDP
rmmat.c
Go to the documentation of this file.
1 #include "dsdpdatamat_impl.h"
2 #include "dsdpsys.h"
8 typedef struct {
9  double ev;
10  const double *spval;
11  const int *spai;
12  int nnz;
13  int n;
14  int ishift;
15  char UPLQ;
16 } r1mat;
17 
18 static int R1MatDestroy(void*);
19 static int R1MatView(void*);
20 static int R1MatVecVec(void*, double[], int, double *);
21 static int R1MatDotP(void*, double[],int,int,double *);
22 static int R1MatDotU(void*, double[],int,int,double *);
23 static int R1MatGetRank(void*, int*, int);
24 static int R1MatFactor(void*);
25 static int R1MatGetEig(void*, int, double*, double[], int,int[],int*);
26 static int R1MatRowNnz(void*, int, int[], int*, int);
27 static int R1MatAddRowMultiple(void*, int, double, double[], int);
28 static int R1MatAddMultipleP(void*, double, double[], int,int);
29 static int R1MatAddMultipleU(void*, double, double[], int,int);
30 
31 static struct DSDPDataMat_Ops r1matopsP;
32 static struct DSDPDataMat_Ops r1matopsU;
33 static int R1MatOpsInitializeP(struct DSDPDataMat_Ops*);
34 static int R1MatOpsInitializeU(struct DSDPDataMat_Ops*);
35 
36 
37 #undef __FUNCT__
38 #define __FUNCT__ "DSDPGetR1Mat"
39 int DSDPGetR1Mat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, char UPLQ, void**mmat){
40  int i;
41  r1mat*AA;
42  DSDPFunctionBegin;
43  for (i=0;i<nnz;i++){
44  if (spai[i]-ishift<0 || spai[i]-ishift >=n){
45  printf("Invalid entry: Entry %d . Is %d <= %d < %d?\n",i,ishift,spai[i],n+ishift);
46  return 1;
47  }
48  }
49  AA=(r1mat*) malloc(1*sizeof(r1mat));
50  if (AA==NULL) return 1;
51  AA->n=n;
52  AA->UPLQ=UPLQ;
53  AA->spval=spval;
54  AA->spai=spai;
55  AA->nnz=nnz;
56  AA->ev=ev;
57  AA->ishift=ishift;
58  if (mmat){*mmat=(void*)AA;}
59  DSDPFunctionReturn(0);
60 }
61 
62 #undef __FUNCT__
63 #define __FUNCT__ "DSDPGetR1PMat"
64 
77 int DSDPGetR1PMat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, struct DSDPDataMat_Ops**mops, void**mmat){
78  int info;
79  DSDPFunctionBegin;
80  info=DSDPGetR1Mat(n,ev,ishift,spai,spval,nnz,'P',mmat);
81  info=R1MatOpsInitializeP(&r1matopsP); if(info){return 1;}
82  if (mops){*mops=&r1matopsP;}
83  DSDPFunctionReturn(0);
84 }
85 
86 #undef __FUNCT__
87 #define __FUNCT__ "DSDPGetR1UMat"
88 
101 int DSDPGetR1UMat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, struct DSDPDataMat_Ops**mops, void**mmat){
102  int info;
103  DSDPFunctionBegin;
104  info=DSDPGetR1Mat(n,ev,ishift,spai,spval,nnz,'U',mmat);
105  info=R1MatOpsInitializeU(&r1matopsU); if(info){return 1;}
106  if (mops){*mops=&r1matopsU;}
107  DSDPFunctionReturn(0);
108 }
109 
110 static int R1MatDotP(void* A, double x[], int nn, int n, double *v){
111  r1mat* AA = (r1mat*)A;
112  int i,i2,i3,j,j2;
113  int nnz=AA->nnz,ishift=AA->ishift;
114  const int *ai=AA->spai;
115  double dtmp=0.0,d3;
116  const double *val=AA->spval;
117  for (i=0;i<nnz;i++){
118  d3=val[i];
119  i2=ai[i]-ishift;
120  i3=(i2+1)*i2/2;
121  for (j=0;j<nnz;j++){
122  j2=ai[j]-ishift;
123  if (j2<=i2){
124  dtmp+=2*x[i3+j2]*d3*val[j];
125  }
126  }
127  }
128  *v=dtmp*AA->ev;
129  return 0;
130 }
131 
132 static int R1MatDotU(void* A, double x[], int nn, int n, double *v){
133  r1mat* AA = (r1mat*)A;
134  int i,i2,i3,j,j2;
135  int nnz=AA->nnz,ishift=AA->ishift;
136  const int *ai=AA->spai;
137  const double *val=AA->spval;
138  double dtmp=0.0,d3;
139 
140  for (i=0;i<nnz;i++){
141  d3=val[i];
142  i2=ai[i]-ishift;
143  i3=i2*n;
144  for (j=0;j<nnz;j++){
145  j2=ai[j]-ishift;
146  if (j2<=i2){
147  dtmp+=2*x[i3+j2]*d3*val[j];
148  }
149  }
150  }
151  *v=dtmp*AA->ev;
152  return 0;
153 }
154 
155 static int R1MatVecVec(void* A, double x[], int n, double *v){
156 
157  r1mat* AA = (r1mat*)A;
158  double dtmp=0.0;
159  const double *val=AA->spval;
160  int i,ishift=AA->ishift,nnz=AA->nnz;
161  const int *ai=AA->spai;
162  for (i=0; i<nnz; i++){
163  dtmp+=val[i] * x[ai[i]-ishift];
164  }
165  *v=dtmp*dtmp*AA->ev;
166  return 0;
167 }
168 
169 static int R1MatAddMultipleP(void*A, double dd, double vv[], int nn, int n){
170  r1mat* AA = (r1mat*)A;
171  int i,i2,i3,j,j2;
172  int nnz=AA->nnz,ishift=AA->ishift;
173  const int *ai=AA->spai;
174  const double *val=AA->spval;
175  double d3,ddd=dd*AA->ev;
176  for (i=0;i<nnz;i++){
177  d3=ddd*val[i];
178  i2=ai[i]-ishift;
179  i3=(i2+1)*i2/2;
180  for (j=0;j<nnz;j++){
181  j2=ai[j]-ishift;
182  if (j2<=i2){
183  vv[i3+j2]+=d3*val[j];
184  }
185  }
186  }
187  return 0;
188 }
189 static int R1MatAddMultipleU(void*A, double dd, double vv[], int nn, int n){
190  r1mat* AA = (r1mat*)A;
191  int i,i2,i3,j,j2;
192  int nnz=AA->nnz,ishift=AA->ishift;
193  const int *ai=AA->spai;
194  const double *val=AA->spval;
195  double d3,ddd=dd*AA->ev;
196  for (i=0;i<nnz;i++){
197  d3=ddd*val[i];
198  i2=ai[i]-ishift;
199  i3=i2*n;
200  for (j=0;j<nnz;j++){
201  j2=ai[j]-ishift;
202  if (j2<=i2){
203  vv[i3+j2]+=d3*val[j];
204  }
205  }
206  }
207  return 0;
208 }
209 
210 static int R1MatAddRowMultiple(void*A, int nrow, double dd, double row[], int n){
211  r1mat* AA = (r1mat*)A;
212  int nnz=AA->nnz,ishift=AA->ishift;
213  const int *ai=AA->spai;
214  const double *val=AA->spval;
215  double ddd=dd*AA->ev;
216  int i,j;
217  for (i=0;i<nnz;i++){
218  if (ai[i]-ishift==nrow){
219  for (j=0;j<nnz;j++){
220  row[ai[j]-ishift]+= ddd*val[i]*val[j];
221  }
222  }
223  }
224  return 0;
225 }
226 
227 
228 static int R1MatFactor(void*A){
229  return 0;
230 }
231 
232 
233 static int R1MatGetRank(void *A, int*rank, int n){
234  *rank=1;
235  return 0;
236 }
237 
238 static int R1MatGetEig(void*A, int neig, double *eig, double v[], int n, int indx[], int*nind){
239  r1mat* AA = (r1mat*)A;
240  int i,aii,ishift=AA->ishift,nnz=AA->nnz;
241  const int *ai=AA->spai;
242  const double *val=AA->spval;
243  for (i=0;i<n;i++){ v[i]=0.0; }
244  *eig=0; *nind=0;
245  if (neig==0){
246  for (i=0;i<nnz;i++){
247  aii=ai[i]-ishift;
248  v[aii]=val[i];
249  indx[i]=aii;
250  }
251  *eig=AA->ev; *nind=AA->nnz;
252  }
253  return 0;
254 }
255 
256 
257 static int R1MatRowNnz(void*A, int row, int nz[], int *rnnz, int n){
258  r1mat* AA = (r1mat*)A;
259  int i,j;
260  int nnz=AA->nnz,ishift=AA->ishift;
261  const int *ai=AA->spai;
262  *rnnz=0;
263  for (i=0;i<nnz;i++){
264  if (ai[i]-ishift==row){
265  for (j=0;j<nnz;j++){
266  nz[ai[j]-ishift]++;
267  }
268  }
269  *rnnz=nnz;
270  }
271  return 0;
272 }
273 
274 static int R1MatFNorm2(void*A, int n, double *fnorm2){
275  r1mat* AA = (r1mat*)A;
276  double dd=0;
277  const double *val=AA->spval;
278  int i,nnz=AA->nnz;
279  for (i=0;i<nnz;i++){
280  dd+=val[i]*val[i];
281  }
282  *fnorm2=dd*dd*AA->ev*AA->ev;
283  return 0;
284 }
285 
286 static int R1MatCountNonzeros(void*A, int *nnz, int n){
287  r1mat* AA = (r1mat*)A;
288  *nnz=AA->nnz*AA->nnz;
289  return 0;
290 }
291 
292 
293 static int R1MatView(void* A){
294  int i;
295  r1mat* AA = (r1mat*)A;
296  printf("This matrix is %4.8e times the outer product of \n",AA->ev);
297  for (i=0;i<AA->nnz;i++){
298  printf("%d %4.8e \n",AA->spai[i],AA->spval[i]);
299  }
300  return 0;
301 }
302 
303 
304 static int R1MatDestroy(void* A){
305  if (A) free(A);
306  return 0;
307 }
308 
309 static const char *datamatname="RANK 1 Outer Product";
310 static int R1MatOpsInitializeP(struct DSDPDataMat_Ops* r1matops){
311  int info;
312  if (r1matops==NULL) return 0;
313  info=DSDPDataMatOpsInitialize(r1matops); DSDPCHKERR(info);
314  r1matops->matfactor1=R1MatFactor;
315  r1matops->matgetrank=R1MatGetRank;
316  r1matops->matgeteig=R1MatGetEig;
317  r1matops->matvecvec=R1MatVecVec;
318  r1matops->matdot=R1MatDotP;
319  r1matops->mataddrowmultiple=R1MatAddRowMultiple;
320  r1matops->mataddallmultiple=R1MatAddMultipleP;
321  r1matops->matdestroy=R1MatDestroy;
322  r1matops->matview=R1MatView;
323  r1matops->matrownz=R1MatRowNnz;
324  r1matops->matfnorm2=R1MatFNorm2;
325  r1matops->matnnz=R1MatCountNonzeros;
326  r1matops->id=15;
327  r1matops->matname=datamatname;
328  return 0;
329 }
330 static int R1MatOpsInitializeU(struct DSDPDataMat_Ops* r1matops){
331  int info;
332  if (r1matops==NULL) return 0;
333  info=DSDPDataMatOpsInitialize(r1matops); DSDPCHKERR(info);
334  r1matops->matfactor1=R1MatFactor;
335  r1matops->matgetrank=R1MatGetRank;
336  r1matops->matgeteig=R1MatGetEig;
337  r1matops->matvecvec=R1MatVecVec;
338  r1matops->matdot=R1MatDotU;
339  r1matops->mataddrowmultiple=R1MatAddRowMultiple;
340  r1matops->mataddallmultiple=R1MatAddMultipleU;
341  r1matops->matdestroy=R1MatDestroy;
342  r1matops->matview=R1MatView;
343  r1matops->matrownz=R1MatRowNnz;
344  r1matops->matfnorm2=R1MatFNorm2;
345  r1matops->matnnz=R1MatCountNonzeros;
346  r1matops->id=15;
347  r1matops->matname=datamatname;
348  return 0;
349 }
350 
dsdpdatamat_impl.h
Structure of function pointers that each SDP data matrix type (sparse, dense, constant,...
DSDPGetR1PMat
int DSDPGetR1PMat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, struct DSDPDataMat_Ops **mops, void **mmat)
Create a rank one matrix usuable by DSDP in packed symmetric format.
Definition: rmmat.c:77
DSDPDataMatOpsInitialize
int DSDPDataMatOpsInitialize(struct DSDPDataMat_Ops *dops)
Initialize the table of function pointers for SDP Data matrices.
Definition: dsdpdatamat.c:47
dsdpsys.h
Error handling, printing, and profiling.
DSDPDataMat_Ops
Table of function pointers that operate on the data matrix.
Definition: dsdpdatamat_impl.h:15
DSDPGetR1UMat
int DSDPGetR1UMat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, struct DSDPDataMat_Ops **mops, void **mmat)
Create a rank one matrix usuable by DSDP in full symmetric format.
Definition: rmmat.c:101