pmatlib.c

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

{
   int i,j;
   
   for(i = 0; i < m; i++) {
      for(j=0,pd=c[i],pd1=a[i],pd2=b[i]; j < n; j++) {
         pd[j] = pd1[j] + pd2[j];
      }
   }
}

/*****************************************************************/
void pmatsublf(double **a,double **b,double ** c,int m,int n)
{
   int i,j;
   
   for(i = 0; i < m; i++) {
      for(j=0,pd=c[i],pd1=a[i],pd2=b[i]; j < n; j++) {
         pd[j] = pd1[j] - pd2[j];
      }
   }
}

void pmatsubf(float **a,float **b,float ** c,int m,int n)
{
   int i,j;
   
   for(i = 0; i < m; i++) {
      for(j=0,pf=c[i],pf1=a[i],pf2=b[i]; j < n; j++) {
         pf[j] = pf1[j] - pf2[j];
      }
   }
}


/*****************************************************************/
void pmatmatmultlf(double **a,double **b,double ** c,int m,int n)
{
   int i,j;
   
   for(i = 0; i < m; i++) {
      for(j=0,pd=c[i],pd1=a[i],pd2=b[i]; j < n; j++) {
         pd[j] = pd1[j] * pd2[j];
      }
   }
}

void pmatmatmultf(float **a,float **b, float ** c,int m,int n)
{
   int i,j;
   
   for(i = 0; i < m; i++) {
      for(j=0,pf=c[i],pf1=a[i],pf2=b[i]; j < n; j++) {
         pf[j] = pf1[j] * pf2[j];
      }
   }
}

/*****************************************************************/
void pmataddandscalelf(double **a,double **b,double d,double ** c,int m,int n)
{
   int i,j;
   
   for(i = 0; i < m; i++) {
      for(j=0,pd=c[i],pd1=a[i],pd2=b[i]; j < n; j++) {
         pd[j] = d*(pd1[j] + pd2[j]);
      }
   }
}

void pmataddandscalef(float **a,float **b,float d,float ** c,int m,int n)
{
   int i,j;
   
   for(i = 0; i < m; i++) {
      for(j=0,pf=c[i],pf1=a[i],pf2=b[i]; j < n; j++) {
         pf[j] = d*(pf1[j] + pf2[j]);
      }
   }
}

/*****************************************************************/
/* Copy matrices: a --> B */
void pmatcopylf(double **a,double **b,int m,int n)
{
   int i,j;
   for(i = 0; i < m; i++)
      for(j=0,pd=b[i],pd1=a[i]; j < n; j++)
         pd[j] = pd1[j];
}

void pmatcopyf(float **a,float **b,int m,int n)
{
   int i,j;
   for(i = 0; i < m; i++)
      for(j=0,pf=b[i],pf1=a[i]; j < n; j++)
         pf[j] = pf1[j];
}

/*****************************************************************/
void pmatscalef(float **a, float b,float **c, int n, int m)
{
   int i,j;
        
   for(i = 0; i < n; i++) {
      for(j=0,pf=c[i],pf1=a[i]; j < m; j++)
         pf[j] = pf1[j] * b;
   }
}


void pmatscalelf(double **a, double b,double **c, int n, int m)
{
   int i,j;
   
   for(i = 0; i < n; i++) {
      for(j=0,pd=c[i],pd1=a[i]; j < m; j++)
         pd[j] = pd1[j] * b;
   }
}


/*****************************************************************/
void pmattpoself(double **a,double **b,int m,int n)
   /* transpose a matrix (m x n) into b matrix, assumed to be declared as */
         /* (n x m) */
{
   int i,j;
   for(i = 0; i < m; i++) {
      for(j = 0; j < n; j++) {
         b[j][i] = a[i][j];
      }
   }
}


void pmattposef(float **a,float **b,int m,int n)
{
   int i,j;
   for(i = 0; i < m; i++) {
      for(j = 0; j < n; j++) {
         b[j][i] = a[i][j];
      }
   }
}


/*****************************************************************/
void pmatmultlf(double **a, double **b, double **c, int n, int m, int q)
/* a is nxm, b is mxq, c is nxq */
{
   int i,j,k;
   for(i=0 ; i<n ; i++){
      for(k=0,pd1=c[i] ; k<q ; k++)
         pd1[k] = 0.0;
      for(k=0,pd2=a[i] ; k<m ; k++)
         if(pd2[k] != 0.0)
            for(j=0,pd1=c[i],pd3=b[k] ; j<q ; j++)
               pd1[j] += pd2[k] * pd3[j];
   }
}


void pmatmultf(float **a, float **b, float **c, int n, int m, int q)
/* a is nxm, b is mxq, c is nxq */
{
   int i,j,k;
   for(i=0 ; i<n ; i++){
      for(k=0,pf1=c[i] ; k<q ; k++)
         pf1[k] = 0.0;
      for(k=0,pf2=a[i] ; k<m ; k++)
         if(pf2[k] != 0.0)
            for(j=0,pf1=c[i],pf3=b[k] ; j<q ; j++)
               pf1[j] += pf2[k] * pf3[j];
   }
}


/*****************************************************************/
void pmatmultTlf(double **a,double **b,double **c,int n,int m,int q,int type)
   /* type 0: c = a*b
      a is n x m
      b is m x q
      c is n x q
      
      type 1: c = a*b'
      a is n x m
      b is q x m
      c is n x q
      
      type 2: c = a'*b
      a is m x n
      b is m x q
      c is n x q
      
      type 3: c = a'*b'
      a is m x n
      b is q x m
      c is n x q
   */
{
   int i,j,k;
   double temp;
   
   switch(type) {
   case 0:              /* c = a*b  */
      for(i=0 ; i<n ; i++){
         for(k=0,pd1=c[i] ; k<q ; k++)
            pd1[k] = 0.0;
         for(k=0,pd2=a[i] ; k<m ; k++)
            if(pd2[k] != 0.0)
               for(j=0,pd1=c[i],pd3=b[k] ; j<q ; j++)
                  pd1[j] += pd2[k] * pd3[j];
      }
      break;
   case 1:                      /* c = a*b' */
      for(i=0; i < n; i++) {
         for(k=0,pd=c[i]; k < q; k++) {
            temp = 0.0;
            for(j=0,pd1=a[i],pd2=b[k]; j < m; j++) {
               temp += pd1[j] * pd2[j];
            }
            pd[k] = temp;
         }
      }
      break;
   case 2:                      /* c = a'*b  */
      for(i = 0; i < n; i++)
         for(k = 0,pd=c[i]; k < q; k++)
            pd[k] = 0.0;
      for(j = 0; j < m; j++)
         for(i = 0,pd2=a[j] ; i < n; i++)
            if(pd2[i] != 0.0)
               for(k = 0,pd=c[i],pd1=b[j] ; k < q; k++)
                  pd[k] += pd2[i] * pd1[k];
      break;
   case 3:                      /* c = a'*b' */

/* type 3: c = a'*b'
a is m x n
b is q x m
c is n x q
*/
      pd2 = a[0];
      for(i = 0; i < n; i++) {
         for(k=0,pd=c[i]; k < q; k++) {
            temp = 0;
            for(j=0,pd1=b[k]; j < m; j++,pd2+=n)
               temp += pd2[i] * pd1[j];
            pd2 = a[0];
            pd[k] = temp;
         }
      }
      break;
   default:                     /* invalid type */
      pmaterr("invalid multiplication type\n");
   }
}

void pmatmultTf(float **a,float **b,float **c,int n,int m,int q,int type)
   /* type 0: c = a*b
      a is n x m
      b is m x q
      c is n x q
      
      type 1: c = a*b'
      a is n x m
      b is q x m
      c is n x q
      
      type 2: c = a'*b
      a is m x n
      b is m x q
      c is n x q
      
      type 3: c = a'*b'
      a is m x n
      b is q x m
      c is n x q
   */
{
   int i,j,k;
   float temp;
   
   switch(type) {
   case 0:              /* c = a*b  */
      for(i=0 ; i<n ; i++){
         for(k=0,pf1=c[i] ; k<q ; k++)
            pf1[k] = 0.0;
         for(k=0,pf2=a[i] ; k<m ; k++)
            if(pf2[k] != 0.0)
               for(j=0,pf1=c[i],pf3=b[k] ; j<q ; j++)
                  pf1[j] += pf2[k] * pf3[j];
      }
      break;
   case 1:                      /* c = a*b' */
      for(i=0; i < n; i++) {
         for(k=0,pf=c[i]; k < q; k++) {
            temp = 0.0;
            for(j=0,pf1=a[i],pf2=b[k]; j < m; j++) {
               temp += pf1[j] * pf2[j];
            }
            pf[k] = temp;
         }
      }
      break;
   case 2:                      /* c = a'*b  */
      for(i = 0; i < n; i++)
         for(k = 0,pf=c[i]; k < q; k++)
            pf[k] = 0.0;
      for(j = 0; j < m; j++)
         for(i = 0,pf2=a[j] ; i < n; i++)
            if(pf2[i] != 0.0)
               for(k = 0,pf=c[i],pf1=b[j] ; k < q; k++)
                  pf[k] += pf2[i] * pf1[k];
      break;
   case 3:                      /* c = a'*b' */
      pf2 = a[0];
      for(i = 0; i < n; i++) {
         for(k=0,pf=c[i]; k < q; k++) {
            temp = 0;
            for(j=0,pf1=b[k]; j < m; j++,pf2+=n)
               temp += pf2[i] * pf1[j];
            pf2 = a[0];
            pf[k] = temp;
         }
      }
      break;
   default:                     /* invalid type */
      pmaterr("invalid multiplication type\n");
   }
}



/*****************************************************************/
void pmatvecmultf(float **a, float *b,float *c, int n, int m)
   /* c = A*b, A is (nxm), b is(mx1), c is (nx1) */
{
   int i,j;
   float temp;
   
   for(i=0; i < n; i++) {
      temp = 0.0;
      for(j=0,pf=a[i]; j < m; j++)
         temp += pf[j] * b[j];
      c[i] = temp;
   }
}

void pmatvecmultlf(double **a, double *b,double *c, int n, int m)
   /* c = A*b, A is (nxm), b is(mx1), c is (nx1) */
{
   int i,j;
   double temp;
   
   
   for(i=0; i < n; i++) {
      temp = 0;
      for(j=0,pd2=a[i]; j < m; j++)
         temp += pd2[j] * b[j];
      c[i] = temp;
   }
}


/*****************************************************************/
/* Multiply a*b', where a is mx1 and b is nx1, 
   to obtain the mxn outer product c */
void pvecouterlf(double *a, double *b, double **c, int m, int n)
{
   int i, j;
   for(i = 0; i < m; i++) {
      for(j = 0; j < n; j++) {
         c[i][j] = a[i]*b[j];
      }
   }
}

void pvecouterf(float *a, float *b, float **c, int m, int n)
{
   int i, j;
   for(i = 0; i < m; i++) {
      for(j = 0; j < n; j++) {
         c[i][j] = a[i]*b[j];
      }
   }
}


/*****************************************************************/
/* Compute c = v'Mv (a generalized vector inner product */
double pvecinnerMlf(double *v, double **M, double *w, int n)
{
   double sum = 0;
   int i,j;
   
   for(i=0; i < n; i++)
      for(j=0,pd2=M[i]; j < n; j++)
         sum += v[i] * pd2[j] * w[j];
   return(sum);
}

float pvecinnerMf(float *v, float **M, float *w, int n)
{
   float sum = 0;
   int i,j;
   
   for(i=0; i < n; i++)
      for(j=0,pf2=M[i]; j < n; j++)
         sum += v[i] * pf2[j] * w[j];
   return(sum);
}

/*****************************************************************/
/* Multiply a vector (transpose) times a matrix: c = a'*B */
/* a is (nx1), b is (nxm), c is (1xm) (a vector) */
void pvecmatmultf(float *a, float **b,float *c, int n, int m)
{
   int i,j;
   
   for(j=0; j < n; j++)
      c[j] = 0.0;
   for(j=0; j < n; j++){
      for(i=0,pf=b[j]; i < m; i++)
         c[i] += pf[i] * a[j];
   }
}

void pvecmatmultlf(double *a, double **b,double *c, int n, int m)
{
   int i,j;
        
   for(j=0; j < n; j++)
      c[j] = 0.0;
   for(j=0; j < n; j++){
      for(i=0,pd=b[j]; i < m; i++)
         c[i] += pd[i] * a[j];
   }
}


/*****************************************************************/
/* Multiply a vector (transpose) times a matrix transpose: c = a'*B' */
/* a is (nx1), b is (mxn), c is (1xm) (a vector) */
void pvecmatTmultlf(double *a, double **b,double *c, int n, int m)
{
   int i,j;
   double temp;
   
   for(i=0; i < m; i++) {
      temp = 0;
      for(j=0,pd=b[i]; j < n; j++)
         temp += a[j] * pd[j];
      c[i] = temp;
   }
}

void pvecmatTmultf(float *a, float **b,float *c, int n, int m)
{
   int i,j;
   float temp;
   
   for(i=0; i < m; i++) {
      temp = 0;
      for(j=0,pf=b[i]; j < n; j++)
         temp += a[j] * pf[j];
      c[i] = temp;
   }
}



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

/* Multiply a matrix transpose times a vector: c = A'*b, 
   A is (nxm), b is(nx1), c is (mx1) */
void pmatTvecmultf(float **a, float *b,float *c, int n, int m)
{
   int i,j;
   
   for(j=0; j < m; j++)
      c[j] = 0.0;
   for(j=0; j < n; j++){
      for(i=0,pf=a[j]; i < m; i++)
         c[i] += pf[i] * b[j];
   }
}

void pmatTvecmultlf(double **a, double *b,double *c, int n, int m)
{
   int i,j;
   
   for(j=0; j < m; j++)
      c[j] = 0;
   for(j=0; j < n; j++) {
      for(i = 0, pd=a[j]; i < m; i++)
         c[i] += pd[i] * b[j];
   }
}


/* form the product c = a * diag(d), where a is nxn */
void pmatdiagmultf(float **a, float *d, float **c, int m, int n)
{
   int i,j;
   for(i = 0; i < m; i++) {
      for(j = 0,pf1 = a[i],pf2 = c[i]; j < n; j++,pf1++,pf2++) {
         *pf2 = *pf1 * d[j];
      }
   }
}

void pmatdiagmultlf(double **a, double *d, double **c, int m, int n)
{
   int i,j;
   for(i = 0; i < m; i++) {
      for(j = 0,pd1 = a[i],pd2 = c[i]; j < n; j++,pd1++,pd2++) {
         *pd2 = *pd1 * d[j];
      }
   }
}


/* Find the trace of a matrix. */
float pmattracef(float **A, int M)
{
   int i;
   float trace = 0;
   
   for(i = 0; i < M; i++) trace += A[i][i];
   return(trace);
}

double pmattracelf(double **A, int M)
{
   int i;
   double trace = 0;
   
   for(i = 0; i < M; i++) trace += A[i][i];
   return(trace);
}

/*********************************************************
Finds the frobenius norm of a matrix.
Where: