DSDP
|
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, /* fail to allocate required space */ 00015 Inverted=3, /* indefinity is detected */ 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 printf("LDA: %d %d %d \n",n,nlda,(int)sizeof(double)); 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 /* Faster, but returns more error codes. ie. INFO>0 sometimes*/ 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 printf("n: %d, ni: %d, nd: %d\n",n,ni/n,nd/n); 00112 printf("SLOW VERSION\n"); 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 /* ABSTOL=dlamch_("Safe minimum" ); */ 00122 if (0==1){ 00123 dtrumat*M; 00124 DTRUMatCreateWData(n,n,A,n*n,&M); 00125 DTRUMatView((void*)M); 00126 } 00127 /* 00128 printf("N: %d N2: %d , ",n,n2); 00129 */ 00130 /* 00131 LWORK=26*n; LIWORK=10*n; 00132 */ 00133 /* 00134 printf("JOBZ: %c, RANGE: %c, UPLO: %c, N: %d LDA: %d \n", 00135 JOBZ,RANGE,UPLO, N,LDA); 00136 printf("VL: %4.4e, VU: %4.4e, IL: %d, IU: %d, ABSTOL: %4.4e, LDZ: %d\n", 00137 VL,VU,IL,IU,ABSTOL,LDZ); 00138 printf("LWORK: %d, LIWORK: %d\n",LWORK,LIWORK); 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 /* Works and uses less memory, but slower by a factor of about 2 or 3 */ 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 /* dsymv(&UPLO,&N,&ALPHA,AP,&LDA,X,&INCX,&BETA,Y,&INCY); */ 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 /* char UPLO=A->UPLO,TRANS='N',DIAG='N'; */ 00498 00499 if (x==0 && N>0) return 3; 00500 /* 00501 memcpy(y,x,N*sizeof(double)); 00502 dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,Y,&INCY); 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 /* char UPLO=A->UPLO,TRANS='N',DIAG='N'; */ 00521 00522 if (x==0 && N>0) return 3; 00523 /* 00524 memcpy(y,x,N*sizeof(double)); 00525 dtrmv(&UPLO,&TRANS,&DIAG,&N,AP,&LDA,Y,&INCY); 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 static void printrow(double r[], int n){int i; 00627 for (i=0;i<n;i++){printf(" %4.2e",r[i]);} printf("\n"); } 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 ffinteger LIWORK=nn1; 00922 dsyevd(&JOBZ,&UPLO,&N,AP,&LDA,W,WORK,&LWORK,IWORK,&LIWORK,&INFO); 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 static int DvecumatNormF2(void* AA, int n, double *v){ 01151 dvecumat* A=(dvecumat*)AA; 01152 return(DTRUMatNormF2((void*)(A->AA), n,v)); 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 /* create a dense array in which to put numbers */ 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 /* Call LAPACK to compute the eigenvalues */ 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 /* Count the nonzero eigenvalues */ 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 /* Copy into structure */ 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