pmatlib.c

压缩文件中是Error Correction Coding - Mathematical Methods and Algorithms(Wiley 2005)作者:(Todd K. Moon )的配

   fprintf(fout,"%s",matleftdelim);
   for(i = 0; i < n; i++)
      fprintf(fout,"%*.*g\t",flwid1,flwid2,v[i]);
   fprintf(fout,"%s",matrightdelim);
   if(matprintnlafter) fprintf(fout,"\n");
}




/*****************************************************************/

void ptenprint(char *name, void ***ten,int nstack,int n, int m, int type)
{
   switch(type) {
   case TKCHAR:
      ptenprints(name,(char ***)ten,nstack,n);
      break;
   case TKSCHAR:
      ptenprintsc(name,(signed char ***)ten,nstack,n,m);
      break;
   case TKUCHAR: 
      ptenprintuc(name,(unsigned char ***)ten,nstack,n,m);
      break;
   case TKSHINT:
      ptenprinthd(name,(short ***)ten,nstack,n,m);
      break;
   case TKINT:
      ptenprintd(name,(int ***)ten,nstack,n,m);
      break;
   case TKFLOAT:
      ptenprintf(name,(float ***)ten,nstack,n,m);
      break;
   case TKDOUBLE:
      ptenprintlf(name,(double ***)ten,nstack,n,m);
      break;
   case TKLINT:
      ptenprintld(name,(long ***)ten,nstack,n,m);
      break;
   case TKHEX:
      ptenprintx(name,(int ***)ten,nstack,n,m);
      break;
   default:
      fprintf(stderr,"Type %d not implemented for ptenprint\n", type);
   }
}


void ptenprints(char *name,char ***ten,int nstack,int m)
{
   int i,j;

   if(fout == NULL) fout = stdout;   
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < m; j++) {
	 fprintf(fout,"%s\n",ten[i][j]);
      }
      fprintf(fout,"\n");
   }
}

/*****************************************************************/
void ptenprintsc(char *name,signed char ***ten,int nstack,int n,int m)
{
   int i,j,k;
   

   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*d\t",intwid1,intwid2,(int)ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}

/*****************************************************************/
void ptenprintuc(char *name,unsigned char ***ten,int nstack,int n,int m)
{
   int i,j,k;
   

   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*d\t",intwid1,intwid2,(int)ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}


/*****************************************************************/
void ptenprinthd(char *name,short ***ten,int nstack,int n,int m)
{
   int i,j,k;
   
   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*hd\t",intwid1,intwid2,ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}


/*****************************************************************/
void ptenprintd(char *name,int ***ten,int nstack,int n,int m)
{
   int i,j,k;
   
   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*d\t",intwid1,intwid2,ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}

/*****************************************************************/
void ptenprintx(char *name,int ***ten,int nstack,int n,int m)
{
   int i,j,k;
   
   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*x\t",intwid1,intwid2,ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}


/*****************************************************************/
void ptenprintf(char *name,float ***ten,int nstack,int n,int m)
{
   int i,j,k;
   
   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*g\t",flwid1,flwid2,ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}


/*****************************************************************/
void ptenprintlf(char *name,double ***ten,int nstack,int n,int m)
{
   int i,j,k;
   
   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*g\t",flwid1,flwid2,ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}



/*****************************************************************/
void ptenprintld(char *name,long ***ten,int nstack,int n,int m)
{
   int i,j,k;

   if(fout == NULL) fout = stdout;
   if(strlen(name)) fprintf(fout,"%s",name);
   if(matprintnamewnl) fprintf(fout,"\n");
   if(ten == NULL) {
      fprintf(fout,"(null)\n");
      return;
   }
   for(i = 0; i < nstack; i++) {
      for(j = 0; j < n; j++) {
	 for(k = 0; k < m; k++) {
	    fprintf(fout,"%*.*ld\t",intwid1,intwid2,ten[i][j][k]);
	 }
	 fprintf(fout,"\n");
      }
      fprintf(fout,"\n");
   }
}


/*****************************************************************/
/*****************************************************************/
/* compute the determinant */
float pdetf(float **a,int n)
{
   float d;
   int j;
   
   TEMPVEC(indx,n,mi,int);
   TEMPMAT(tmat1f,n,n, mm1f, nm1f, float);

   pludcmpf(a,tmat1f,n,indx,&d,lueps);	/* Lu decompose the matrix */
   
   /* for compatibility with the previous library, also compute the */
   /* determinant */
   for(j = 0; j < n; j++) {
      d *= tmat1f[j][j];
   }
   return(d);
}

double pdetlf(double **a,int n)
{
   double d;
   int j;
   
   TEMPVEC(indx,n,mi,int);
   TEMPMAT(tmat1lf,n,n, mm1lf, nm1lf, double);

   pludcmplf(a,tmat1lf,n,indx,&d,lueps);	/* Lu decompose the matrix */
   
   /* for compatibility with the previous library, also compute the */
   /* determinant */
   for(j = 0; j < n; j++) {
      d *= tmat1lf[j][j];
   }
   return(d);
}



/*****************************************************************/
void pmatsolvef(float **A, float *x, float *b, int n)
/* solve Ax = b */
{
   float d;
   int i;
   for(i = 0; i < n; i++) {
      x[i] = b[i];
   }
   TEMPMAT(tmat1f,n,n,mm1f,nm1f,float);
   TEMPVEC(indx,n,mi,int);
   pludcmpf(A,tmat1f,n,indx,&d,lueps);
   plubksubf(tmat1f,n,indx,x);
}


void pmatresolvef(float *x, float *b, int n)
/* Solve Ax = b, where A is the same as the last call to pmatsolve
   (and tmat1 has not been modified) */
{
   int i;
   for(i = 0; i < n; i++) {
      x[i] = b[i];
   }
   plubksubf(tmat1f,n,indx,x);
}


void pmatsolvelf(double **A, double *x, double *b, int n)
/* solve Ax = b */
{
   double d;
   int i;
   for(i = 0; i < n; i++) {
      x[i] = b[i];
   }
   TEMPMAT(tmat1lf,n,n,mm1lf,nm1lf,double);
   TEMPVEC(indx,n,mi,int);
   pludcmplf(A,tmat1lf,n,indx,&d,lueps);
   plubksublf(tmat1lf,n,indx,x);
}


void pmatresolvelf(double *x, double *b, int n)
/* Solve Ax = b, where A is the same as the last call to pmatsolve
   (and tmat1 has not been modified) */
{
   int i;
   for(i = 0; i < n; i++) {
      x[i] = b[i];
   }
   plubksublf(tmat1lf,n,indx,x);
}

/*****************************************************************/
float pmatinvf(float **a,float **ainv, int n)
{
   float d;
   int i,j;
   
   TEMPVEC(indx,n,mi,int);
   /* TEMPVEC(tvec1,n,mv1f,float);  <-- done by pludcmp */
   TEMPMAT(tmat1f,n,n,mm1f,nm1f,float);

   pludcmpf(a,ainv,n,indx,&d,lueps);	/* Lu decompose the matrix */
   
   /* for compatibility with the previous library, also compute the */
   /* determinant */
   for(j = 0; j < n; j++) {
      d *= ainv[j][j];
   }
   if(d) {
      for(j = 0; j < n; j++) {
	 for(i = 0; i < n; i++) tvec1f[i] = 0;
	 tvec1f[j] = 1.; 
	 plubksubf(ainv,n,indx,tvec1f);
	 for(i = 0; i < n; i++)
	    tmat1f[i][j] = tvec1f[i];
      }
      /* now copy the inverse into the called matrix */
      for(i = 0; i < n; i++) 
	 for(j = 0,pf=ainv[i],pf1=tmat1f[i]; j < n; j++)
	    pf[j] = pf1[j];
   }
   return(d);
} /* end of inverse */

void psetlueps(double luepsin)
{  lueps = luepsin;
}

double pmatinvlf(double **a,double **ainv, int n)
{
   double d;
   int i,j;
   register double *p,*p1;

   TEMPVEC(indx,n,mi,int);
   /* TEMPVEC(tvec1lf,n,mv1lf,double);  <-- done by pludcmp */
   TEMPMAT(tmat1lf,n,n,mm1lf,nm1lf,double);
   
   pludcmplf(a,ainv,n,indx,&d,lueps);	/* Lu decompose the matrix */
   
   /* for compatibility with the previous library, also compute the */
   /* determinant */
   for(j = 0; j < n; j++) {
      d *= ainv[j][j];
   }
   for(j=0,p1=tvec1lf; j < n; j++,p1++) {
      for(i=0,p=tvec1lf; i < n; i++,p++) *p = 0;
      *p1 = 1.; 
      plubksublf(ainv,n,indx,tvec1lf);
      for(i=0,p=tvec1lf; i < n; i++,p++)
	 tmat1lf[i][j] = *p;
   }
   /* now copy the inverse into the called matrix */
   for(i = 0; i < n; i++) 
      for(j=0,p=ainv[i],p1=tmat1lf[i]; j < n; j++,p++,p1++)
	 *p = *p1;
   return(d);
} /* end of inverse */


/*****************************************************************/



void pludcmpf(float **a,float **LUa,int n,int *indx,float *d,float TINY)
{
   int i,imax,j,k;
   float big,dum,sum,temp;
   
   TEMPVEC(tvec1f,n,mv1f,float);
   
   imax = 0;
   *d=1.0;
   for (i=0,pf2=tvec1f;i<n;i++,pf2++) {
      big=0.0;
      for (j=0,pf=LUa[i],pf1=a[i];j<n;j++,pf++,pf1++) {
		 *pf = *pf1;
		 if ((temp=fabs(*pf1)) > big) big=temp;
      }
      if (big  < TINY) pmaterr("Error: singular matrix in pludcmpf");
      *pf2 = 1.0/big;
   }
   for (j=0;j<n;j++) {
      for (i=0;i<j;i++) {
	 pf = LUa[i];
	 sum=pf[j];
	 for (k=0;k<i;k++) sum -= pf[k]*LUa[k][j];
	 pf[j]=sum;
      }
      big=0.0;
      for (i=j;i<n;i++) {
	 pf = LUa[i];
	 sum=pf[j];
	 for (k=0;k<j;k++)
	    sum -= pf[k]*LUa[k][j];
	 pf[j]=sum;
	 if ( (dum=tvec1f[i]*fabs(sum)) >= big) {
	    big=dum;
	    imax=i;
	 }
      }
      if (j != imax) {
	 for (k=0,pf=LUa[j];k<n;k++,pf++) {
	    dum=LUa[imax][k];
	    LUa[imax][k] = *pf;
	    *pf = dum;
	 }
	 *d = -(*d);
	 tvec1f[imax]=tvec1f[j];
      }
      indx[j]=imax;
      if (LUa[j][j] == 0) {
	 LUa[j][j] = TINY;
      }
      if (j != n - 1) {
	 dum=1.0/(LUa[j][j]);
	 for (i=j+1;i<n;i++) LUa[i][j] *= dum;
      }
   }
}



void pludcmplf(double **a,double **LUa,int n,int *indx,double *d,double TINY)
{
   int i,imax,j,k;
   double big,dum,sum,temp;
   
   TEMPVEC(tvec1lf,n,mv1lf,double);
   
   imax = 0;
   *d=1.0;
   for (i=0,pd2=tvec1lf;i<n;i++,pd2++) {
      big=0.0;
      for (j=0,pd=LUa[i],pd1=a[i];j<n;j++,pd++,pd1++) {
	 *pd = *pd1;
	 if ((temp=fabs(*pd1)) > big) big=temp;
      }
      if (big  < TINY) pmaterr("Error: singular matrix in pludcmplf");
      *pd2 = 1.0/big;
   }
   for (j=0;j<n;j++) {
      for (i=0;i<j;i++) {
	 pd = LUa[i];
	 sum=pd[j];
	 for (k=0;k<i;k++) sum -= pd[k]*LUa[k][j];
	 pd[j]=sum;
      }
      big=0.0;
      for (i=j;i<n;i++) {
	 pd = LUa[i];
	 sum=pd[j];
	 for (k=0;k<j;k++)
	    sum -= pd[k]*LUa[k][j];
	 pd[j]=sum;
	 if ( (dum=tvec1lf[i]*fabs(sum)) >= big) {
	    big=dum;
	    imax=i;
	 }
      }
      if (j != imax) {
	 for (k=0,pd=LUa[j];k<n;k++,pd++) {
	    dum=LUa[imax][k];
	    LUa[imax][k] = *pd;
	    *pd = dum;
	 }
	 *d = -(*d);
	 tvec1lf[imax]=tvec1lf[j];
      }
      indx[j]=imax;
      if (LUa[j][j] == 0) {
	 LUa[j][j] = TINY;
      }
      if (j != n - 1) {
	 dum=1.0/(LUa[j][j]);
	 for (i=j+1;i<n;i++) LUa[i][j] *= dum;
      }
   }
}