DSDP
rmmat.c
Go to the documentation of this file.
1#include "dsdpdatamat_impl.h"
2#include "dsdpsys.h"
8typedef 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
18static int R1MatDestroy(void*);
19static int R1MatView(void*);
20static int R1MatVecVec(void*, double[], int, double *);
21static int R1MatDotP(void*, double[],int,int,double *);
22static int R1MatDotU(void*, double[],int,int,double *);
23static int R1MatGetRank(void*, int*, int);
24static int R1MatFactor(void*);
25static int R1MatGetEig(void*, int, double*, double[], int,int[],int*);
26static int R1MatRowNnz(void*, int, int[], int*, int);
27static int R1MatAddRowMultiple(void*, int, double, double[], int);
28static int R1MatAddMultipleP(void*, double, double[], int,int);
29static int R1MatAddMultipleU(void*, double, double[], int,int);
30
31static struct DSDPDataMat_Ops r1matopsP;
32static struct DSDPDataMat_Ops r1matopsU;
33static int R1MatOpsInitializeP(struct DSDPDataMat_Ops*);
34static int R1MatOpsInitializeU(struct DSDPDataMat_Ops*);
35
36
37#undef __FUNCT__
38#define __FUNCT__ "DSDPGetR1Mat"
39int 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"
77int 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"
101int 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
110static 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
132static 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
155static 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
169static 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}
189static 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
210static 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
228static int R1MatFactor(void*A){
229 return 0;
230}
231
232
233static int R1MatGetRank(void *A, int*rank, int n){
234 *rank=1;
235 return 0;
236}
237
238static 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
257static 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
274static 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
286static 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
293static 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
304static int R1MatDestroy(void* A){
305 if (A) free(A);
306 return 0;
307}
308
309static const char *datamatname="RANK 1 Outer Product";
310static 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}
330static 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 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,...
Error handling, printing, and profiling.
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
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
Table of function pointers that operate on the data matrix.