matrix.cpp

fisher判别函数分类的实现

// Matrix.cpp: implementation of the CMatrix class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Fisher.h"
#include "Matrix.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CMatrix::CMatrix()
{

}

CMatrix::~CMatrix()
{

}
/*js()函数用于计算行列式,通过递归算法实现*/ 
float  CMatrix::js(double aa[],int n)  
{
//	int z,j,k;
//	static int cnt = 0;
//    float r,total = 0; 
//    float b[ITEMNUM][ITEMNUM];/*b[N][N]用于存放，在矩阵s[N][N]中元素s[0]的余子式*/ 
//    if(n>2) 
//	{
//		for(z=0; z<n; z++) 
//		{
//			for(j=0; j<n-1; j++) 
//                for(k=0; k<n-1; k++) 
//                    if(k >= z) 
//						b[j][k] = s[j+1][k+1]; 
//                    else 
//						b[j][k] = s[j+1][k]; 
//            if(z%2 == 0) 
//				r = s[0][z]*js(b,n-1);  /*递归调用*/ 
//            else 
//				r = (-1)*s[0][z]*js(b,n-1); 
//           total = total+r; 
//		} 
//	} 
//   else if(n == 2) 
//       total = s[0][0]*s[1][1]-s[0][1]*s[1][0]; 
//   else if(n == 1)
//	   total = s[0][0];
//   return total; 
	int i,j,k,is,js,l,u,v;
    double f,det,q,d;
    f=1.0; det=1.0;
	double a[ITEMNUM*ITEMNUM];
	for(i=0;i<ITEMNUM*ITEMNUM;i++)
		a[i] = aa[i];
    for (k=0; k<=n-2; k++)
      { q=0.0;
        for (i=k; i<=n-1; i++)
        for (j=k; j<=n-1; j++)
          { l=i*n+j; d=fabs(a[l]);
	    if (d>q) { q=d; is=i; js=j;}
          }
        if (q+1.0==1.0)
          { det=0.0; return(det);}
        if (is!=k)
          { f=-f;
            for (j=k; j<=n-1; j++)
              { u=k*n+j; v=is*n+j;
                d=a[u]; a[u]=a[v]; a[v]=d;
              }
          }
        if (js!=k)
          { f=-f;
            for (i=k; i<=n-1; i++)
              { u=i*n+js; v=i*n+k;
                d=a[u]; a[u]=a[v]; a[v]=d;
              }
          }
        l=k*n+k;
        det=det*a[l];
        for (i=k+1; i<=n-1; i++)
          { d=a[i*n+k]/a[l];
            for (j=k+1; j<=n-1; j++)
              { u=i*n+j;
                a[u]=a[u]-d*a[k*n+j];
              }
          }
      }
    det=f*det*a[n*n-1];
    return(det);
} 

/*n_1()函数用于求原矩阵各元素对应的余子式，存放在数组b[N][N]中，定义为float型*/ 
void CMatrix::n_1(double s[][ITEMNUM],double b[][ITEMNUM],int n) 
{
	int z,j,k,l,m,g;
	double a[ITEMNUM][ITEMNUM]; 
    for(z=0; z<n; z++) 
	{
		l = z; 
       for(j=0; j<n; j++) 
	   { 
		   m = j; 
           for(k=0; k<n-1; k++) 
             for(g=0; g<n-1; g++) 
			 { 
				 if(g>=m && k<l) 
					 a[k][g] = s[k][g+1]; 
                 else if(k>=l && g<m) 
					 a[k][g] = s[k+1][g]; 
                 else if(k>=l && g>=m) 
					 a[k][g] = s[k+1][g+1]; 
                 else a[k][g] = s[k][g]; 
			 } 
		    b[z][j] = js(&a[0][0],n-1); 
	   } 
   } 
} 

bool CMatrix::ReMatrix(double a[][ITEMNUM])
{
	double r,temp;
	double b[ITEMNUM][ITEMNUM];
	int n = ITEMNUM, z, j;
	r = js(&a[0][0],n);
	if( r < 0.01 && r > -0.01)
		return false;
	n_1(a,b,n);
	for(z=0; z<n; z++)                                  /*求代数余子式，此时b[N][N]中存放的为原矩阵各元素对应的"代数余子式"*/ 
          for(j=0; j<n; j++) 
              if((z+j)%2!=0 && b[z][j]!=0) 
				  b[z][j] = -b[z][j]; 
    for(z=0; z<n; z++)        /*对b[N][N]转置,此时b[N][N]中存放的为原矩阵的伴随矩阵*/ 
       for(j=0; j<z; j++) 
	   {
	       temp = b[z][j]; 
           b[z][j] = b[j][z]; 
           b[j][z] = temp; 
	   } 
     for(z=0; z<n; z++)                                  /*求逆矩阵，此时b[N][N]中存放的是原矩阵的逆矩阵*/ 
          for(j=0; j<n; j++) 
             b[z][j] = b[z][j]/r;
	for(z=0; z<n; z++)
		for(j=0; j<n; j++)
			reveMatrix[z][j] = b[z][j];
	return true;
 }