DSDP
dufull.c
Go to the documentation of this file.
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