00001 #include "dsdpsys.h"
00002 #include "dsdpvec.h"
00003 #include "dsdplapack.h"
00004 #include "dsdpdatamat_impl.h"
00005
00011 typedef enum {
00012 Init=0,
00013 Assemble=1,
00014 Factored=2,
00015 Inverted=3,
00016 ISymmetric=4
00017 } MatStatus;
00018
00019 typedef struct{
00020 char UPLO;
00021 int LDA;
00022 double *val,*v2;
00023 double *sscale;
00024 double *workn;
00025 int scaleit;
00026 int n;
00027 int owndata;
00028 MatStatus status;
00029 } dtrumat;
00030
00031 static int DTRUMatView(void*);
00032
00033
00034 #undef __FUNCT__
00035 #define __FUNCT__ "DSDPLAPACKROUTINE"
00036 static void dtruscalevec(double alpha, double v1[], double v2[], double v3[], int n){
00037 int i;
00038 for (i=0;i<n;i++){
00039 v3[i] = (alpha*v1[i]*v2[i]);
00040 }
00041 return;
00042 }
00043
00044 static void dsydadd(double x[], double s[], double y[],int n){
00045 int i;
00046 for (i=0;i<n;i++){
00047 y[i] += x[i]*(1/(s[i]*s[i])-2);
00048 }
00049 return;
00050 }
00051
00052
00053 static void dtruscalemat(double vv[], double ss[], int n,int LDA){
00054 int i;
00055 for (i=0;i<n;i++,vv+=LDA){
00056 dtruscalevec(ss[i],vv,ss,vv,i+1);
00057 }
00058 return;
00059 }
00060
00061 static void DTRUMatOwnData(void* A, int owndata){
00062 dtrumat* AA=(dtrumat*)A;
00063 AA->owndata=owndata;
00064 return;
00065 }
00066
00067 static int SUMatGetLDA(int n){
00068 int nlda=n;
00069 if (n>8 && nlda%2==1){ nlda++;}
00070 if (n>100){
00071 while (nlda%8!=0){ nlda++;}
00072 }
00073
00074
00075
00076 return (nlda);
00077 }
00078
00079 static int DTRUMatCreateWData(int n,int LDA,double nz[], int nnz, dtrumat**M){
00080 int i,info;
00081 dtrumat*M23;
00082 if (nnz<n*n){DSDPSETERR1(2,"Array must have length of : %d \n",n*n);}
00083 DSDPCALLOC1(&M23,dtrumat,&info);DSDPCHKERR(info);
00084 DSDPCALLOC2(&M23->sscale,double,n,&info);DSDPCHKERR(info);
00085 DSDPCALLOC2(&M23->workn,double,n,&info);DSDPCHKERR(info);
00086 M23->owndata=0; M23->val=nz; M23->n=n; M23->UPLO='U';M23->LDA=n;
00087 M23->status=Init;
00088 for (i=0;i<n;i++)M23->sscale[i]=1.0;
00089 M23->scaleit=1;
00090 M23->LDA=LDA;
00091 if (n<=0){M23->LDA=1;}
00092 *M=M23;
00093 return 0;
00094 }
00095
00096
00097
00098 #undef __FUNCT__
00099 #define __FUNCT__ "DSDPGetEigs"
00100 int DSDPGetEigs(double A[],int n, double AA[], int nn0, long int IA[], int nn1,
00101 double W[],int n2,
00102 double WORK[],int nd, int LLWORK[], int ni){
00103 ffinteger N=n,LDA=n,LDZ=n,LWORK=nd,INFO=0;
00104 char UPLO='U',JOBZ='V',RANGE='A';
00105
00106
00107 LDA=DSDPMax(1,n);
00108 LDZ=DSDPMax(1,n);
00109 if ( n2/2.5 > n || (ni<10*n+1) || (nd<26*n+1) || (nn0 < n*LDA) || (nn1<LDZ*n) ){
00110
00111
00112
00113
00114 dsyev(&JOBZ,&UPLO,&N,A,&LDA,W,WORK,&LWORK,&INFO);
00115 } else {
00116
00117 int i;
00118 ffinteger M,IL=1,IU=n,*ISUPPZ=(ffinteger*)IA;
00119 ffinteger *IWORK=(ffinteger*)LLWORK,LIWORK=(ffinteger)ni;
00120 double *Z=AA,VL=-1e10,VU=1e10,ABSTOL=0;
00121
00122 if (0==1){
00123 dtrumat*M;
00124 DTRUMatCreateWData(n,n,A,n*n,&M);
00125 DTRUMatView((void*)M);
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 dsyevr(&JOBZ,&RANGE,&UPLO,&N,A,&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,Z,&LDZ,ISUPPZ,WORK,&LWORK,IWORK,&LIWORK,&INFO);
00142 for (i=0;i<N*N;i++){A[i]=Z[i];}
00143
00144 }
00145 return INFO;
00146 }
00147
00148 #undef __FUNCT__
00149 #define __FUNCT__ "DSDPGetEigsSTEP"
00150 int DSDPGetEigsSTEP(double A[],int n, double AA[], int nn0, long int IA[], int nn1,
00151 double W[],int n2,
00152 double WORK[],int nd, int LLWORK[], int ni){
00153 int info;
00154 info=DSDPGetEigs(A,n,AA,nn0,IA,nn1,W,n2,WORK,nd,LLWORK,ni);
00155 return info;
00156 }
00157
00158 #undef __FUNCT__
00159 #define __FUNCT__ "DSDPGetEigs2"
00160 int DSDPGetEigs2(double A[],int n, double AA[], int nn0, long int IA[], int nn1,
00161 double W[],int n2,
00162 double WORK[],int nd, int LLWORK[], int ni){
00163 ffinteger N=n,LDA=n,LDZ=n,LWORK=nd,INFO=0;
00164 char UPLO='U',JOBZ='V';
00165
00166 LDA=DSDPMax(1,n);
00167 LDZ=DSDPMax(1,n);
00168 if (0==1){
00169 dtrumat*M;
00170 DTRUMatCreateWData(n,n,A,n*n,&M);
00171 DTRUMatView((void*)M);
00172 }
00173 dsyev(&JOBZ,&UPLO,&N,A,&LDA,W,WORK,&LWORK,&INFO);
00174 return INFO;
00175 }
00176
00177
00178 #undef __FUNCT__
00179 #define __FUNCT__ "DSDP FULL SYMMETRIC U LAPACK ROUTINE"
00180
00181 static int DTRUMatMult(void* AA, double x[], double y[], int n){
00182 dtrumat* A=(dtrumat*) AA;
00183 ffinteger N=n,LDA=A->LDA,INCX=1,INCY=1;
00184 double BETA=0.0,ALPHA=1.0;
00185 double *AP=A->val,*Y=y,*X=x;
00186 char UPLO=A->UPLO,TRANS='N';
00187
00188 if (A->n != n) return 1;
00189 if (x==0 && n>0) return 3;
00190 if (0==1){
00191 dgemv(&TRANS,&N,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
00192 } else {
00193 dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
00194 }
00195 return 0;
00196 }
00197
00198
00199 static int DTRUMatMultR(void* AA, double x[], double y[], int n){
00200 dtrumat* A=(dtrumat*) AA;
00201 ffinteger N=n,LDA=A->LDA,INCX=1,INCY=1;
00202 double ALPHA=1.0;
00203 double *AP=A->val,*Y=y,*X=x,*ss=A->sscale,*W=A->workn;
00204 char UPLO=A->UPLO,TRANS='N',DIAG='U';
00205
00206 UPLO='L';
00207 if (A->n != n) return 1;
00208 if (x==0 && n>0) return 3;
00209
00210
00211 memset(y,0,n*sizeof(double));
00212
00213 memcpy(W,X,n*sizeof(double));
00214 TRANS='N'; UPLO='L';
00215 dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,W,&INCY);
00216 daxpy(&N,&ALPHA,W,&INCX,Y,&INCY);
00217
00218 memcpy(W,X,n*sizeof(double));
00219 TRANS='T'; UPLO='L';
00220 dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,W,&INCY);
00221 daxpy(&N,&ALPHA,W,&INCX,Y,&INCY);
00222
00223 dsydadd(x,ss,y,n);
00224 return 0;
00225 }
00226
00227
00228 static void DTRUMatScale(void*AA){
00229 dtrumat* A=(dtrumat*) AA;
00230 int i,N=A->n,LDA=A->LDA;
00231 double *ss=A->sscale,*v=A->val;
00232 for (i=0;i<N;i++){ ss[i]=*v;v+=(LDA+1);}
00233 for (i=0;i<N;i++){ if (ss[i]!=0.0){ss[i]=1.0/sqrt(fabs(ss[i]));} else {ss[i]=1.0; }}
00234 dtruscalemat(A->val,ss,N,LDA);
00235 }
00236
00237 static int DTRUMatCholeskyFactor(void* AA, int *flag){
00238 dtrumat* A=(dtrumat*) AA;
00239 ffinteger INFO,LDA=A->LDA,N=A->n;
00240 double *AP=A->val;
00241 char UPLO=A->UPLO;
00242
00243 if (A->scaleit){ DTRUMatScale(AA);}
00244 dpotrf(&UPLO, &N, AP, &LDA, &INFO );
00245 *flag=INFO;
00246 A->status=Factored;
00247 return 0;
00248 }
00249
00250
00251 int dtrsm2(char*,char*,char*,char*,ffinteger*,ffinteger*,double*,double*,ffinteger*,double*,ffinteger*);
00252
00253 static int DTRUMatSolve(void* AA, double b[], double x[],int n){
00254 dtrumat* A=(dtrumat*) AA;
00255 ffinteger ierr,INFO=0,NRHS=1,LDA=A->LDA,LDB=A->LDA,N=A->n;
00256 double *AP=A->val,*ss=A->sscale,ONE=1.0;
00257 char SIDE='L',UPLO=A->UPLO,TRANSA='N',DIAG='N';
00258
00259 dtruscalevec(1.0,ss,b,x,n);
00260
00261 if (0==1){
00262 dpotrs(&UPLO, &N, &NRHS, AP, &LDA, x, &LDB, &INFO );
00263 } else {
00264 TRANSA='T';
00265 ierr=dtrsm2(&SIDE,&UPLO,&TRANSA,&DIAG,&N, &NRHS, &ONE, AP, &LDA, x, &LDB );
00266 TRANSA='N';
00267 ierr=dtrsm2(&SIDE,&UPLO,&TRANSA,&DIAG,&N, &NRHS, &ONE, AP, &LDA, x, &LDB );
00268 }
00269
00270 dtruscalevec(1.0,ss,x,x,n);
00271 return INFO;
00272 }
00273
00274
00275 static int DTRUMatShiftDiagonal(void* AA, double shift){
00276 dtrumat* A=(dtrumat*) AA;
00277 int i,n=A->n, LDA=A->LDA;
00278 double *v=A->val;
00279 for (i=0; i<n; i++){
00280 v[i*LDA+i]+=shift;
00281 }
00282 return 0;
00283 }
00284
00285 static int DTRUMatAddRow(void* AA, int nrow, double dd, double row[], int n){
00286 dtrumat* A=(dtrumat*) AA;
00287 ffinteger ione=1,LDA=A->LDA,nn,INCX=1,INCY=A->LDA;
00288 double *vv=A->val;
00289
00290 nn=nrow;
00291 daxpy(&nn,&dd,row,&INCX,vv+nrow,&INCY);
00292 nn=nrow+1;
00293 daxpy(&nn,&dd,row,&ione,vv+nrow*LDA,&ione);
00294
00295 return 0;
00296 }
00297
00298 static int DTRUMatZero(void* AA){
00299 dtrumat* A=(dtrumat*) AA;
00300 int mn=A->n*(A->LDA);
00301 double *vv=A->val;
00302 memset((void*)vv,0,mn*sizeof(double));
00303 A->status=Assemble;
00304 return 0;
00305 }
00306
00307 static int DTRUMatGetSize(void *AA, int *n){
00308 dtrumat* A=(dtrumat*) AA;
00309 *n=A->n;
00310 return 0;
00311 }
00312
00313 static int DTRUMatDestroy(void* AA){
00314 int info;
00315 dtrumat* A=(dtrumat*) AA;
00316 if (A && A->owndata){DSDPFREE(&A->val,&info);DSDPCHKERR(info);}
00317 if (A && A->sscale) {DSDPFREE(&A->sscale,&info);DSDPCHKERR(info);}
00318 if (A && A->workn) {DSDPFREE(&A->workn,&info);DSDPCHKERR(info);}
00319 if (A) {DSDPFREE(&A,&info);DSDPCHKERR(info);}
00320 return 0;
00321 }
00322
00323 static int DTRUMatView(void* AA){
00324 dtrumat* M=(dtrumat*) AA;
00325 int i,j;
00326 double *val=M->val;
00327 ffinteger LDA=M->LDA;
00328 for (i=0; i<M->n; i++){
00329 for (j=0; j<=i; j++){
00330 printf(" %9.2e",val[i*LDA+j]);
00331 }
00332 for (j=i+1; j<M->LDA; j++){
00333 printf(" %9.1e",val[i*LDA+j]);
00334 }
00335 printf("\n");
00336 }
00337 return 0;
00338 }
00339
00340 static int DTRUMatView2(void* AA){
00341 dtrumat* M=(dtrumat*) AA;
00342 int i,j;
00343 double *val=M->v2;
00344 ffinteger LDA=M->LDA;
00345 for (i=0; i<M->n; i++){
00346 for (j=0; j<=i; j++){
00347 printf(" %9.2e",val[i*LDA+j]);
00348 }
00349 for (j=i+1; j<M->LDA; j++){
00350 printf(" %9.2e",val[i*LDA+j]);
00351 }
00352 printf("\n");
00353 }
00354 return 0;
00355 }
00356
00357
00358 #include "dsdpschurmat_impl.h"
00359 #include "dsdpdualmat_impl.h"
00360 #include "dsdpdatamat_impl.h"
00361 #include "dsdpxmat_impl.h"
00362 #include "dsdpdsmat_impl.h"
00363 #include "dsdpsys.h"
00364
00365 #undef __FUNCT__
00366 #define __FUNCT__ "Tassemble"
00367 static int DTRUMatAssemble(void*M){
00368 int info;
00369 double shift=1.0e-15;
00370 DSDPFunctionBegin;
00371 info= DTRUMatShiftDiagonal(M, shift); DSDPCHKERR(info);
00372 DSDPFunctionReturn(0);
00373 }
00374
00375 static int DTRUMatRowNonzeros(void*M, int row, double cols[], int *ncols,int nrows){
00376 int i;
00377 DSDPFunctionBegin;
00378 *ncols = row+1;
00379 for (i=0;i<=row;i++){
00380 cols[i]=1.0;
00381 }
00382 memset((void*)(cols+row+1),0,(nrows-row-1)*sizeof(int));
00383 DSDPFunctionReturn(0);
00384 }
00385
00386
00387 #undef __FUNCT__
00388 #define __FUNCT__ "TAddDiag"
00389 static int DTRUMatAddDiag(void*M, int row, double dd){
00390 int ii;
00391 dtrumat* ABA=(dtrumat*)M;
00392 ffinteger LDA=ABA->LDA;
00393 DSDPFunctionBegin;
00394 ii=row*LDA+row;
00395 ABA->val[ii] +=dd;
00396 DSDPFunctionReturn(0);
00397 }
00398 #undef __FUNCT__
00399 #define __FUNCT__ "TAddDiag2"
00400 static int DTRUMatAddDiag2(void*M, double diag[], int m){
00401 int row,ii;
00402 dtrumat* ABA=(dtrumat*)M;
00403 ffinteger LDA=ABA->LDA;
00404 DSDPFunctionBegin;
00405 for (row=0;row<m;row++){
00406 ii=row*LDA+row;
00407 ABA->val[ii] +=diag[row];
00408 }
00409 DSDPFunctionReturn(0);
00410 }
00411 static struct DSDPSchurMat_Ops dsdpmmatops;
00412 static const char* lapackname="DENSE,SYMMETRIC U STORAGE";
00413
00414 static int DSDPInitSchurOps(struct DSDPSchurMat_Ops* mops){
00415 int info;
00416 DSDPFunctionBegin;
00417 info=DSDPSchurMatOpsInitialize(mops);DSDPCHKERR(info);
00418 mops->matrownonzeros=DTRUMatRowNonzeros;
00419 mops->matscaledmultiply=DTRUMatMult;
00420 mops->matmultr=DTRUMatMultR;
00421 mops->mataddrow=DTRUMatAddRow;
00422 mops->mataddelement=DTRUMatAddDiag;
00423 mops->matadddiagonal=DTRUMatAddDiag2;
00424 mops->matshiftdiagonal=DTRUMatShiftDiagonal;
00425 mops->matassemble=DTRUMatAssemble;
00426 mops->matfactor=DTRUMatCholeskyFactor;
00427 mops->matsolve=DTRUMatSolve;
00428 mops->matdestroy=DTRUMatDestroy;
00429 mops->matzero=DTRUMatZero;
00430 mops->matview=DTRUMatView;
00431 mops->id=1;
00432 mops->matname=lapackname;
00433 DSDPFunctionReturn(0);
00434 }
00435
00436
00437 #undef __FUNCT__
00438 #define __FUNCT__ "DSDPGetLAPACKSUSchurOps"
00439 int DSDPGetLAPACKSUSchurOps(int n,struct DSDPSchurMat_Ops** sops, void** mdata){
00440 int info,nn,LDA;
00441 double *vv;
00442 dtrumat *AA;
00443 DSDPFunctionBegin;
00444
00445 LDA=SUMatGetLDA(n);
00446 nn=n*LDA;
00447 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
00448 info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
00449 AA->owndata=1;
00450 info=DSDPInitSchurOps(&dsdpmmatops); DSDPCHKERR(info);
00451 *sops=&dsdpmmatops;
00452 *mdata=(void*)AA;
00453 DSDPFunctionReturn(0);
00454 }
00455
00456
00457 static int DTRUMatCholeskyBackward(void* AA, double b[], double x[], int n){
00458 dtrumat* M=(dtrumat*) AA;
00459 ffinteger N=M->n,INCX=1,LDA=M->LDA;
00460 double *AP=M->val,*ss=M->sscale;
00461 char UPLO=M->UPLO,TRANS='N',DIAG='N';
00462 memcpy(x,b,N*sizeof(double));
00463 dtrsv(&UPLO,&TRANS,&DIAG, &N, AP, &LDA, x, &INCX );
00464 dtruscalevec(1.0,ss,x,x,n);
00465 return 0;
00466 }
00467
00468
00469 static int DTRUMatCholeskyForward(void* AA, double b[], double x[], int n){
00470 dtrumat* M=(dtrumat*) AA;
00471 ffinteger N=M->n,INCX=1,LDA=M->LDA;
00472 double *AP=M->val,*ss=M->sscale;
00473 char UPLO=M->UPLO,TRANS='T',DIAG='N';
00474 dtruscalevec(1.0,ss,b,x,n);
00475 dtrsv(&UPLO,&TRANS,&DIAG, &N, AP, &LDA, x, &INCX );
00476 return 0;
00477 }
00478
00479 static int DTRUMatLogDet(void* AA, double *dd){
00480 dtrumat* A=(dtrumat*) AA;
00481 int i,n=A->n,LDA=A->LDA;
00482 double d=0,*v=A->val,*ss=A->sscale;
00483 for (i=0; i<n; i++){
00484 if (*v<=0) return 1;
00485 d+=2*log(*v/ss[i]);
00486 v+=LDA+1;
00487 }
00488 *dd=d;
00489 return 0;
00490 }
00491
00492
00493 static int DTRUMatCholeskyForwardMultiply(void* AA, double x[], double y[], int n){
00494 dtrumat* A=(dtrumat*) AA;
00495 ffinteger i,j,N=A->n,LDA=A->LDA;
00496 double *AP=A->val,*ss=A->sscale;
00497
00498
00499 if (x==0 && N>0) return 3;
00500
00501
00502
00503
00504 for (i=0;i<N;i++)y[i]=0;
00505 for (i=0;i<N;i++){
00506 for (j=0;j<=i;j++){
00507 y[i]+=AP[j]*x[j];
00508 }
00509 AP=AP+LDA;
00510 }
00511
00512 for (i=0;i<N;i++){ y[i]=y[i]/ss[i];}
00513 return 0;
00514 }
00515
00516 static int DTRUMatCholeskyBackwardMultiply(void* AA, double x[], double y[], int n){
00517 dtrumat* A=(dtrumat*) AA;
00518 ffinteger i,j,N=A->n,LDA=A->LDA;
00519 double *AP=A->val,*ss=A->sscale;
00520
00521
00522 if (x==0 && N>0) return 3;
00523
00524
00525
00526
00527 for (i=0;i<N;i++)y[i]=0;
00528 for (i=0;i<N;i++){
00529 for (j=0;j<=i;j++){
00530 y[j]+=AP[j]*x[i]/ss[i];
00531 }
00532 AP=AP+LDA;
00533 }
00534
00535 for (i=0;i<-N;i++){ y[i]=y[i]/ss[i];}
00536 return 0;
00537 }
00538
00539 static int DTRUMatInvert(void* AA){
00540 dtrumat* A=(dtrumat*) AA;
00541 ffinteger INFO,LDA=A->LDA,N=A->n,nn=LDA*N;
00542 double *v=A->val,*AP=A->v2,*ss=A->sscale;
00543 char UPLO=A->UPLO;
00544 memcpy((void*)AP,(void*)v,nn*sizeof(double));
00545 dpotri(&UPLO, &N, AP, &LDA, &INFO );
00546 if (INFO){
00547 INFO=DTRUMatShiftDiagonal(AA,1e-11); INFO=0;
00548 memcpy((void*)AP,(void*)v,nn*sizeof(double));
00549 dpotri(&UPLO, &N, AP, &LDA, &INFO );
00550 }
00551 if (A->scaleit){
00552 dtruscalemat(AP,ss,N,LDA);
00553 }
00554 A->status=Inverted;
00555 return INFO;
00556
00557 }
00558
00559 static void DTRUSymmetrize(dtrumat* A){
00560 int i,j,n=A->n,row,LDA=A->LDA;
00561 double *v2=A->v2,*r1=A->v2,*r2=A->v2+LDA,*c1;
00562 for (i=0;i<n/2;i++){
00563 row=2*i;
00564 r1=v2+LDA*(row);
00565 r2=v2+LDA*(row+1);
00566 c1=v2+LDA*(row+1);
00567 r1[row+1]=c1[row];
00568 c1+=LDA;
00569 for (j=row+2;j<n;j++){
00570 r1[j]=c1[row];
00571 r2[j]=c1[row+1];
00572 c1+=LDA;
00573 }
00574 r1+=LDA*2;
00575 r2+=LDA*2;
00576 }
00577
00578 for (row=2*n/2;row<n;row++){
00579 r1=v2+LDA*(row);
00580 c1=v2+LDA*(row+1);
00581 for (j=row+1;j<n;j++){
00582 r1[j]=c1[row];
00583 c1+=LDA;
00584 }
00585 r1+=LDA;
00586 }
00587 A->status=ISymmetric;
00588 return;
00589 }
00590
00591 static int DTRUMatInverseAdd(void* AA, double alpha, double y[], int nn, int n){
00592 dtrumat* A=(dtrumat*) AA;
00593 ffinteger i,LDA=A->LDA,N,ione=1;
00594 double *v2=A->v2;
00595 for (i=0;i<n;i++){
00596 N=i+1;
00597 daxpy(&N,&alpha,v2,&ione,y,&ione);
00598 v2+=LDA; y+=n;
00599 }
00600 return 0;
00601 }
00602
00603 static int DTRUMatInverseAddP(void* AA, double alpha, double y[], int nn, int n){
00604 dtrumat* A=(dtrumat*) AA;
00605 ffinteger N,LDA=A->LDA,i,ione=1;
00606 double *v2=A->v2;
00607 for (i=0;i<n;i++){
00608 N=i+1;
00609 daxpy(&N,&alpha,v2,&ione,y,&ione);
00610 v2+=LDA; y+=i+1;
00611 }
00612 return 0;
00613 }
00614
00615 static void daddrow(double *v, double alpha, int i, double row[], int n){
00616 double *s1;
00617 ffinteger j,nn=n,ione=1;
00618 if (alpha==0){return;}
00619 nn=i+1; s1=v+i*n;
00620 daxpy(&nn,&alpha,s1,&ione,row,&ione);
00621 s1=v+i*n+n;
00622 for (j=i+1;j<n;j++){ row[j]+=alpha*s1[i]; s1+=n; }
00623 return;
00624 }
00625
00626
00627
00628
00629 static int DTRUMatInverseMultiply(void* AA, int indx[], int nind, double x[], double y[],int n){
00630 dtrumat* A=(dtrumat*) AA;
00631 ffinteger nn=n,LDA=A->LDA,N=A->n,INCX=1,INCY=1;
00632 double *AP=A->v2,*s1=A->v2,*s2,*X=x,*Y=y,ALPHA=1.0,BETA=0.0;
00633 char UPLO=A->UPLO,TRANS='N';
00634 int i,ii,usefull=1;
00635
00636 if (usefull){
00637 if (A->status==Inverted){
00638 DTRUSymmetrize(A);
00639 }
00640 if (nind < n/4){
00641 memset((void*)y,0,n*sizeof(double));
00642 for (ii=0;ii<nind;ii++){
00643 i=indx[ii]; nn=n; ALPHA=x[i];s2=s1+i*LDA;
00644 daxpy(&nn,&ALPHA,s2,&INCY,y,&INCX);
00645 }
00646 } else{
00647 ALPHA=1.0;
00648 dgemv(&TRANS,&N,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
00649 }
00650
00651 } else {
00652 if (nind<n/4 ){
00653 memset((void*)y,0,n*sizeof(double));
00654 for (ii=0;ii<nind;ii++){
00655 i=indx[ii]; ALPHA=x[i];
00656 daddrow(s1,ALPHA,i,y,n);
00657 }
00658 } else {
00659 ALPHA=1.0;
00660 dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY);
00661 }
00662 }
00663 return 0;
00664 }
00665
00666 static int DTRUMatSetXMat(void*A, double v[], int nn, int n){
00667 dtrumat* ABA=(dtrumat*)A;
00668 int i,LDA=ABA->LDA;
00669 double *vv=ABA->val;
00670 if (vv!=v){
00671 for (i=0;i<n;i++){
00672 memcpy((void*)vv,(void*)v,(i+1)*sizeof(double));
00673 vv+=LDA; v+=n;
00674 }
00675 }
00676 ABA->status=Assemble;
00677 return 0;
00678 }
00679 static int DTRUMatSetXMatP(void*A, double v[], int nn, int n){
00680 dtrumat* ABA=(dtrumat*)A;
00681 int i,LDA=ABA->LDA;
00682 double *vv=ABA->val;
00683 if (vv!=v){
00684 for (i=0;i<n;i++){
00685 memcpy((void*)vv,(void*)v,(i+1)*sizeof(double));
00686 v+=(i+1); vv+=LDA;
00687 }
00688 }
00689 ABA->status=Assemble;
00690 return 0;
00691 }
00692
00693 static int DTRUMatFull(void*A, int*full){
00694 *full=1;
00695 return 0;
00696 }
00697
00698 static int DTRUMatGetArray(void*A,double **v,int *n){
00699 dtrumat* ABA=(dtrumat*)A;
00700 *n=ABA->n*ABA->LDA;
00701 *v=ABA->val;
00702 return 0;
00703 }
00704
00705 static struct DSDPDualMat_Ops sdmatops;
00706 static int SDualOpsInitialize(struct DSDPDualMat_Ops* sops){
00707 int info;
00708 if (sops==NULL) return 0;
00709 info=DSDPDualMatOpsInitialize(sops);DSDPCHKERR(info);
00710 sops->matseturmat=DTRUMatSetXMat;
00711 sops->matgetarray=DTRUMatGetArray;
00712 sops->matcholesky=DTRUMatCholeskyFactor;
00713 sops->matsolveforward=DTRUMatCholeskyForward;
00714 sops->matsolvebackward=DTRUMatCholeskyBackward;
00715 sops->matinvert=DTRUMatInvert;
00716 sops->matinverseadd=DTRUMatInverseAdd;
00717 sops->matinversemultiply=DTRUMatInverseMultiply;
00718 sops->matforwardmultiply=DTRUMatCholeskyForwardMultiply;
00719 sops->matbackwardmultiply=DTRUMatCholeskyBackwardMultiply;
00720 sops->matfull=DTRUMatFull;
00721 sops->matdestroy=DTRUMatDestroy;
00722 sops->matgetsize=DTRUMatGetSize;
00723 sops->matview=DTRUMatView;
00724 sops->matlogdet=DTRUMatLogDet;
00725 sops->matname=lapackname;
00726 sops->id=1;
00727 return 0;
00728 }
00729
00730
00731 #undef __FUNCT__
00732 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate"
00733 static int DSDPLAPACKSUDualMatCreate(int n,struct DSDPDualMat_Ops **sops, void**smat){
00734 dtrumat *AA;
00735 int info,nn,LDA=n;
00736 double *vv;
00737 DSDPFunctionBegin;
00738 LDA=SUMatGetLDA(n);
00739 nn=n*LDA;
00740 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
00741 info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
00742 AA->owndata=1;
00743 info=SDualOpsInitialize(&sdmatops);DSDPCHKERR(info);
00744 *sops=&sdmatops;
00745 *smat=(void*)AA;
00746 DSDPFunctionReturn(0);
00747 }
00748
00749
00750 static int switchptr(void *SD,void *SP){
00751 dtrumat *s1,*s2;
00752 s1=(dtrumat*)(SD);
00753 s2=(dtrumat*)(SP);
00754 s1->v2=s2->val;
00755 s2->v2=s1->val;
00756 return 0;
00757 }
00758
00759
00760 #undef __FUNCT__
00761 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate2"
00762 int DSDPLAPACKSUDualMatCreate2(int n,
00763 struct DSDPDualMat_Ops **sops1, void**smat1,
00764 struct DSDPDualMat_Ops **sops2, void**smat2){
00765 int info;
00766 DSDPFunctionBegin;
00767 info=DSDPLAPACKSUDualMatCreate(n,sops1,smat1);DSDPCHKERR(info);
00768 info=DSDPLAPACKSUDualMatCreate(n,sops2,smat2);DSDPCHKERR(info);
00769 info=switchptr(*smat1,*smat2);DSDPCHKERR(info);
00770 DSDPFunctionReturn(0);
00771 }
00772
00773 static struct DSDPDualMat_Ops sdmatopsp;
00774 static int SDualOpsInitializeP(struct DSDPDualMat_Ops* sops){
00775 int info;
00776 if (sops==NULL) return 0;
00777 info=DSDPDualMatOpsInitialize(sops);DSDPCHKERR(info);
00778 sops->matseturmat=DTRUMatSetXMatP;
00779 sops->matgetarray=DTRUMatGetArray;
00780 sops->matcholesky=DTRUMatCholeskyFactor;
00781 sops->matsolveforward=DTRUMatCholeskyForward;
00782 sops->matsolvebackward=DTRUMatCholeskyBackward;
00783 sops->matinvert=DTRUMatInvert;
00784 sops->matinverseadd=DTRUMatInverseAddP;
00785 sops->matinversemultiply=DTRUMatInverseMultiply;
00786 sops->matforwardmultiply=DTRUMatCholeskyForwardMultiply;
00787 sops->matbackwardmultiply=DTRUMatCholeskyBackwardMultiply;
00788 sops->matfull=DTRUMatFull;
00789 sops->matdestroy=DTRUMatDestroy;
00790 sops->matgetsize=DTRUMatGetSize;
00791 sops->matview=DTRUMatView;
00792 sops->matlogdet=DTRUMatLogDet;
00793 sops->matname=lapackname;
00794 sops->id=1;
00795 return 0;
00796 }
00797
00798 #undef __FUNCT__
00799 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate"
00800 static int DSDPLAPACKSUDualMatCreateP(int n, struct DSDPDualMat_Ops **sops, void**smat){
00801 dtrumat *AA;
00802 int info,nn,LDA;
00803 double *vv;
00804 DSDPFunctionBegin;
00805 LDA=SUMatGetLDA(n);
00806 nn=LDA*n;
00807 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
00808 info=DTRUMatCreateWData(n,LDA,vv,nn,&AA); DSDPCHKERR(info);
00809 AA->owndata=1;
00810 info=SDualOpsInitializeP(&sdmatopsp);DSDPCHKERR(info);
00811 *sops=&sdmatopsp;
00812 *smat=(void*)AA;
00813 DSDPFunctionReturn(0);
00814 }
00815
00816
00817 #undef __FUNCT__
00818 #define __FUNCT__ "DSDPLAPACKSUDualMatCreate2P"
00819 int DSDPLAPACKSUDualMatCreate2P(int n,
00820 struct DSDPDualMat_Ops* *sops1, void**smat1,
00821 struct DSDPDualMat_Ops* *sops2, void**smat2){
00822 int info;
00823 DSDPFunctionBegin;
00824 info=DSDPLAPACKSUDualMatCreateP(n,sops1,smat1);
00825 info=DSDPLAPACKSUDualMatCreateP(n,sops2,smat2);
00826 info=switchptr(*smat1,*smat2);
00827 DSDPFunctionReturn(0);
00828 }
00829
00830 static int DTRUMatScaleDiagonal(void* AA, double dd){
00831 dtrumat* A=(dtrumat*) AA;
00832 ffinteger LDA=A->LDA;
00833 int i,n=A->n;
00834 double *v=A->val;
00835 for (i=0; i<n; i++){
00836 *v*=dd;
00837 v+=LDA+1;
00838 }
00839 return 0;
00840 }
00841
00842 static int DTRUMatOuterProduct(void* AA, double alpha, double x[], int n){
00843 dtrumat* A=(dtrumat*) AA;
00844 ffinteger ione=1,N=n,LDA=A->LDA;
00845 double *v=A->val;
00846 char UPLO=A->UPLO;
00847 dsyr(&UPLO,&N,&alpha,x,&ione,v,&LDA);
00848 return 0;
00849 }
00850
00851 static int DenseSymPSDNormF2(void* AA, int n, double *dddot){
00852 dtrumat* A=(dtrumat*) AA;
00853 ffinteger ione=1,nn=A->n*A->n;
00854 double dd,tt=sqrt(0.5),*val=A->val;
00855 int info;
00856 info=DTRUMatScaleDiagonal(AA,tt);
00857 dd=dnrm2(&nn,val,&ione);
00858 info=DTRUMatScaleDiagonal(AA,1.0/tt);
00859 *dddot=dd*dd*2;
00860 return 0;
00861 }
00862
00863 static int DTRUMatGetDenseArray(void* A, double *v[], int*n){
00864 dtrumat* ABA=(dtrumat*)A;
00865 *v=ABA->val;
00866 *n=ABA->n*ABA->LDA;
00867 return 0;
00868 }
00869
00870 static int DTRUMatRestoreDenseArray(void* A, double *v[], int *n){
00871 *v=0;*n=0;
00872 return 0;
00873 }
00874
00875 static int DDenseSetXMat(void*A, double v[], int nn, int n){
00876 dtrumat* ABA=(dtrumat*)A;
00877 int i,LDA=ABA->LDA;
00878 double *vv=ABA->val;
00879 if (v!=vv){
00880 for (i=0;i<n;i++){
00881 memcpy((void*)vv,(void*)v,nn*sizeof(double));
00882 v+=n;vv+=LDA;
00883 }
00884 }
00885 ABA->status=Assemble;
00886 return 0;
00887 }
00888
00889 static int DDenseVecVec(void* A, double x[], int n, double *v){
00890 dtrumat* ABA=(dtrumat*)A;
00891 int i,j,k=0,LDA=ABA->LDA;
00892 double dd=0,*val=ABA->val;
00893 *v=0.0;
00894 for (i=0; i<n; i++){
00895 for (j=0;j<i;j++){
00896 dd+=2*x[i]*x[j]*val[j];
00897 k++;
00898 }
00899 dd+=x[i]*x[i]*val[i];
00900 k+=LDA;
00901 }
00902 *v=dd;
00903 return 0;
00904 }
00905
00906
00907
00908 static int DTRUMatEigs(void*AA,double W[],double IIWORK[], int nn1, double *mineig){
00909 dtrumat* AAA=(dtrumat*) AA;
00910 ffinteger info,INFO=0,M,N=AAA->n;
00911 ffinteger IL=1,IU=1,LDA=AAA->LDA,LDZ=LDA,LWORK,IFAIL;
00912 ffinteger *IWORK=(ffinteger*)IIWORK;
00913 double *AP=AAA->val;
00914 double Z=0,VL=-1e10,VU=1e10,*WORK,ABSTOL=1e-13;
00915 char UPLO=AAA->UPLO,JOBZ='N',RANGE='I';
00916 DSDPCALLOC2(&WORK,double,8*N,&info);
00917 DSDPCALLOC2(&IWORK,ffinteger,5*N,&info);
00918 LWORK=8*N;
00919 dsyevx(&JOBZ,&RANGE,&UPLO,&N,AP,&LDA,&VL,&VU,&IL,&IU,&ABSTOL,&M,W,&Z,&LDZ,WORK,&LWORK,IWORK,&IFAIL,&INFO);
00920
00921
00922
00923
00924 DSDPFREE(&WORK,&info);
00925 DSDPFREE(&IWORK,&info);
00926 *mineig=W[0];
00927 return INFO;
00928 }
00929
00930
00931 static struct DSDPVMat_Ops turdensematops;
00932
00933 static int DSDPDenseXInitializeOps(struct DSDPVMat_Ops* densematops){
00934 int info;
00935 if (!densematops) return 0;
00936 info=DSDPVMatOpsInitialize(densematops); DSDPCHKERR(info);
00937 densematops->matview=DTRUMatView;
00938 densematops->matscalediagonal=DTRUMatScaleDiagonal;
00939 densematops->matshiftdiagonal=DTRUMatShiftDiagonal;
00940 densematops->mataddouterproduct=DTRUMatOuterProduct;
00941 densematops->matmult=DTRUMatMult;
00942 densematops->matdestroy=DTRUMatDestroy;
00943 densematops->matfnorm2=DenseSymPSDNormF2;
00944 densematops->matgetsize=DTRUMatGetSize;
00945 densematops->matzeroentries=DTRUMatZero;
00946 densematops->matgeturarray=DTRUMatGetDenseArray;
00947 densematops->matrestoreurarray=DTRUMatRestoreDenseArray;
00948 densematops->matmineig=DTRUMatEigs;
00949 densematops->id=1;
00950 densematops->matname=lapackname;
00951 return 0;
00952 }
00953
00954 #undef __FUNCT__
00955 #define __FUNCT__ "DSDPXMatUCreateWithData"
00956 int DSDPXMatUCreateWithData(int n,double nz[],int nnz,struct DSDPVMat_Ops* *xops, void * *xmat){
00957 int i,info;
00958 double dtmp;
00959 dtrumat*AA;
00960 DSDPFunctionBegin;
00961 if (nnz<n*n){DSDPSETERR1(2,"Array must have length of : %d \n",n*n);}
00962 for (i=0;i<n*n;i++) dtmp=nz[i];
00963 info=DTRUMatCreateWData(n,n,nz,nnz,&AA); DSDPCHKERR(info);
00964 AA->owndata=0;
00965 info=DSDPDenseXInitializeOps(&turdensematops); DSDPCHKERR(info);
00966 *xops=&turdensematops;
00967 *xmat=(void*)AA;
00968 DSDPFunctionReturn(0);
00969 }
00970
00971 #undef __FUNCT__
00972 #define __FUNCT__ "DSDPXMatUCreate"
00973 int DSDPXMatUCreate(int n,struct DSDPVMat_Ops* *xops, void * *xmat){
00974 int info,nn=n*n;
00975 double *vv;
00976 DSDPFunctionBegin;
00977 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
00978 info=DSDPXMatUCreateWithData(n,vv,nn,xops,xmat);DSDPCHKERR(info);
00979 DTRUMatOwnData((dtrumat*)(*xmat),1);
00980 DSDPFunctionReturn(0);
00981 }
00982
00983 static struct DSDPDSMat_Ops tdsdensematops;
00984 static int DSDPDSDenseInitializeOps(struct DSDPDSMat_Ops* densematops){
00985 int info;
00986 if (!densematops) return 0;
00987 info=DSDPDSMatOpsInitialize(densematops); DSDPCHKERR(info);
00988 densematops->matseturmat=DDenseSetXMat;
00989 densematops->matview=DTRUMatView;
00990 densematops->matdestroy=DTRUMatDestroy;
00991 densematops->matgetsize=DTRUMatGetSize;
00992 densematops->matzeroentries=DTRUMatZero;
00993 densematops->matmult=DTRUMatMult;
00994 densematops->matvecvec=DDenseVecVec;
00995 densematops->id=1;
00996 densematops->matname=lapackname;
00997 return 0;
00998 }
00999
01000 #undef __FUNCT__
01001 #define __FUNCT__ "DSDPCreateDSMatWithArray2"
01002 int DSDPCreateDSMatWithArray2(int n,double vv[],int nnz,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
01003 int info;
01004 dtrumat*AA;
01005 DSDPFunctionBegin;
01006 info=DTRUMatCreateWData(n,n,vv,nnz,&AA); DSDPCHKERR(info);
01007 AA->owndata=0;
01008 info=DSDPDSDenseInitializeOps(&tdsdensematops); DSDPCHKERR(info);
01009 *dsmatops=&tdsdensematops;
01010 *dsmat=(void*)AA;
01011 DSDPFunctionReturn(0);
01012 }
01013
01014 #undef __FUNCT__
01015 #define __FUNCT__ "DSDPCreateXDSMat2"
01016 int DSDPCreateXDSMat2(int n,struct DSDPDSMat_Ops* *dsmatops, void**dsmat){
01017 int info,nn=n*n;
01018 double *vv;
01019 DSDPFunctionBegin;
01020 DSDPCALLOC2(&vv,double,nn,&info);DSDPCHKERR(info);
01021 info=DSDPCreateDSMatWithArray2(n,vv,nn,dsmatops,dsmat);DSDPCHKERR(info);
01022 DTRUMatOwnData((dtrumat*)(*dsmat),1);
01023 DSDPFunctionReturn(0);
01024 }
01025
01026
01027 typedef struct {
01028 int neigs;
01029 double *eigval;
01030 double *an;
01031 } Eigen;
01032
01033 typedef struct {
01034 dtrumat* AA;
01035 Eigen *Eig;
01036 } dvecumat;
01037
01038 #undef __FUNCT__
01039 #define __FUNCT__ "CreateDvecumatWData"
01040 static int CreateDvecumatWdata(int n, double vv[], dvecumat **A){
01041 int info,nnz=n*n;
01042 dvecumat* V;
01043 DSDPCALLOC1(&V,dvecumat,&info);DSDPCHKERR(info);
01044 info=DTRUMatCreateWData(n,n,vv,nnz,&V->AA); DSDPCHKERR(info);
01045 V->Eig=0;
01046 *A=V;
01047 return 0;
01048 }
01049
01050
01051 static int DvecumatGetRowNnz(void* AA, int trow, int nz[], int *nnzz,int n){
01052 int k;
01053 *nnzz=n;
01054 for (k=0;k<n;k++) nz[k]++;
01055 return 0;
01056 }
01057
01058 static int DTRUMatGetRowAdd(void* AA, int nrow, double ytmp, double row[], int n){
01059 dtrumat* A=(dtrumat*) AA;
01060 ffinteger i,nnn=n;
01061 double *v=A->val;
01062
01063 nnn=nrow*n;
01064 for (i=0;i<=nrow;i++){
01065 row[i]+=ytmp*v[nnn+i];
01066 }
01067 for (i=nrow+1;i<n;i++){
01068 nnn+=nrow;
01069 row[i]+=ytmp*v[nrow];
01070 }
01071 return 0;
01072 }
01073
01074 static int DvecumatGetRowAdd(void* AA, int trow, double scl, double r[], int m){
01075 int info;
01076 dvecumat* A=(dvecumat*)AA;
01077 info=DTRUMatGetRowAdd((void*)A->AA ,trow,scl,r,m);
01078 return 0;
01079 }
01080
01081 static int DvecumatAddMultiple(void* AA, double alpha, double r[], int nnn, int n){
01082 dvecumat* A=(dvecumat*)AA;
01083 ffinteger nn=nnn, ione=1;
01084 double *val=A->AA->val;
01085 daxpy(&nn,&alpha,val,&ione,r,&ione);
01086 return 0;
01087 }
01088
01089
01090 static int DvecuEigVecVec(void*, double[], int, double*);
01091 static int DvecumatVecVec(void* AA, double x[], int n, double *v){
01092 dvecumat* A=(dvecumat*)AA;
01093 int i,j,k=0,LDA=A->AA->LDA;
01094 double dd=0,*val=A->AA->val;
01095 *v=0.0;
01096 if (A->Eig && A->Eig->neigs<n/5){
01097 i=DvecuEigVecVec(AA,x,n,v);
01098 } else {
01099 for (i=0; i<n; i++){
01100 for (j=0;j<i;j++){
01101 dd+=2*x[i]*x[j]*val[j];
01102 }
01103 dd+=x[i]*x[i]*val[i];
01104 k+=LDA;
01105 }
01106 *v=dd;
01107 }
01108 return 0;
01109 }
01110
01111
01112 static int DvecumatFNorm2(void* AA, int n, double *v){
01113 dvecumat* A=(dvecumat*)AA;
01114 long int i,j,k=0,LDA=A->AA->LDA;
01115 double dd=0,*x=A->AA->val;
01116 for (i=0; i<n; i++){
01117 for (j=0;j<i;j++){
01118 dd+=2*x[j]*x[j];
01119 }
01120 dd+=x[i]*x[i];
01121 k+=LDA;
01122 }
01123 *v=dd;
01124 return 0;
01125 }
01126
01127
01128 static int DvecumatCountNonzeros(void* AA, int *nnz, int n){
01129 *nnz=n*(n+1)/2;
01130 return 0;
01131 }
01132
01133
01134 static int DvecumatDot(void* AA, double x[], int nn, int n, double *v){
01135 dvecumat* A=(dvecumat*)AA;
01136 double d1,dd=0,*v1=x,*v2=A->AA->val;
01137 ffinteger i,n2,ione=1,LDA=A->AA->LDA;
01138
01139 for (i=0;i<n;i++){
01140 n2=i+1;
01141 d1=ddot(&n2,v1,&ione,v2,&ione);
01142 v1+=n; v2+=LDA;
01143 dd+=d1;
01144 }
01145 *v=2*dd;
01146 return 0;
01147 }
01148
01149
01150
01151
01152
01153
01154
01155 #undef __FUNCT__
01156 #define __FUNCT__ "DvecumatDestroy"
01157 static int DvecumatDestroy(void* AA){
01158 dvecumat* A=(dvecumat*)AA;
01159 int info;
01160 info=DTRUMatDestroy((void*)(A->AA));
01161 if (A->Eig){
01162 DSDPFREE(&A->Eig->an,&info);DSDPCHKERR(info);
01163 DSDPFREE(&A->Eig->eigval,&info);DSDPCHKERR(info);
01164 }
01165 DSDPFREE(&A->Eig,&info);DSDPCHKERR(info);
01166 DSDPFREE(&A,&info);DSDPCHKERR(info);
01167 return 0;
01168 }
01169
01170
01171 static int DvecumatView(void* AA){
01172 dvecumat* A=(dvecumat*)AA;
01173 dtrumat* M=A->AA;
01174 int i,j,LDA=M->LDA;
01175 double *val=M->val;
01176 for (i=0; i<M->n; i++){
01177 for (j=0; j<M->n; j++){
01178 printf(" %4.2e",val[j]);
01179 }
01180 val+=LDA;
01181 }
01182 return 0;
01183 }
01184
01185
01186 #undef __FUNCT__
01187 #define __FUNCT__ "DSDPCreateDvecumatEigs"
01188 static int CreateEigenLocker(Eigen **EE,int neigs, int n){
01189 int info;
01190 Eigen *E;
01191
01192 DSDPCALLOC1(&E,Eigen,&info);DSDPCHKERR(info);
01193 DSDPCALLOC2(&E->eigval,double,neigs,&info);DSDPCHKERR(info);
01194 DSDPCALLOC2(&E->an,double,n*neigs,&info);DSDPCHKERR(info);
01195 E->neigs=neigs;
01196 *EE=E;
01197 return 0;
01198 }
01199
01200
01201 static int EigMatSetEig(Eigen* A,int row, double eigv, double v[], int n){
01202 double *an=A->an;
01203 A->eigval[row]=eigv;
01204 memcpy((void*)(an+n*row),(void*)v,n*sizeof(double));
01205 return 0;
01206 }
01207
01208
01209 static int EigMatGetEig(Eigen* A,int row, double *eigenvalue, double eigenvector[], int n){
01210 double* an=A->an;
01211 *eigenvalue=A->eigval[row];
01212 memcpy((void*)eigenvector,(void*)(an+n*row),n*sizeof(double));
01213 return 0;
01214 }
01215
01216
01217 static int DvecumatComputeEigs(dvecumat*,double[],int,double[],int,double[],int,int[],int);
01218
01219 static int DvecumatFactor(void*AA, double dmatp[], int nn0, double dwork[], int n, double ddwork[], int n1, int iptr[], int n2){
01220
01221 int info;
01222 dvecumat* A=(dvecumat*)AA;
01223 if (A->Eig) return 0;
01224 info=DvecumatComputeEigs(A,dmatp,nn0,dwork,n,ddwork,n1,iptr,n2);DSDPCHKERR(info);
01225 return 0;
01226 }
01227
01228 static int DvecumatGetRank(void *AA,int *rank, int n){
01229 dvecumat* A=(dvecumat*)AA;
01230 if (A->Eig){
01231 *rank=A->Eig->neigs;
01232 } else {
01233 DSDPSETERR(1,"Vecu Matrix not factored yet\n");
01234 }
01235 return 0;
01236 }
01237
01238 static int DvecumatGetEig(void* AA, int rank, double *eigenvalue, double vv[], int n, int indz[], int *nind){
01239 dvecumat* A=(dvecumat*)AA;
01240 int i,info;
01241 if (A->Eig){
01242 info=EigMatGetEig(A->Eig,rank,eigenvalue,vv,n);DSDPCHKERR(info);
01243 *nind=n;
01244 for (i=0;i<n;i++){ indz[i]=i;}
01245 } else {
01246 DSDPSETERR(1,"Vecu Matrix not factored yet\n");
01247 }
01248 return 0;
01249 }
01250
01251 static int DvecuEigVecVec(void* AA, double v[], int n, double *vv){
01252 dvecumat* A=(dvecumat*)AA;
01253 int i,rank,neigs;
01254 double *an,dd,ddd=0,*eigval;
01255 if (A->Eig){
01256 an=A->Eig->an;
01257 neigs=A->Eig->neigs;
01258 eigval=A->Eig->eigval;
01259 for (rank=0;rank<neigs;rank++){
01260 for (dd=0,i=0;i<n;i++){
01261 dd+=v[i]*an[i];
01262 }
01263 an+=n;
01264 ddd+=dd*dd*eigval[rank];
01265 }
01266 *vv=ddd;
01267 } else {
01268 DSDPSETERR(1,"Vecu Matrix not factored yet\n");
01269 }
01270 return 0;
01271 }
01272
01273
01274 static struct DSDPDataMat_Ops dvecumatops;
01275 static const char *datamatname="STANDARD VECU MATRIX";
01276
01277 static int DvecumatOpsInitialize(struct DSDPDataMat_Ops *sops){
01278 int info;
01279 if (sops==NULL) return 0;
01280 info=DSDPDataMatOpsInitialize(sops); DSDPCHKERR(info);
01281 sops->matvecvec=DvecumatVecVec;
01282 sops->matdot=DvecumatDot;
01283 sops->mataddrowmultiple=DvecumatGetRowAdd;
01284 sops->mataddallmultiple=DvecumatAddMultiple;
01285 sops->matview=DvecumatView;
01286 sops->matdestroy=DvecumatDestroy;
01287 sops->matfactor2=DvecumatFactor;
01288 sops->matgetrank=DvecumatGetRank;
01289 sops->matgeteig=DvecumatGetEig;
01290 sops->matrownz=DvecumatGetRowNnz;
01291 sops->matfnorm2=DvecumatFNorm2;
01292 sops->matnnz=DvecumatCountNonzeros;
01293 sops->id=1;
01294 sops->matname=datamatname;
01295 return 0;
01296 }
01297
01298 #undef __FUNCT__
01299 #define __FUNCT__ "DSDPGetDUmat"
01300 int DSDPGetDUMat(int n,double *val, struct DSDPDataMat_Ops**sops, void**smat){
01301 int info,k;
01302 double dtmp;
01303 dvecumat* A;
01304 DSDPFunctionBegin;
01305
01306 for (k=0;k<n*n;++k) dtmp=val[k];
01307 info=CreateDvecumatWdata(n,val,&A); DSDPCHKERR(info);
01308 A->Eig=0;
01309 info=DvecumatOpsInitialize(&dvecumatops); DSDPCHKERR(info);
01310 if (sops){*sops=&dvecumatops;}
01311 if (smat){*smat=(void*)A;}
01312 DSDPFunctionReturn(0);
01313 }
01314
01315
01316 #undef __FUNCT__
01317 #define __FUNCT__ "DvecumatComputeEigs"
01318 static int DvecumatComputeEigs(dvecumat* AA,double DD[], int nn0, double W[], int n, double WORK[], int n1, int iiptr[], int n2){
01319
01320 int i,neigs,info;
01321 long int *i2darray=(long int*)DD;
01322 int ownarray1=0,ownarray2=0,ownarray3=0;
01323 double *val=AA->AA->val;
01324 double *dmatarray=0,*dworkarray=0,eps=1.0e-12;
01325 int nn1=0,nn2=0;
01326
01327
01328 if (n*n>nn1){
01329 DSDPCALLOC2(&dmatarray,double,(n*n),&info); DSDPCHKERR(info);
01330 ownarray1=1;
01331 }
01332 memcpy((void*)dmatarray,(void*)val,n*n*sizeof(double));
01333
01334 if (n*n>nn2){
01335 DSDPCALLOC2(&dworkarray,double,(n*n),&info); DSDPCHKERR(info);
01336 ownarray2=1;
01337 }
01338
01339 if (n*n*sizeof(long int)>nn0*sizeof(double)){
01340 DSDPCALLOC2(&i2darray,long int,(n*n),&info); DSDPCHKERR(info);
01341 ownarray3=1;
01342 }
01343
01344
01345
01346 info=DSDPGetEigs(dmatarray,n,dworkarray,n*n,i2darray,n*n,
01347 W,n,WORK,n1,iiptr,n2);
01348 if (info){
01349 memcpy((void*)dmatarray,(void*)val,n*n*sizeof(double));
01350 info=DSDPGetEigs2(dmatarray,n,dworkarray,n*n,i2darray,n*n,
01351 W,n,WORK,n1,iiptr+3*n,n2-3*n); DSDPCHKERR(info);
01352 }
01353
01354
01355 for (neigs=0,i=0;i<n;i++){
01356 if (fabs(W[i])> eps ){ neigs++;}
01357 }
01358
01359 info=CreateEigenLocker(&AA->Eig,neigs,n);DSDPCHKERR(info);
01360
01361
01362 for (neigs=0,i=0; i<n; i++){
01363 if (fabs(W[i]) > eps){
01364 info=EigMatSetEig(AA->Eig,neigs,W[i],dmatarray+n*i,n);DSDPCHKERR(info);
01365 neigs++;
01366 }
01367 }
01368
01369 if (ownarray1){ DSDPFREE(&dmatarray,&info);DSDPCHKERR(info);}
01370 if (ownarray2){ DSDPFREE(&dworkarray,&info);DSDPCHKERR(info);}
01371 if (ownarray3){ DSDPFREE(&i2darray,&info);DSDPCHKERR(info);}
01372 return 0;
01373 }
01374