1 #include "dsdpdatamat_impl.h"
2 #include "dsdpsys.h"
3 /*! \file vech.c
4 \brief DSDPDataMat for sparse matrices in upper packed symmetric format
5 */
6
7 typedef struct {
8 int neigs;
9 double *eigval;
10 double *an;
11 int *cols,*nnz;
12 } Eigen;
13
14 /*!
15 struct vechmat ;
16
17 \brief A sparse symmetric data matrix is packed format.
18 */
19 typedef struct {
20 int nnzeros;
21 const int *ind;
22 const double *val;
23 int ishift;
24 double alpha;
25
26 Eigen *Eig;
27 int factored;
28 int owndata;
29 int n;
30 } vechmat;
31
32 #define GETI(a) (int)(sqrt(2*(a)+0.25)-0.5)
33 #define GETJ(b,c) ((b)-((c)*((c)+1))/2)
34
getij(int k,int * i,int * j)35 static void getij(int k, int *i,int *j){
36 *i=GETI(k);
37 *j=GETJ(k,*i);
38 return;
39 }
40 /*
41 static void geti2(int k, int *i,int *j){
42 int r,c;
43 double rr=sqrt(2*k+0.25)-0.5;
44 r=(int)rr;
45 c=k-(r*(r+1))/2;
46 *i=r;*j=c;
47 return;
48 }
49 */
50 #undef __FUNCT__
51 #define __FUNCT__ "CreateVechMatWData"
CreateVechMatWdata(int n,int ishift,double alpha,const int * ind,const double * vals,int nnz,vechmat ** A)52 static int CreateVechMatWdata(int n, int ishift, double alpha,const int *ind, const double *vals, int nnz, vechmat **A){
53 int info;
54 vechmat* V;
55 DSDPCALLOC1(&V,vechmat,&info);DSDPCHKERR(info);
56 V->n=n; V->ishift=ishift, V->ind=ind; V->val=vals;V->nnzeros=nnz;
57 V->alpha=alpha;
58 V->owndata=0;
59 V->Eig=0;
60 *A=V;
61 return 0;
62 }
63
64
VechMatAddMultiple(void * AA,double scl,double r[],int nn,int n)65 static int VechMatAddMultiple(void* AA, double scl, double r[], int nn, int n){
66 vechmat* A=(vechmat*)AA;
67 int k,nnz=A->nnzeros;
68 const int* ind=A->ind;
69 const double *val=A->val;
70 double *rr=r-A->ishift;
71 scl*=A->alpha;
72 for (k=0; k<nnz; ++k,++ind,++val){
73 *(rr+((*ind))) +=scl*(*val);
74 }
75 return 0;
76 }
77
VechMatDot(void * AA,double x[],int nn,int n,double * v)78 static int VechMatDot(void* AA, double x[], int nn, int n, double *v){
79 vechmat* A=(vechmat*)AA;
80 int k,nnz=A->nnzeros;
81 const int *ind=A->ind;
82 double vv=0, *xx=x-A->ishift;
83 const double *val=A->val;
84 for (k=0;k<nnz;++k,++ind,++val){
85 vv+=(*val)*(*(xx+(*ind)));
86 }
87 *v=2*vv*A->alpha;
88 return 0;
89 }
90
91 static int EigMatVecVec(Eigen*, double[], int, double*);
92 static int VechMatGetRank(void*,int*,int);
93
VechMatVecVec(void * AA,double x[],int n,double * v)94 static int VechMatVecVec(void* AA, double x[], int n, double *v){
95 vechmat* A=(vechmat*)AA;
96 int info,rank=n,i=0,j,k,kk;
97 const int *ind=A->ind,ishift=A->ishift;
98 double vv=0,dd;
99 const double *val=A->val,nnz=A->nnzeros;
100
101 if (A->factored==3){
102 info=VechMatGetRank(AA,&rank,n);
103 if (nnz>3 && rank<nnz){
104 info=EigMatVecVec(A->Eig,x,n,&vv);
105 *v=vv*A->alpha;
106 return 0;
107 }
108 }
109
110 for (k=0; k<nnz; ++k,++ind,++val){
111 kk=*ind-ishift;
112 i=GETI(kk);
113 j=GETJ(kk,i);
114 dd=x[i]*x[j]*(*val);
115 vv+=2*dd;
116 if (i==j){ vv-=dd; }
117 }
118 *v=vv*A->alpha;
119
120 return 0;
121 }
122
123
VechMatGetRowNnz(void * AA,int trow,int nz[],int * nnzz,int nn)124 static int VechMatGetRowNnz(void* AA, int trow, int nz[], int *nnzz,int nn){
125 vechmat* A=(vechmat*)AA;
126 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros;
127 const int *ind=A->ind;
128 *nnzz=0;
129 for (k=0; k<nnz; ++k,++ind){
130 getij((*ind)-ishift,&i,&j);
131 if (i==trow){
132 nz[j]++;(*nnzz)++;
133 } else if (j==trow){
134 nz[i]++;(*nnzz)++;
135 }
136 }
137 return 0;
138 }
139
VechMatFNorm2(void * AA,int n,double * fnorm2)140 static int VechMatFNorm2(void* AA, int n, double *fnorm2){
141 vechmat* A=(vechmat*)AA;
142 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros;
143 const int *ind=A->ind;
144 double fn2=0;
145 const double *val=A->val;
146 for (k=0; k<nnz; ++k,++ind){
147 getij((*ind)-ishift,&i,&j);
148 if (i==j){
149 fn2+= val[k]*val[k];
150 } else {
151 fn2+= 2*val[k]*val[k];
152 }
153 }
154 *fnorm2=fn2*A->alpha*A->alpha;
155 return 0;
156 }
157
VechMatAddRowMultiple(void * AA,int trow,double scl,double r[],int m)158 static int VechMatAddRowMultiple(void* AA, int trow, double scl, double r[], int m){
159 vechmat* A=(vechmat*)AA;
160 int i=0,j,k,ishift=A->ishift,nnz=A->nnzeros;
161 const int *ind=A->ind;
162 const double *val=A->val;
163 scl*=A->alpha;
164 for (k=0; k<nnz; ++k,++ind){
165 getij((*ind)-ishift,&i,&j);
166 if (i==trow){
167 /* if (i==j){ r[j]+=scl*val[k];} */
168 r[j]+=scl*val[k];
169 } else if (j==trow){
170 r[i]+=scl*val[k];
171 }
172 }
173 return 0;
174 }
175
VechMatCountNonzeros(void * AA,int * nnz,int n)176 static int VechMatCountNonzeros(void* AA, int*nnz, int n){
177 vechmat* A=(vechmat*)AA;
178 *nnz=A->nnzeros;
179 return 0;
180 }
181
182 #undef __FUNCT__
183 #define __FUNCT__ "VechMatDestroy"
VechMatDestroy(void * AA)184 static int VechMatDestroy(void* AA){
185 vechmat* A=(vechmat*)AA;
186 int info;
187 if (A->owndata){
188 /* Never happens
189 if (A->ind){ DSDPFREE(&A->ind,&info);DSDPCHKERR(info);}
190 if (A->val){ DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
191 */
192 return 1;
193 }
194 if (A->Eig){
195 DSDPFREE(&A->Eig->eigval,&info);DSDPCHKERR(info);
196 DSDPFREE(&A->Eig->an,&info);DSDPCHKERR(info);
197 if (A->Eig->cols){DSDPFREE(&A->Eig->cols,&info);DSDPCHKERR(info);}
198 if (A->Eig->nnz){DSDPFREE(&A->Eig->nnz,&info);DSDPCHKERR(info);}
199 DSDPFREE(&A->Eig,&info);DSDPCHKERR(info);
200 }
201 DSDPFREE(&A,&info);DSDPCHKERR(info);
202 return 0;
203 }
204
205
206
207 #undef __FUNCT__
208 #define __FUNCT__ "DSDPCreateVechMatEigs"
CreateEigenLocker(Eigen ** EE,int iptr[],int neigs,int n)209 static int CreateEigenLocker(Eigen **EE,int iptr[], int neigs, int n){
210 int i,k,info;
211 Eigen *E;
212
213 for (k=0,i=0;i<neigs;i++) k+=iptr[i];
214 if (k>n*neigs/4){
215
216 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
217 DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info);
218 DSDPCALLOC2(&E->an,double,n*neigs,&info);DSDPCHKERR(info);
219 E->neigs=neigs;
220 E->cols=0;
221 E->nnz=0;
222
223 } else {
224
225 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
226 DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info);
227 DSDPCALLOC2(&E->nnz,int,neigs,&info);DSDPCHKERR(info);
228 DSDPCALLOC2(&E->an,double,k,&info);DSDPCHKERR(info);
229 DSDPCALLOC2(&E->cols,int,k,&info);DSDPCHKERR(info);
230 E->neigs=neigs;
231
232 if (neigs>0) E->nnz[0]=iptr[0];
233 for (i=1;i<neigs;i++){E->nnz[i]=E->nnz[i-1]+iptr[i];}
234 }
235 *EE=E;
236 return 0;
237 }
238
239
EigMatSetEig(Eigen * A,int row,double eigv,int idxn[],double v[],int nsub,int n)240 static int EigMatSetEig(Eigen* A,int row, double eigv, int idxn[], double v[], int nsub,int n){
241 int j,k,*cols=A->cols;
242 double *an=A->an;
243
244 A->eigval[row]=eigv;
245 if (cols){
246 k=0; if (row>0){ k=A->nnz[row-1];}
247 cols+=k; an+=k;
248 for (k=0,j=0; j<nsub; j++){
249 if (v[j]==0.0) continue;
250 cols[k]=idxn[j]; an[k]=v[j]; k++;
251 }
252 } else {
253 an+=n*row;
254 for (j=0; j<nsub; j++){
255 if (v[j]==0.0) continue;
256 an[idxn[j]]=v[j];
257 }
258 }
259 return 0;
260 }
261
262
EigMatGetEig(Eigen * A,int row,double * eigenvalue,double eigenvector[],int n,int spind[],int * nind)263 static int EigMatGetEig(Eigen* A,int row, double *eigenvalue, double eigenvector[], int n, int spind[], int *nind){
264 int i,*cols=A->cols,bb,ee;
265 double* an=A->an;
266 *eigenvalue=A->eigval[row];
267 *nind=0;
268 if (cols){
269 memset((void*)eigenvector,0,n*sizeof(double));
270 if (row==0){ bb=0;} else {bb=A->nnz[row-1];} ee=A->nnz[row];
271 for (i=bb;i<ee;i++){
272 eigenvector[cols[i]]=an[i];
273 spind[i-bb]=cols[i]; (*nind)++;
274 }
275 } else {
276 memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double));
277 for (i=0;i<n;i++)spind[i]=i;
278 *nind=n;
279 }
280 return 0;
281 }
282
EigMatVecVec(Eigen * A,double v[],int n,double * vv)283 static int EigMatVecVec(Eigen* A, double v[], int n, double *vv){
284 int i,rank,*cols=A->cols,neigs=A->neigs,*nnz=A->nnz,bb,ee;
285 double* an=A->an,*eigval=A->eigval,dd,ddd=0;
286
287 if (cols){
288 for (rank=0;rank<neigs;rank++){
289 if (rank==0){ bb=0;} else {bb=nnz[rank-1];} ee=nnz[rank];
290 for (dd=0,i=bb;i<ee;i++){
291 dd+=an[i]*v[cols[i]];
292 }
293 ddd+=dd*dd*eigval[rank];
294 }
295 } else {
296 for (rank=0;rank<neigs;rank++){
297 for (dd=0,i=0;i<n;i++){
298 dd+=an[i]*v[i];
299 }
300 an+=n;
301 ddd+=dd*dd*eigval[rank];
302 }
303 }
304 *vv=ddd;
305 return 0;
306 }
307
308
309 static int VechMatComputeEigs(vechmat*,double[],int,double[],int,double[],int,int[],int,double[],int,double[],int);
310
VechMatFactor(void * AA,double dmatp[],int nn0,double dwork[],int n,double ddwork[],int n1,int iptr[],int n2)311 static int VechMatFactor(void*AA, double dmatp[], int nn0, double dwork[], int n, double ddwork[], int n1, int iptr[], int n2){
312 vechmat* A=(vechmat*)AA;
313 int i,j,k,ishift=A->ishift,isdiag,nonzeros=A->nnzeros,info;
314 int nn1=0,nn2=0;
315 double *ss1=0,*ss2=0;
316 const int *ind=A->ind;
317
318 if (A->factored) return 0;
319 memset((void*)iptr,0,3*n*sizeof(int));
320 /* Find number of nonzeros in each row */
321 for (isdiag=1,k=0; k<nonzeros; k++){
322 getij(ind[k]-ishift,&i,&j);
323 iptr[i]++;
324 if (i!=j) {isdiag=0;iptr[j]++;}
325 }
326
327 if (isdiag){ A->factored=1; return 0;}
328 /* Find most nonzeros per row */
329 for (j=0,i=0; i<n; i++){ if (iptr[i]>j) j=iptr[i]; }
330 if (j<2){ A->factored=2; return 0; }
331 info=VechMatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2,ss1,nn1,ss2,nn2);DSDPCHKERR(info);
332 A->factored=3;
333 return 0;
334 }
335
VechMatGetRank(void * AA,int * rank,int n)336 static int VechMatGetRank(void *AA,int *rank,int n){
337 vechmat* A=(vechmat*)AA;
338 switch (A->factored){
339 case 1:
340 *rank=A->nnzeros;
341 break;
342 case 2:
343 *rank=2*A->nnzeros;
344 break;
345 case 3:
346 *rank=A->Eig->neigs;
347 break;
348 default:
349 DSDPSETERR(1,"Vech Matrix not factored yet\n");
350 }
351 return 0;
352 }
353
VechMatGetEig(void * AA,int rank,double * eigenvalue,double vv[],int n,int indx[],int * nind)354 static int VechMatGetEig(void* AA, int rank, double *eigenvalue, double vv[], int n, int indx[], int *nind){
355 vechmat* A=(vechmat*)AA;
356 const double *val=A->val,tt=sqrt(0.5);
357 int info,i,j,k,ishift=A->ishift;
358 const int *ind=A->ind;
359
360 *nind=0;
361 switch (A->factored){
362 case 1:
363 memset(vv,0,n*sizeof(double));
364 getij(ind[rank]-ishift,&i,&j);
365 vv[i]=1.0;
366 *eigenvalue=val[rank]*A->alpha;
367 *nind=1;
368 indx[0]=i;
369 break;
370 case 2:
371 memset(vv,0,n*sizeof(double));
372 k=rank/2;
373 getij(ind[k]-ishift,&i,&j);
374 if (i==j){
375 if (k*2==rank){
376 vv[i]=1.0; *eigenvalue=val[k]*A->alpha;
377 *nind=1;
378 indx[0]=i;
379 } else {
380 *eigenvalue=0;
381 }
382 } else {
383 if (k*2==rank){
384 vv[i]=tt; vv[j]=tt; *eigenvalue=val[k]*A->alpha;
385 *nind=2;
386 indx[0]=i; indx[1]=j;
387 } else {
388 vv[i]=-tt; vv[j]=tt; *eigenvalue=-val[k]*A->alpha;
389 *nind=2;
390 indx[0]=i; indx[1]=j;
391 }
392 }
393 break;
394 case 3:
395 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n,indx,nind);DSDPCHKERR(info);
396 *eigenvalue=*eigenvalue*A->alpha;
397 break;
398 default:
399 DSDPSETERR(1,"Vech Matrix not factored yet\n");
400 }
401
402 return 0;
403 }
404
VechMatView(void * AA)405 static int VechMatView(void* AA){
406 vechmat* A=(vechmat*)AA;
407 int info,i=0,j,k,rank=0,ishift=A->ishift,nnz=A->nnzeros;
408 const int *ind=A->ind;
409 const double *val=A->val;
410 for (k=0; k<nnz; k++){
411 getij(ind[k]-ishift,&i,&j);
412 printf("Row: %d, Column: %d, Value: %10.8f \n",i,j,A->alpha*val[k]);
413 }
414 if (A->factored>0){
415 info=VechMatGetRank(AA,&rank,A->n);DSDPCHKERR(info);
416 printf("Detected Rank: %d\n",rank);
417 }
418 return 0;
419 }
420
421
422 static struct DSDPDataMat_Ops vechmatops;
423 static const char *datamatname="STANDARD VECH MATRIX";
424
VechMatOpsInitialize(struct DSDPDataMat_Ops * sops)425 static int VechMatOpsInitialize(struct DSDPDataMat_Ops *sops){
426 int info;
427 if (sops==NULL) return 0;
428 info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info);
429 sops->matvecvec=VechMatVecVec;
430 sops->matdot=VechMatDot;
431 sops->matfnorm2=VechMatFNorm2;
432 sops->mataddrowmultiple=VechMatAddRowMultiple;
433 sops->mataddallmultiple=VechMatAddMultiple;
434 sops->matview=VechMatView;
435 sops->matdestroy=VechMatDestroy;
436 sops->matfactor2=VechMatFactor;
437 sops->matgetrank=VechMatGetRank;
438 sops->matgeteig=VechMatGetEig;
439 sops->matrownz=VechMatGetRowNnz;
440 sops->matnnz=VechMatCountNonzeros;
441 sops->id=3;
442 sops->matname=datamatname;
443 return 0;
444 }
445
446 /*!
447 \fn int DSDPGetVechMat(int n,int ishift,double alpha, const int ind[], const double val[],int nnz, struct DSDPDataMat_Ops**sops, void**smat)
448 \brief Given data in packed symmetric format, create a sparse matrix usuable by DSDP.
449 \param n number of rows and columns of the matrix
450 \param ishift the index of the first element in the matrix (usually 0)
451 \param alpha the multiple of these matrix.
452 \param ind array of indices for matrix.
453 \param val array of matrix values.
454 \param nnz size of arrays.
455 \param sops address of a pointer to a table of function pointers
456 \param smat address of a pointer to an opaque data type.
457 */
458 #undef __FUNCT__
459 #define __FUNCT__ "DSDPGetVechMat"
DSDPGetVechMat(int n,int ishift,double alpha,const int ind[],const double val[],int nnz,struct DSDPDataMat_Ops ** sops,void ** smat)460 int DSDPGetVechMat(int n,int ishift,double alpha, const int ind[], const double val[],int nnz, struct DSDPDataMat_Ops**sops, void**smat){
461 int info,i,j,k,itmp,nn=n*(n+1)/2;
462 double dtmp;
463 vechmat* AA;
464 DSDPFunctionBegin;
465 for (k=0;k<nnz;++k){
466 itmp=ind[k]-ishift;
467 if (itmp>=nn){
468 getij(itmp,&i,&j);
469 /*
470 DSDPSETERR(2,"Illegal index value: Element %d in array has row %d (>0) or column %d (>0) is greater than %d. \n",k+1,i+1,j+1,n);
471 */
472 DSDPSETERR3(2,"Illegal index value: Element %d in array has index %d greater than or equal to %d. \n",k,itmp,nn);
473 } else if (itmp<0){
474 DSDPSETERR1(2,"Illegal index value: %d. Must be >= 0\n",itmp);
475 }
476 }
477 for (k=0;k<nnz;++k) dtmp=val[k];
478 info=CreateVechMatWdata(n,ishift,alpha,ind,val,nnz,&AA); DSDPCHKERR(info);
479 AA->factored=0;
480 AA->Eig=0;
481 info=VechMatOpsInitialize(&vechmatops); DSDPCHKERR(info);
482 if (sops){*sops=&vechmatops;}
483 if (smat){*smat=(void*)AA;}
484 DSDPFunctionReturn(0);
485 }
486
487
488 #undef __FUNCT__
489 #define __FUNCT__ "VechMatComputeEigs"
VechMatComputeEigs(vechmat * AA,double DD[],int nn0,double W[],int n,double WORK[],int n1,int iiptr[],int n2,double ss1[],int nn1,double ss2[],int nn2)490 static int VechMatComputeEigs(vechmat* AA,double DD[], int nn0, double W[], int n, double WORK[], int n1, int iiptr[], int n2, double ss1[],int nn1, double ss2[], int nn2){
491
492 int i,j,k,nsub,neigs,info,*iptr,*perm,*invp;
493 long int *i2darray=(long int*)DD;
494 int ownarray1=0,ownarray2=0,ownarray3=0;
495 int ishift=AA->ishift,nonzeros=AA->nnzeros;
496 const int *ind=AA->ind;
497 const double *val=AA->val;
498 double *dmatarray=ss1,*dworkarray=ss2,maxeig,eps=1.0e-12,eps2=1.0e-12;
499
500 iptr=iiptr; perm=iptr+n; invp=perm+n;
501 /* These operations were done before calling this routine * /
502 / * Integer arrays corresponding to rows with nonzeros and inverse map * /
503 memset((void*)iiptr,0,3*n*sizeof(int));
504
505 / * Find number of nonzeros in each row * /
506 for (i=0,k=0; k<nonzeros; k++){
507 getij(ind[k],i,n,&i,&j);
508 iptr[i]++; iptr[j]++;
509 }
510 */
511 /* Number of rows with a nonzero. Order the rows with nonzeros. */
512 for (nsub=0,i=0; i<n; i++){
513 if (iptr[i]>0){ invp[nsub]=i; perm[i]=nsub; nsub++;}
514 }
515
516 /* create a dense array in which to put numbers */
517 if (nsub*nsub>nn1){
518 DSDPCALLOC2(&dmatarray,double,(nsub*nsub),&info); DSDPCHKERR(info);
519 ownarray1=1;
520 }
521 memset((void*)dmatarray,0,nsub*nsub*sizeof(double));
522 if (nsub*nsub>nn2){
523 DSDPCALLOC2(&dworkarray,double,(nsub*nsub),&info); DSDPCHKERR(info);
524 ownarray2=1;
525 }
526
527 if (nsub*nsub*sizeof(long int)>nn0*sizeof(double)){
528 DSDPCALLOC2(&i2darray,long int,(nsub*nsub),&info); DSDPCHKERR(info);
529 ownarray3=1;
530 }
531
532 for (i=0,k=0; k<nonzeros; k++){
533 getij(ind[k]-ishift,&i,&j);
534 dmatarray[perm[i]*nsub+perm[j]] += val[k];
535 if (i!=j){
536 dmatarray[perm[j]*nsub+perm[i]] += val[k];
537 }
538 }
539 /* Call LAPACK to compute the eigenvalues */
540 info=DSDPGetEigs(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
541 W,n,WORK,n1,iiptr+3*n,n2-3*n);
542 if (info){
543 memset((void*)dmatarray,0,nsub*nsub*sizeof(double));
544 for (i=0,k=0; k<nonzeros; k++){
545 getij(ind[k]-ishift,&i,&j);
546 dmatarray[perm[i]*nsub+perm[j]] += val[k];
547 if (i!=j){dmatarray[perm[j]*nsub+perm[i]] += val[k];}
548 }
549 info=DSDPGetEigs2(dmatarray,nsub,dworkarray,nsub*nsub,i2darray,nsub*nsub,
550 W,nsub,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
551 }
552 for (maxeig=0,i=0;i<nsub;i++){
553 if (fabs(W[i])>maxeig){ maxeig=fabs(W[i]); }
554 }
555 memset((void*)iptr,0,nsub*sizeof(int));
556 /* Compute sparsity pattern for eigenvalue and eigenvector structures */
557 /* Count the nonzero eigenvalues */
558 for (neigs=0,k=0; k<nsub; k++){
559 if (fabs(W[k]) /* /maxeig */ > eps){
560 for (j=0;j<nsub;j++){
561 if (fabs(dmatarray[nsub*k+j]) >= eps2){iptr[neigs]++;
562 } else { dmatarray[nsub*k+j]=0.0;}
563 }
564 neigs++;
565 /*
566 } else if (fabs(W[k])>1.0e-100){
567 printf("SKIPPING EIGENVALUE: %4.4e, max is : %4.4e\n",W[k],maxeig);
568 */
569 }
570 }
571
572 info=CreateEigenLocker(&AA->Eig,iptr,neigs,n);DSDPCHKERR(info);
573 DSDPLogInfo(0,49," Data Mat has %d eigenvalues: \n",neigs);
574 /* Copy into structure */
575 for (neigs=0,i=0; i<nsub; i++){
576 if (fabs(W[i]) > eps){
577 info=EigMatSetEig(AA->Eig,neigs,W[i],invp,dmatarray+nsub*i,nsub,n);DSDPCHKERR(info);
578 neigs++;
579 }
580 }
581
582 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
583 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
584 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
585 return 0;
586 }
587
588