pmatlib.c

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

/*********************
*
*  pmatlib.c -- linear algebra library 
*  using pointer-allocated matrices and vectors
*
*  Todd K. Moon, Utah State University  (and other sources)
*********************************************/
/* Copyright 2004 by Todd K. Moon
 Permission is granted to use this program/data
 for educational/research only
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include "pmatlib.h"

#define MAXSTR 5000
#define isCidstart(c) (isalpha(c) || (c == '_'))
#define isCid(c) (isalnum(c) || (c == '_'))
#define isnum(c) (isdigit(c) || (c=='-') || (c=='.'))
#define isintnum(c) (isdigit(c) || (c=='-'))
#define isfnum(c) (isdigit(c) || (c=='-') || (c=='.') || (c=='e') || (c=='E'))

#define ishexnum(c) (isdigit(c) || (c=='a') || (c=='A') || (c=='b') || 	(c=='B') || (c=='c') || (c=='C') || (c=='d') || (c=='D') || (c=='e') ||  (c=='E') || (c=='f') || (c=='F') || (c=='x') || (c== 'X'))

#define iscommentstart(c) (c=='%' || c == '#' || c == ';')
#define isquote(c) (c=='"' || c=='\'')
#define LPR '('
#define RPR ')'

static int *indx = NULL;        /* (nx1) used to record index changes in LU*/
static int mi=0;                /* length of indx */

static float *tvec1f = NULL;    /* (nx1) tempoarary vect. */
static float *tvec2f = NULL;    /* (nx1) tempoarary vect. */
static float *tvec3f = NULL;    /* (nx1) tempoarary vect. */
static double *tvec1lf = NULL;  /* (nx1) tempoarary vect. */
static double *tvec2lf = NULL;  /* (nx1) tempoarary vect. */
static double *tvec3lf = NULL;  /* (nx1) tempoarary vect. */

static int mv1f=0,mv2f=0,mv3f=0;  /* size of tvec1f, tvec2f, tvec3f */
static int mv1lf=0,mv2lf=0,mv3lf=0;  /* size of tvec1lf, tvec2lf,tvev3lf */

static int mm1f=0, nm1f=0;   /* size of tmat1f */
static int mm2f=0, nm2f=0;   /* size of tmat2f */
static int mm3f=0, nm3f=0;   /* size of tmat3f */
static int mm4f=0, nm4f=0;   /* size of tmat4f */
static int mm5f=0, nm5f=0;   /* size of tmat5f */
static int mm1lf=0, nm1lf=0;   /* size of tmat1lf */
static int mm2lf=0, nm2lf=0;   /* size of tmat2lf */
static int mm3lf=0, nm3lf=0;   /* size of tmat3lf */
static int mm4lf=0, nm4lf=0;   /* size of tmat4lf */
static int mm5lf=0, nm5lf=0;   /* size of tmat5lf */


static float **tmat1f = NULL;   /* (nxn) */
static float **tmat2f = NULL;   /* (nxn) */
static float **tmat3f = NULL;   /* (nxn) */
static float **tmat4f = NULL;   /* (mxn) */
static float **tmat5f = NULL;   /* (mxn) */


static double **tmat1lf = NULL; /* (nxn) */
static double **tmat2lf = NULL; /* (nxn) */
static double **tmat3lf = NULL; /* (nxn) */
static double **tmat4lf = NULL; /* (mxn) */
static double **tmat5lf = NULL; /* (mxn) */

static float *pf,*pf1,*pf2,*pf3; /* For pointer indexing. */
static double *pd,*pd1,*pd2,*pd3; /* For pointer indexing. */

static int pbytes(int type);
static char *pdtn(int type);
static void pmaterr(char *error_text);
static int m_foundf_( float a, float b, float e);
static int m_foundlf_( double a, double b, double e);



/* stuff to set print options */
static int intwid1 = 0,intwid2 = 0;
static int flwid1=7,flwid2=4;
static int matprintnamewnl = 0;  /* print a matrix name with a newline */
static char matleftdelim[5] = {0}, matrightdelim[5] = {0};
static char matrowdelim[5] = {'\n'};  /* row delimiter: \n or ; are typical */
static int matprintnlafter = 1;

/* set output file pointer */
static FILE *fout = NULL;

#define TEMPVEC(name,newsize,size,type) if(newsize>size) { pfree_vector((void **)&name); name = (type *)pcalloc_vector2(newsize,sizeof(type)); size = newsize; }
#define TEMPMAT(name,newm,newn,m,n,type) if(newm > m || newn > n) { pfree_matrix((void ***)&name); name = (type **)pcalloc_matrix2(newm,newn,sizeof(type)); m = newm;  n = newn;}
static double lueps = 1e-20;

/***********************************************************************/
static double at,bt,ct;
#define SVDPYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? (ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))

static double maxarg1,maxarg2;
#define SVDMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
#define SVDSIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

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

void **pcalloc_matrix(int type,int length,int width,char *name)
{
   void **array;
   register int i;
   int nbytes,bwidth;
   
   if(length == 0 || width == 0)
      return NULL;
   nbytes = pbytes(type);
   bwidth = width * nbytes;
   array = (void **)calloc(length,sizeof(void *));
   if(array == NULL) {
      fprintf(stderr,"Can't allocate matrix \"%s\" of type %s.\n",name,
              pdtn(type));
      exit(-1);
   }
   
   array[0] = (char *)calloc(length,bwidth);
   if(array[0] == NULL) {
      fprintf(stderr,"Can't allocate matrix \"%s\" of type %s.\n",name,
              pdtn(type));
      exit(-1);
   }
   for(i=1 ; i<length ; i++)
      array[i] = array[0] + i * bwidth;
   return (void **)array;
}

void **pcalloc_matrix2(int length,int width, int sizeoftype)
{
   void **array;
   register int i;
   int nbytes,bwidth;
   
   if(length == 0 || width == 0)
      return NULL;
   array = (void **)calloc(length,sizeof(void *));
   if(array == NULL) return array;
   
   bwidth = width * sizeoftype;

   array[0] = (void *)calloc(length,bwidth);
   if(array[0] == NULL) {
      free(array);
      return NULL;
   }
   for(i=1 ; i<length ; i++)
      array[i] = array[0] + i * bwidth;
   return (void **)array;
}

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

void *pcalloc_vector(int type,int length,char *name)
{
   void *array;
   
   if(length == 0)
      return NULL;
   if(NULL == (array = (void *)calloc(length,pbytes(type)))){
      fprintf(stderr,"Can't allocate array \"%s\" of type %s.\n",name,
              pdtn(type));
      exit(-1);
   }
   return array;
}

void *pcalloc_vector2(int length, int sizeoftype) 
{
   void *array;
   if(length == 0) return NULL;
   if((array = (void *)calloc(length,sizeoftype)) == NULL) {
      return NULL;
   }
   return array;
}


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

void pfree_matrix(void ***array)
{
   if(*array != NULL) {
      free(*(array)[0]);
      free(*array);
      *array = NULL;
   }
}


void pfree_vector(void **array)
{
   if(*array != NULL){
      free(*array);
      *array = NULL;
   }
}

/*****************************************************************/
void ***pcalloc_tensor(int type,int nstack,int nrow,int ncol,char *s)
{
   int i;
   void ***m;
   char errmsg[100];
   
   if(nstack == 0 || ncol == 0 || nrow == 0)
      return NULL;
   m=(void ***)malloc( (unsigned) (nstack*sizeof(void***)));
   if (m==NULL){
      fprintf(stderr,"Can't allocate tensor \"%s\" of type %s.\n",s,
              pdtn(type));
      exit(2);
   }
   errmsg[0] = 0;
   strcat(errmsg,s);
   strcat(errmsg,"  (calloc_tensor)");
   for(i = 0; i < nstack; i++) {
      m[i] = (void **)pcalloc_matrix(type,nrow,ncol,errmsg);
   }
   return m;
}


void ***pcalloc_tensor2(int nmat,int nrow,int ncol,int sizeoftype)
{
   int i;
   void ***m;

   if(nmat == 0 || ncol == 0 || nrow == 0)
      return NULL;
   m=(void ***)malloc( (unsigned) (nmat*sizeof(void***)));
   if (m==NULL) return(m);

   for(i = 0; i < nmat; i++) {
      m[i] = (void **)pcalloc_matrix2(nrow,ncol,sizeoftype);
   }
   return m;
}

void pfree_tensor(void ****m,int nstack)
{
   int i;
   
   if(*m != NULL) {
      for(i = 0; i < nstack; i++) {
         pfree_matrix(*m +i);
      }
      free(*m);
      *m = NULL;
   }
}

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

void pveczerod(int *v,int dim)
{
   int i;
   for(i = 0; i < dim; i++)  v[i] = 0;
}

void pveczerof(float *v,int dim)
{
   int i;
   for(i = 0; i < dim; i++)  v[i] = 0;
}

void pveczerolf(double *v,int dim)
{
   int i;
   for(i = 0; i < dim; i++)  v[i] = 0;
}

void pmatzerof(float **m, int r, int c)
{
   int i,j;
   for(i = 0; i < r; i++) {
      pf = m[i];
      for(j = 0; j < c; j++) pf[j] = 0;
   }
}


void pmatzerolf(double **m, int r, int c)
{
   int i,j;
   for(i = 0; i < r; i++) {
      pd = m[i];
      for(j = 0; j < c; j++) pd[j] = 0;
   }
}


void pvecconstf(float *v,float c,int dim)
{
   int i;
   for(i = 0; i < dim; i++)  v[i] = c;
}


void pvecconstlf(double *v,double c,int dim)
{
   int i;
   for(i = 0; i < dim; i++)  v[i] = c;
}

void pmatconstf(float **mat, float c,int m, int n)
{
   int i,j;
   for(i = 0; i < m; i++) {
      pf=mat[i];
      for(j = 0; j < n; j++) pf[j] = c;
   }
}

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

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

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

/*****************************************************************/
void pvecsubf(float *a,float *b, float *c,int n)
{
   int i;
   for(i =0; i < n; i++) 
      c[i] = a[i] - b[i];
}


void pvecsublf(double *a,double *b, double *c,int n)
{
   int i;
   for(i =0; i < n; i++) 
      c[i] = a[i] - b[i];
}


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

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

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

/*****************************************************************/
/* Inner product: c = a'b */
float pvecinnerf(float *a, float *b, int n)
{
   int i;
   
   float c = 0;
   for(i=0; i < n; i++) c += a[i] * b[i];
   return(c);
}

double pvecinnerlf(double *a, double *b, int n)
{
   int i;
   double c = 0; 
   for(i=0; i < n; i++) c += a[i] * b[i];
   return(c);
}


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

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

/*****************************************************************/
/* copy: a --> b */
void pveccopylf(double *a,double *b,int m)
{
   int i;
   for(i = 0; i < m; i++)
      b[i] = a[i];
}

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

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


void pveccopyc(char *a,char *b,int m)
{
   int i;
   for(i = 0; i < m; i++)
      b[i] = a[i];
}


/*****************************************************************/
/* Vector norm */
double pvecnorm2lf(double *v, int n)
{
   int i;
   double sum = 0;
   for(i=0,pd=v; i < n; i++) 
      sum += v[i] * v[i];
   return(sqrt(sum));
}

float pvecnorm2f(float *v, int n)
{
   int i;
   double sum = 0;
   for(i = 0; i < n; i++) 
      sum += v[i] * v[i];
   return(sqrt(sum));
}

/*****************************************************************/
/* convolve the sequence in1 (of n1 points) with the sequence in2
   (of n2 points) into out
   This does not do "fast" convolution (e.g., using FFTs).
   Return value is the length of the new sequence
*/
int pconvlf(double *in1, int l1, double *in2, int l2, double *out)
{
   int i,j,nout;
   nout = l1+l2-1;
   /* clear out the destination */
   memset(out,0,sizeof(double)*(nout));
   for(i = 0; i < l1; i++) {
      for(j = 0; j < l2; j++) {
         out[i+j] += in1[i]*in2[j];
      }
   }
   return nout;
}

int pconvf(float *in1, int l1, float *in2, int l2, float *out)
{
   int i,j,nout;
   nout = l1+l2-1;
   /* clear out the destination */
   memset(out,0,sizeof(float)*(nout));
   for(i = 0; i < l1; i++) {
      for(j = 0; j < l2; j++) {
         out[i+j] += in1[i]*in2[j];
      }
   }
   return nout;
}


/*****************************************************************/
int pminveclf(double *v, int n)
{
   int i,j;
   for(i=1, j=0; i < n; i++) 
      if(v[i] < v[j]) j = i;
   return j;
}

int pmaxveclf(double *v, int n)
{
   int i,j;
   for(i=1, j=0; i < n; i++) 
      if(v[i] > v[j]) j = i;
   return j;
}

int pminvecf(float *v, int n)
{
   int i,j;
   for(i=1, j=0; i < n; i++) 
      if(v[i] < v[j]) j = i;
   return j;
}

int pmaxvecf(float *v, int n)
{
   int i,j;
   for(i=1, j=0; i < n; i++) 
      if(v[i] > v[j]) j = i;
   return j;
}

/*****************************************************************/
void pmataddf(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 pmataddlf(double **a,double **b,double ** c,int m,int n)