www.pudn.com > ica_C.rar > matrix.c

#include 
#include 
#include 
#include "matrix.h"

void error(char *message)
{
  printf("\n%s\n",message);
  exit(0);
}

vector newvec(int n)
{
  return (SCALAR *)malloc(sizeof(SCALAR)*n);
}

matrix newmat(int nrow, int ncol)
{
  int i;
  matrix a;
  
  a=(SCALAR **)malloc((nrow+1)*sizeof(SCALAR));
  if(a==NULL) return NULL;
  for(i=0;i=0)free(a[i]);
      free(a);return(NULL);
    }
  }
  a[nrow]=NULL;
  return a;
}

vector new_vector(int n)
{
  vector v;
  
  v=newvec(n);
  if(v==NULL)error("No Memory Space\n");
  return(v);
}

matrix new_matrix(int nrow, int ncol)
{
  matrix a;

  /*printf("Request is to malloc %ldbyte.\n", nrow*ncol*sizeof(SCALAR));*/
  a=newmat(nrow, ncol);
  if(a==NULL) error("No Space\n");
  return a;
}

void free_vector(vector v)
{
  free(v);
}

void free_matrix(matrix a)
{
  matrix b;
  
  b=a;
  while(*b!=NULL) free(*b++);
  free(a);
}

void zeros(int n, int m, matrix a)
{
  int i,j;

  for(i=0;iu) u=t;
    }
    if(u==0) goto EXIT;
    weight[k]=1/u;
  }
  det=1;
  for(k=0;ku){u=t;j=i;}
    }
    ik=ip[j];
    if(j!=k){
      ip[j]=ip[k]; ip[k]=ik;
      det=-det;
    }
    u=a[ik][k];det*=u;
    if(u==0) goto EXIT;
    for(i=k+1;i=0;i--){
	t=a_inv[i][k]; ii=ip[i]; 
	for(j=i+1;j=0;i--){
	t=a_inv[i][k]; ii=ip[i]; 
	for(j=i+1;jk; i--){
      p=v[i];q=d[i]-t*p; d[i]=q;
      for(j=i;j=2){ 
    d[n-2]=a[n-2][n-2];
    e[n-2]=a[n-2][n-1];
  }
  if(n>=1)
    d[n-1]=a[n-1][n-1];
  for(k=n-1;k>=0;k--){
    v=a[k];
    if(k0;h--){
    j=h;

    while( fabs(e[j]) > (EPS *(fabs(d[j-1])+fabs(d[j]))) ){
      j--;
      if(j<=0) break;
    }
    if(j==h) continue;
    iter = 0;

    do {
      if(++iter > MAX_ITER){
	printf("Cannot get eigenvalue\n");
	exit(0);
      }
      w=(d[h-1]-d[h])/2;
      t=e[h]*e[h];
      s=(SCALAR)sqrt((double)(w*w+t)); if(w<0)s=-s;
      xx=d[j]-d[h]+t/(w+s); yy=e[j+1];

      for(k=j;k=fabs(yy)){
	  t=-yy/xx; c=(SCALAR)(1/sqrt(t*t+1));
	  s=t*c;
	}
	else {
	  t=-xx/yy; s=(SCALAR)(1/sqrt(t*t+1));
	  c=t*s;
	}
	w=d[k]-d[k+1];
	t=(w*s+2*c*e[k+1])*s;
	d[k]-=t; d[k+1]+=t;
	if(k>j)e[k]=c*e[k]-s*yy;
	e[k+1]+=s*(c*w-2*s*e[k+1]);

	for(i=0;i<3;i++){
	  xx=a[k][i]; yy=a[k+1][i];
	  a[k  ][i]=c*xx-s*yy;
	  a[k+1][i]=s*xx+c*yy;
	}
	if(kEPS*(fabs(d[h-1])+fabs(d[h])));
  }
  
  for(k=0;k<2;k++){
    h=k; t=d[h];
    for(i=k+1;i<3;i++)
      if(d[i]>t){h=i;t=d[h];}
    d[h]=d[k];d[k]=t;
    v=a[h];a[h]=a[k];a[k]=v;
    }

  *x=a[0][0];
  *y=a[0][1];
  *z=a[0][2];
  return(d[0]);
}


/*
  n x n matrix, 
  a's engenvectors will be stored in a as a'. 
  d is eigen values, e is working space.
*/
void eigen(int n, matrix a, vector d, vector e)
{
  int i,j,k,h,iter;
  SCALAR c,s,t,w,x,y;
  vector v;

  tridiagonalize(n,a,d,&e[1]);
  e[0]=0;

  for(h=n-1;h>0;h--){
    j=h;

    while( fabs(e[j]) > (EPS *(fabs(d[j-1])+fabs(d[j]))) ){
      j--;
      if(j<=0) break;
    }
    if(j==h) continue;
    iter = 0;

    do {
      if(++iter > MAX_ITER){
	printf("Cannot get eigenvalue\n");
	exit(0);
      }
      w=(d[h-1]-d[h])/2;
      t=e[h]*e[h];
      s=(SCALAR)sqrt((double)(w*w+t)); if(w<0)s=-s;
      x=d[j]-d[h]+t/(w+s); y=e[j+1];

      for(k=j;k=fabs(y)){
	  t=-y/x; c=(SCALAR)(1/sqrt((double)(t*t+1)));
	  s=t*c;
	}
	else {
	  t=-x/y; s=(SCALAR)(1/sqrt((double)(t*t+1)));
	  c=t*s;
	}
	w=d[k]-d[k+1];
	t=(w*s+2*c*e[k+1])*s;
	d[k]-=t; d[k+1]+=t;
	if(k>j)e[k]=c*e[k]-s*y;
	e[k+1]+=s*(c*w-2*s*e[k+1]);

	for(i=0;iEPS*(fabs(d[h-1])+fabs(d[h])));
  }
  
  for(k=0;kt){h=i;t=d[h];}
    d[h]=d[k];d[k]=t;
    v=a[h];a[h]=a[k];a[k]=v;
    }
}