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