//
//////////////////////////////////////////////////////////////////////
#if !defined(AFX_MATRIXCALCULATE_H__CCBC1D7D_4466_4E8B_87DD_0A98B462C18D__INCLUDED_)
#define AFX_MATRIXCALCULATE_H__CCBC1D7D_4466_4E8B_87DD_0A98B462C18D__INCLUDED_
#if _MSC_VER >1000
#pragma once
#endif // _MSC_VER >1000
////求矩阵matrix的行列式值，n为维数
float CalculateLiner(float *matrix,int n)
////求矩阵matrix的第i行，j列的代数余子式,n为维数
float CalculateCofacter(float *matrix,int i,int j,int n)
////matrixAT=(matrixA)T，m,n为matrixA的行、列数
void CalculateAT(float *matrixA,float *matrixAT,int m,int n)
////matrixAB=matrixA*matrixB，i,j为matrixA的行、列数，j,k为为matrixB的行、列数
void CalculateAB(float *matrixA,float *matrixB,float *matrixAB,int i,int j,int k)
////matrixATA=(matrixA)T*matrixA，m,n分别为matrixA的行、列数
void CalculateATA(float *matrixA,float *matrixATA,int m,int n)
////matrixA_为matrixA的逆，m为matrixA的行、列数
void CalculateA_(float *matrixA,float *matrixA_,int m)
///*矩阵求逆子程序(Good)*/
int Invers_matrix(float *m1,int n)
////求正定矩阵a的逆矩，n为阶数
int MatrixInversion(float *a, int n)
void MatInversion(float *a,int n)
////解矩阵方程matrixA*matrixX=matrixL,m，n分别为matrixA矩阵的行，列数
void EquationResolution(float *matrixA,float *matrixL,float *matrixX,int m,int n)
#endif // !defined(AFX_MATRIXCALCULATE_H__CCBC1D7D_4466_4E8B_87DD_0A98B462C18D__INCLUDED_)
// Matrix.cpp: implementation of the matrix calculation functions.
//
//////////////////////////////////////////////////////////////////////
#include "StdAfx.h"
#include "Matrix.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__
#define new DEBUG_NEW
#endif
////求矩阵matrix的行列式值，n为维数
float CalculateLiner(float *matrix,int n)
{
float liner=0
int i=0,j=0,k=0
int p=0,q=0
if(n==1)
return matrix[0]
else
{
float *tr=(float *)calloc((n-1)*(n-1),sizeof(float))
for(k=0k<nk++)
{
p=0
for(i=0i<ni++)
{
if(i!=0)
{
q=0
for(j=0j<nj++)
{
if(j!=k)
{
tr[p*(n-1)+q]=matrix[i*n+j]
q++
}
}
p++
}
}
liner+=matrix[k]*pow(-1,k)*CalculateLiner(tr,n-1)
}
free(tr)
return liner
}
}
////求矩阵matrix的第i行，j列的代数余子式,n为维数
float CalculateCofacter(float *matrix,int i,int j,int n)
{
int x=0,y=0
int p=0,q=0
if(n==1)
return matrix[0]
else
{
float *tr=(float *)calloc((n-1)*(n-1),sizeof(float))
p=0
for(x=0x<nx++)
{
if(x!=i)
{
q=0
for(y=0y<ny++)
{
if(y!=j)
{
tr[p*(n-1)+q]=matrix[x*n+y]
q++
}
}
p++
}
}
float cc=pow(-1,i+j)*CalculateLiner(tr,n-1)
free(tr)
return cc
}
}
////matrixAT=(matrixA)T，m,n为matrixA的行、列数
void CalculateAT(float *matrixA,float *matrixAT,int m,int n)
{
for (int i=0i<mi++)
{
for (int j=0j<nj++)
{
matrixAT[j*m+i]=matrixA[i*n+j]
}
}
}
////matrixAB=matrixA*matrixB，i,j为matrixA的行、列数，j,k为为matrixB的行、列数
void CalculateAB(float *matrixA,float *matrixB,float *matrixAB,int i,int j,int k)
{
for (int m=0m<im++)
{
for (int n=0n<kn++)
{
matrixAB[m*k+n]=0
for (int s=0s<js++)
{
matrixAB[m*k+n]+=matrixA[m*j+s]*matrixB[s*k+n]
}
}
}
}
////matrixATA=(matrixA)T*matrixA，m,n为分别为matrixA的行、列数
void CalculateATA(float *matrixA,float *matrixATA,int m,int n)
{
float *at=(float *)calloc((m)*(n),sizeof(float))
CalculateAT(matrixA,at,m,n)
CalculateAB(at,matrixA,matrixATA,n,m,n)
free(at)
}
////matrixA_为matrixA的逆，m为matrixA的行、列数
void CalculateA_(float *matrixA,float *matrixA_,int m)
{
float liner=CalculateLiner(matrixA,m)
for(int i=0i<mi++)
{
for(int j=0j<mj++)
matrixA_[j*m+i]=CalculateCofacter(matrixA,i,j,m)/liner
}
}
////////////////////////////////////////////////////////////////////
////求正定矩阵a的逆矩，n为阶数
int MatrixInversion(float *a, int n)
{
int i, j, k, m
float w, g, *b
b = new float [n]
for(k = 0k <= n - 1k++)
{
w = a[0]
w=a[0]+1.0e-15
/*
if(fabs(w)+1.0 == 1.0)
{
delete b
printf("fail\n")
return(-2)
}
*/
m = n - k - 1
for(i = 1i <= n - 1i++)
{
g = a[i * n]
b[i] = g / w
if(i <= m)
b[i] = -b[i]
for(j = 1j <= ij++)
a[(i - 1) * n + j - 1] = a[i * n + j] + g * b[j]
}
a[n * n - 1] = 1.0 / w
for(i = 1i <= n - 1i++)
a[(n - 1) * n + i - 1] = b[i]
}
for(i = 0i <= n - 2i++)
for(j = i + 1j <= n - 1j++)
a[i * n + j] = a[j * n + i]
delete b
return(2)
}
////求正定矩阵a的逆矩，n为阶数
void MatInversion(float *a,int n)
{
int i,j,k
for(k=0k<nk++)
{
for(i=0i<ni++)
{
if(i!=k) *(a+i*n+k) = -*(a+i*n+k)/(*(a+k*n+k))
}
*(a+k*n+k)=1/(*(a+k*n+k))
for(i=0i<ni++)
{
if(i!=k)
{
for(j=0j<nj++)
{
if(j!=k) *(a+i*n+j) += *(a+k*n+j)* *(a+i*n+k)
}
}
}
for(j=0j<nj++)
{
if(j!=k) *(a+k*n+j)*=*(a+k*n+k)
}
}
}
/*矩阵求逆子程序(Good)*/
int Invers_matrix(float *m1,int n)
{
int *is,*js
int i,j,k,l,u,v
float temp,max_v
is=(int *)malloc(n*sizeof(int))
js=(int *)malloc(n*sizeof(int))
if(is==NULL||js==NULL)
{
printf("out of memory!\n")
return(0)
}
for(k=0k<nk++)
{
max_v=0.0
for(i=ki<ni++)
{
for(j=kj<nj++)
{
temp=fabs(m1[i*n+j])
if(temp>max_v)
{
max_v=tempis[k]=ijs[k]=j
}
}
}
if(max_v==0.0)
{
free(is)free(js)
printf("invers is not availble!\n")
return(0)
}
if(is[k]!=k)
{
for(j=0j<nj++)
{
u=k*n+jv=is[k]*n+j
temp=m1[u]m1[u]=m1[v]m1[v]=temp
}
}
if(js[k]!=k)
for(i=0i<ni++)
{
u=i*n+kv=i*n+js[k]
temp=m1[u]m1[u]=m1[v]m1[v]=temp
}
l=k*n+k
m1[l]=1.0/m1[l]
for(j=0j<nj++)
{
if(j!=k)
{
u=k*n+j
m1[u]*=m1[l]
}
}
for(i=0i<ni++)
{
if(i!=k)
{
for(j=0j<nj++)
{
if(j!=k)
{
u=i*n+j
m1[u]-=m1[i*n+k]*m1[k*n+j]
}
}
}
}
for(i=0i<ni++)
{
if(i!=k)
{
u=i*n+k
m1[u]*=-m1[l]
}
}
}
for(k=n-1k>=0k--)
{
if(js[k]!=k)
for(j=0j<nj++)
{
u=k*n+jv=js[k]*n+j
temp=m1[u]m1[u]=m1[v]m1[v]=temp
}
if(is[k]!=k)
for(i=0i<ni++)
{
u=i*n+kv=i*n+is[k]
temp=m1[u]m1[u]=m1[v]m1[v]=temp
}
}
free(is)free(js)
return(1)
}
void EquationResolution(float *matrixA,float *matrixL,float *matrixX,int m,int n)
{
if (m<n) return
float *at=(float *)malloc((m)*(n)*sizeof(float))
float *ata=(float *)malloc((n)*(n)*sizeof(float))
float *atl=(float *)malloc((n)*sizeof(float))
CalculateATA(matrixA,ata,m,n)
MatrixInversion(ata,n)
CalculateAT(matrixA,at,m,n)
CalculateAB(at,matrixL,atl,n,m,1)
CalculateAB(ata,atl,matrixX,n,n,1)
free(at)
free(ata)
free(atl)
}
#include <stdio.h>
void main()
{
int a[4][4],b[4][4],c[4][4]
int i,j,kint s[4][4]={0}
printf("请输入矩阵A:\n")
for(i=0i<=3i++)
{
for(j=0j<=3j++)
{
printf("a[%d][%d]=",i,j)
scanf("%d",&a[i][j])
}
}
printf("\n")
printf("请输入矩阵B:\n")
for(i=0i<=3i++)
{
for(j=0j<=3j++)
{
printf("a[%d][%d]=",i,j)
scanf("%d",&b[i][j])
}
}
printf("矩阵A:\n")
for(i=0i<=3i++)
{
for(j=0j<4j++)
{
printf("%d ",a[i][j])
if(j==3)
printf("\n")
}
}
printf("\n")
printf("矩阵B:\n")
for(i=0i<4i++)
{
for(j=0j<4j++)
{
printf("%d ",b[i][j])
if(j==3)
printf("\n")
}
}
printf("\n")
printf("矩阵C=A*B:\n")
for(i=0i<4i++)
{
for(j=0j<4j++)
{
for(k=0k<4k++)
{
s[i][j]+=a[i][k]*b[k][j]
}
c[i][j]=s[i][j]
printf("%d ",c[i][j])
if(j==3)
printf("\n")
}
}
}
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）