rmmat.c

Go to the documentation of this file.
00001 #include "dsdpdatamat_impl.h"
00002 #include "dsdpsys.h"
00008 typedef struct {
00009   double ev;
00010   const double *spval;
00011   const int    *spai;
00012   int    nnz;
00013   int    n;
00014   int    ishift;
00015   char   UPLQ;
00016 } r1mat;
00017 
00018 static int R1MatDestroy(void*);
00019 static int R1MatView(void*);
00020 static int R1MatVecVec(void*, double[], int, double *);
00021 static int R1MatDotP(void*, double[],int,int,double *);
00022 static int R1MatDotU(void*, double[],int,int,double *);
00023 static int R1MatGetRank(void*, int*, int);
00024 static int R1MatFactor(void*);
00025 static int R1MatGetEig(void*, int, double*, double[], int,int[],int*);
00026 static int R1MatRowNnz(void*, int, int[], int*, int);
00027 static int R1MatAddRowMultiple(void*, int, double, double[], int);
00028 static int R1MatAddMultipleP(void*, double, double[], int,int);
00029 static int R1MatAddMultipleU(void*, double, double[], int,int);
00030 
00031 static struct  DSDPDataMat_Ops r1matopsP;
00032 static struct  DSDPDataMat_Ops r1matopsU;
00033 static int R1MatOpsInitializeP(struct  DSDPDataMat_Ops*);
00034 static int R1MatOpsInitializeU(struct  DSDPDataMat_Ops*);
00035 
00036 
00037 #undef __FUNCT__  
00038 #define __FUNCT__ "DSDPGetR1Mat"
00039 int DSDPGetR1Mat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, char UPLQ, void**mmat){
00040   int i;
00041   r1mat*AA;
00042   DSDPFunctionBegin;
00043   for (i=0;i<nnz;i++){
00044     if (spai[i]-ishift<0 || spai[i]-ishift >=n){
00045       printf("Invalid entry: Entry %d . Is %d <= %d < %d?\n",i,ishift,spai[i],n+ishift);
00046       return 1;
00047     }
00048   }
00049   AA=(r1mat*) malloc(1*sizeof(r1mat));
00050   if (AA==NULL) return 1;
00051   AA->n=n;
00052   AA->UPLQ=UPLQ;
00053   AA->spval=spval;
00054   AA->spai=spai;
00055   AA->nnz=nnz;
00056   AA->ev=ev;
00057   AA->ishift=ishift;
00058   if (mmat){*mmat=(void*)AA;}
00059   DSDPFunctionReturn(0);
00060 }
00061 
00062 #undef __FUNCT__
00063 #define __FUNCT__ "DSDPGetR1PMat"
00064 
00077 int DSDPGetR1PMat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, struct  DSDPDataMat_Ops**mops, void**mmat){
00078   int info;
00079   DSDPFunctionBegin;
00080   info=DSDPGetR1Mat(n,ev,ishift,spai,spval,nnz,'P',mmat);
00081   info=R1MatOpsInitializeP(&r1matopsP); if(info){return 1;}
00082   if (mops){*mops=&r1matopsP;}
00083   DSDPFunctionReturn(0);
00084 }
00085 
00086 #undef __FUNCT__  
00087 #define __FUNCT__ "DSDPGetR1UMat"
00088 
00101 int DSDPGetR1UMat(int n, double ev, int ishift, const int spai[], const double spval[], int nnz, struct  DSDPDataMat_Ops**mops, void**mmat){
00102   int info;
00103   DSDPFunctionBegin;
00104   info=DSDPGetR1Mat(n,ev,ishift,spai,spval,nnz,'U',mmat);
00105   info=R1MatOpsInitializeU(&r1matopsU); if(info){return 1;}
00106   if (mops){*mops=&r1matopsU;}
00107   DSDPFunctionReturn(0);
00108 }
00109 
00110 static int R1MatDotP(void* A, double x[], int nn, int n, double *v){
00111   r1mat* AA = (r1mat*)A;
00112   int i,i2,i3,j,j2;
00113   int nnz=AA->nnz,ishift=AA->ishift;
00114   const int *ai=AA->spai;
00115   double dtmp=0.0,d3;
00116   const double *val=AA->spval;
00117   for (i=0;i<nnz;i++){
00118     d3=val[i];
00119     i2=ai[i]-ishift;
00120     i3=(i2+1)*i2/2;
00121     for (j=0;j<nnz;j++){
00122       j2=ai[j]-ishift;
00123       if (j2<=i2){
00124         dtmp+=2*x[i3+j2]*d3*val[j];
00125       } 
00126     }
00127   }
00128   *v=dtmp*AA->ev;
00129   return 0;
00130 }
00131 
00132 static int R1MatDotU(void* A, double x[], int nn, int n, double *v){
00133   r1mat* AA = (r1mat*)A;
00134   int i,i2,i3,j,j2;
00135   int nnz=AA->nnz,ishift=AA->ishift;
00136   const int *ai=AA->spai;
00137   const double *val=AA->spval;
00138   double dtmp=0.0,d3;
00139 
00140   for (i=0;i<nnz;i++){
00141     d3=val[i];
00142     i2=ai[i]-ishift;
00143     i3=i2*n;
00144     for (j=0;j<nnz;j++){
00145       j2=ai[j]-ishift;
00146       if (j2<=i2){
00147         dtmp+=2*x[i3+j2]*d3*val[j];
00148       } 
00149     }
00150   }
00151   *v=dtmp*AA->ev;
00152   return 0;
00153 }
00154 
00155 static int R1MatVecVec(void* A, double x[], int n, double *v){
00156 
00157   r1mat* AA = (r1mat*)A;
00158   double dtmp=0.0;
00159   const double *val=AA->spval;
00160   int i,ishift=AA->ishift,nnz=AA->nnz;
00161   const int *ai=AA->spai;
00162   for (i=0; i<nnz; i++){
00163     dtmp+=val[i] * x[ai[i]-ishift];
00164   }
00165   *v=dtmp*dtmp*AA->ev;
00166   return 0;
00167 }
00168 
00169 static int R1MatAddMultipleP(void*A, double dd, double vv[], int nn, int n){
00170   r1mat* AA = (r1mat*)A;
00171   int i,i2,i3,j,j2;
00172   int nnz=AA->nnz,ishift=AA->ishift;
00173   const int *ai=AA->spai;
00174   const double *val=AA->spval;
00175   double d3,ddd=dd*AA->ev;
00176   for (i=0;i<nnz;i++){
00177     d3=ddd*val[i];
00178     i2=ai[i]-ishift;
00179     i3=(i2+1)*i2/2;
00180     for (j=0;j<nnz;j++){
00181       j2=ai[j]-ishift;
00182       if (j2<=i2){
00183         vv[i3+j2]+=d3*val[j];
00184       } 
00185     }
00186   }
00187   return 0;
00188 }
00189 static int R1MatAddMultipleU(void*A, double dd, double vv[], int nn, int n){
00190   r1mat* AA = (r1mat*)A;
00191   int i,i2,i3,j,j2;
00192   int nnz=AA->nnz,ishift=AA->ishift;
00193   const int *ai=AA->spai;
00194   const double *val=AA->spval;
00195   double d3,ddd=dd*AA->ev;
00196   for (i=0;i<nnz;i++){
00197     d3=ddd*val[i];
00198     i2=ai[i]-ishift;
00199     i3=i2*n;      
00200     for (j=0;j<nnz;j++){
00201       j2=ai[j]-ishift;
00202       if (j2<=i2){
00203         vv[i3+j2]+=d3*val[j];
00204       } 
00205     }
00206   }
00207   return 0;
00208 }
00209 
00210 static int R1MatAddRowMultiple(void*A, int nrow, double dd, double row[], int n){
00211   r1mat* AA = (r1mat*)A;
00212   int nnz=AA->nnz,ishift=AA->ishift;
00213   const int *ai=AA->spai;
00214   const double *val=AA->spval;
00215   double ddd=dd*AA->ev;
00216   int i,j;
00217   for (i=0;i<nnz;i++){ 
00218     if (ai[i]-ishift==nrow){
00219       for (j=0;j<nnz;j++){
00220         row[ai[j]-ishift]+= ddd*val[i]*val[j];
00221       }
00222     }
00223   }
00224   return 0;
00225 }
00226 
00227 
00228 static int R1MatFactor(void*A){
00229   return 0;
00230 }
00231 
00232 
00233 static int R1MatGetRank(void *A, int*rank, int n){
00234   *rank=1;
00235   return 0;
00236 }
00237 
00238 static int R1MatGetEig(void*A, int neig, double *eig, double v[], int n, int  indx[], int*nind){
00239   r1mat* AA = (r1mat*)A;
00240   int i,aii,ishift=AA->ishift,nnz=AA->nnz;
00241   const int *ai=AA->spai;
00242   const double *val=AA->spval;
00243   for (i=0;i<n;i++){ v[i]=0.0; }
00244   *eig=0; *nind=0;
00245   if (neig==0){
00246     for (i=0;i<nnz;i++){ 
00247       aii=ai[i]-ishift;
00248       v[aii]=val[i]; 
00249       indx[i]=aii;
00250     }
00251     *eig=AA->ev; *nind=AA->nnz;
00252   }
00253   return 0;
00254 }
00255 
00256 
00257 static int R1MatRowNnz(void*A, int row, int nz[], int *rnnz, int n){
00258   r1mat* AA = (r1mat*)A;
00259   int i,j;
00260   int nnz=AA->nnz,ishift=AA->ishift;
00261   const int *ai=AA->spai;
00262   *rnnz=0;
00263   for (i=0;i<nnz;i++){ 
00264     if (ai[i]-ishift==row){
00265       for (j=0;j<nnz;j++){
00266         nz[ai[j]-ishift]++;
00267       }
00268     }
00269     *rnnz=nnz;
00270   }
00271   return 0;
00272 }
00273 
00274 static int R1MatFNorm2(void*A, int n, double *fnorm2){
00275   r1mat* AA = (r1mat*)A;
00276   double dd=0;
00277   const double *val=AA->spval;
00278   int i,nnz=AA->nnz;
00279   for (i=0;i<nnz;i++){
00280     dd+=val[i]*val[i];
00281   }
00282   *fnorm2=dd*dd*AA->ev*AA->ev;
00283   return 0;
00284 }
00285 
00286 static int R1MatCountNonzeros(void*A, int *nnz, int n){
00287   r1mat* AA = (r1mat*)A;
00288   *nnz=AA->nnz*AA->nnz;
00289   return 0;
00290 }
00291 
00292 
00293 static int R1MatView(void* A){
00294   int i;
00295   r1mat* AA = (r1mat*)A;
00296   printf("This matrix is %4.8e times the outer product of \n",AA->ev);
00297   for (i=0;i<AA->nnz;i++){
00298     printf("%d  %4.8e \n",AA->spai[i],AA->spval[i]);
00299   }
00300   return 0;
00301 }
00302 
00303 
00304 static int R1MatDestroy(void* A){
00305   if (A) free(A);
00306   return 0;
00307 }
00308 
00309 static const char *datamatname="RANK 1 Outer Product";
00310 static int R1MatOpsInitializeP(struct  DSDPDataMat_Ops* r1matops){
00311   int info;
00312   if (r1matops==NULL) return 0;
00313   info=DSDPDataMatOpsInitialize(r1matops); DSDPCHKERR(info);
00314   r1matops->matfactor1=R1MatFactor;
00315   r1matops->matgetrank=R1MatGetRank;
00316   r1matops->matgeteig=R1MatGetEig;
00317   r1matops->matvecvec=R1MatVecVec;
00318   r1matops->matdot=R1MatDotP;
00319   r1matops->mataddrowmultiple=R1MatAddRowMultiple;
00320   r1matops->mataddallmultiple=R1MatAddMultipleP;
00321   r1matops->matdestroy=R1MatDestroy;
00322   r1matops->matview=R1MatView;
00323   r1matops->matrownz=R1MatRowNnz;
00324   r1matops->matfnorm2=R1MatFNorm2;
00325   r1matops->matnnz=R1MatCountNonzeros;
00326   r1matops->id=15;
00327   r1matops->matname=datamatname;
00328   return 0;
00329 }
00330 static int R1MatOpsInitializeU(struct  DSDPDataMat_Ops* r1matops){
00331   int info;
00332   if (r1matops==NULL) return 0;
00333   info=DSDPDataMatOpsInitialize(r1matops); DSDPCHKERR(info);
00334   r1matops->matfactor1=R1MatFactor;
00335   r1matops->matgetrank=R1MatGetRank;
00336   r1matops->matgeteig=R1MatGetEig;
00337   r1matops->matvecvec=R1MatVecVec;
00338   r1matops->matdot=R1MatDotU;
00339   r1matops->mataddrowmultiple=R1MatAddRowMultiple;
00340   r1matops->mataddallmultiple=R1MatAddMultipleU;
00341   r1matops->matdestroy=R1MatDestroy;
00342   r1matops->matview=R1MatView;
00343   r1matops->matrownz=R1MatRowNnz;
00344   r1matops->matfnorm2=R1MatFNorm2;
00345   r1matops->matnnz=R1MatCountNonzeros;
00346   r1matops->id=15;
00347   r1matops->matname=datamatname;
00348   return 0;
00349 }
00350 

Generated on Wed May 21 19:31:19 2008 for DSDP by  doxygen 1.5.5