作者markoo (上行下笑阿..)
看板NCTU-STAT98G
标题0723 统算 file
时间Thu Jul 23 16:48:53 2009
#include <stdlib.h>
#include <math.h>
void matrix_print(int m,int n, double **tmp)
{
int i,j;
for (i=0;i<m;i++)
{
for (j=0;j<n;j++) printf("%8.4lf ",tmp[i][j]);
printf("\n");
}
}
double **matrix_zero(int m,int n)//生成m列n行的0矩阵
{
int i,j;
double **ans;
ans = (double **) malloc (m*sizeof(*ans));
*ans = (double *) malloc (n*m*sizeof(**ans));
for (i=0;i<m;i++) *(ans+i) = *ans+n*i;
for (i=0;i<m;i++)
{
for (j=0;j<n;j++) ans[i][j] = 0;
}
return ans;
}
double **matrix_plus(int m,int n, double **tmp1, double **tmp2)
{
int i,j;
double **ans;
ans = matrix_zero(m,n);
for (i=0;i<m;i++)
{
for (j=0;j<n;j++) ans[i][j] = tmp1[i][j]+tmp2[i][j];
}
return ans;
}
double **matrix_minus(int m,int n,double **tmp1,double **tmp2)
{
int i,j;
double **ans;
ans = matrix_zero(m,n);
for (i=0;i<m;i++)
{
for (j=0;j<n;j++) ans[i][j] = tmp1[i][j]-tmp2[i][j];
}
return ans;
}
double **matrix_trans(int m,int n,double **tmp)
{
int i,j;
double **ans;
ans = matrix_zero(n,m);
for (i=0;i<n;i++)
{
for (j=0;j<m;j++) ans[i][j] = tmp[j][i];
}
return ans;
}
double **matrix_product(int m, int n, double **tmp1, int a,int b,double **tmp2)
{
if (n!=a)
{
printf("Operational error!\n");
exit(0);
}
int i,j,k;
double **ans;
ans = matrix_zero(m,b);
for (i=0;i<m;i++)
{
for (j=0;j<b;j++)
{
for (k=0;k<n;k++) ans[i][j]+=tmp1[i][k]*tmp2[k][j];
}
}
return ans;
}
double **matrix_gauss(int m,int n,double **tmp)
{
int i,j,k;
double **ans;
double temp;
ans = matrix_zero(m,n);
for (i=0;i<m;i++)
{
for (j=0;j<n;j++) ans[i][j] = tmp[i][j];
}
for (i=0;i<m;i++)
{
for (j=i+1;j<m;j++)
{
temp = ans[j][i]/ans[i][i];
for (k=i;k<n;k++) ans[j][k] = ans[j][k] - temp*ans[i][k];
}
}
return ans;
}
double matrix_det(int m,int n,double **tmp)
{
if (m!=n)
{
printf("Not a square matrix:\n");
exit(0);
}
int i;
double **temp;
double ans=1;
temp = matrix_gauss(m,n,tmp);
matrix_print(m,n,temp);
for (i=0;i<m;i++) ans*=temp[i][i];
return ans;
}
double **matrix_gauss_jordan(int m,int n,double **tmp)
{
int i,j,k;
double temp;
double **ans;
ans = matrix_gauss(m,n,tmp);
//printf("G(X):\n");
//matrix_print(m,n,ans);
for (i=0;i<m;i++)
{
temp = ans[i][i];
for (j=0;j<n;j++) ans[i][j] = ans[i][j]/temp;
}
//printf("IG(X):\n");
//matrix_print(m,n,ans);
for (i=m-1;i>0;i--)
{
for (j=i-1;j>=0;j--)
{
temp = ans[j][i];
for (k=i;k<n;k++) ans[j][k] = ans[j][k] -ans[i][k]*temp;
}
}
return ans;
}
double **matrix_identity(int m)
{
int i;
double **ans;
ans = matrix_zero(m,m);
for (i=0;i<m;i++) ans[i][i] =1;
return ans;
}
double **matrix_inverse(int m,int n, double **tmp)
{
if (m!=n)
{
printf("Not a square matrix\n");
exit(0);
}
int i,j;
double **I;
I = matrix_identity(m);
double **B;
B = matrix_zero(m,2*m);
for (i=0;i<m;i++)
{
for (j=0;j<m;j++) B[i][j] = tmp[i][j];
for (j=m;j<2*m;j++) B[i][j] = I[i][j-m];
}
matrix_print(m,2*m,B);
double **temp;
temp = matrix_gauss_jordan(m,2*m,B);
matrix_print(m,2*m,temp);
double **ans;
ans = matrix_zero(m,m);
for (i=0;i<m;i++)
{
for (j=0;j<m;j++) ans[i][j] = temp[i][j+m];
}
return ans;
}
double **matrix_delet(int m,int a,int b,double **tmp)//m阶方阵,拿掉a-th列,b-th行
{
if (a>m || b>m)
{
printf("Operational error!\n");
exit(0);
}
int i,j;
double **ans;
ans = matrix_zero(m-1,m-1);
for (i=0;i<m-1;i++)
{
for (j=0;j<m-1;j++)
{
if (i<(a-1) && j<(b-1)) ans[i][j] = tmp[i][j];
if (i>=(a-1) && j<(b-1)) ans[i][j] = tmp[i+1][j];
if (i<(a-1) && j>=(b-1)) ans[i][j] = tmp[i][j+1];
if (i>=(a-1) && j>=(b-1)) ans[i][j] = tmp[i+1][j+1];
}
}
return ans;
}
double **matrix_inverse2(int m,int n,double **tmp)
{
if (m!=n)
{
printf("Not a square matrix\n");
exit(0);
}
int i,j;
double **ans;
double det;
ans = matrix_zero(m,m);
det = matrix_det(m,m,tmp);
for (i=0;i<m;i++)
{
for (j=0;j<m;j++)
{
ans[i][j] = pow(-1.0,(i+j)*1.0)*matrix_det(m-1,m-1,matrix_delet(m,j+1,i+1,tmp))/det;
}
}
return ans;
}
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.113.7.248