Main Page | Modules | Alphabetical List | Data Structures | Directories | File List | Globals | Related Pages

dufull.c

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

Generated on Sat Oct 15 11:05:42 2005 for DSDP by  doxygen 1.4.2