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