同學，以下是反矩陣補完的程式 主要是增加一個全域變數pivot儲存每次做gauss消去法時換行的次數 然後拿來計算行列式值 那這邊有問題歡迎上課或來信發問:) #include <stdlib.h> #include <math.h> int pivot; 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"); } } void freem(double **tmp) { free(*tmp); free(tmp); } 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]; } pivot=0; for (i=0;i<m-1;i++) { j=1; while (fabs(ans[i][i]) <= 0.00000001 ) { pivot++; for (k=i;k<n;k++) { temp = ans[i][k]; ans[i][k] = ans[i+j][k]; ans[i+j][k] = temp; } j++; } 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; pivot=0; temp = matrix_gauss(m,n,tmp); //matrix_print(m,n,temp); ans = pow(-1.0,pivot*1.0); 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: 118.166.52.192