++++++++++ 5阶行列式计算 +++++++

[minghui000]

#include <iostream.h>
#define max 5
void main()
{
int a,b,c,num[max][max],m=0,k=1,oh=4;
//输入数据
for(a=0;a<max;a++)
{
for(b=0;b<max;b++)
{
cout<<"请输入第"<<a+1<<"行第"<<b+1<<"个";
cin>>num[a][b];
}
}

cout<<"你的数目如下";
for(a=0;a<max;a++)
{
cout<<endl;
for(b=0;b<max;b++)
{
cout<<num[b][a]<<" ";
}
}

//消零法
for(b=0;b<max;b++)
{
for(a=b+1;a<oh;a++)
{
m=(-1)*(num[b][a]/num[b][b]);
if(num[b][b]==0)
{
k=0;
break;
}
//问题所在
cout<<" "<<m<<endl;
for(c=0;c<max;c++)
{
num[c][a]=num[c][a-1]*m+num[c][a];
}
}
oh=oh-1;
}

cout<<endl<<"变化后"<<endl;
for(a=0;a<max;a++)
{
cout<<endl;
for(b=0;b<max;b++)
{
cout<<num[b][a]<<" ";
}
}

for(a=0;a<max;a++)
{
k=k*num[a][a];
}
cout<<k;
}



我解决不了，请大家有空的话写一个给我参考参考，我上面那个太混乱了，也不能正确运算，谢谢

[ttlb]

up

[ttlb]

呵呵，minghui000，你的短信收到了，我以前写过一个Matrix类，把求行列式的部分贴出来显现丑
template<typename T>
double Matrix<T>::formulate()
{
assert(is_formulatable());
Matrix<T> mat(*this);
// 结果
double iResult = 1;
// 从上到下的行变换
for (int i = 0; i < m_iRow - 1; ++i)
{
// 如果主对角线上是 0， 则通过行交换使之不为零
if (0 == mat.m_pptMat[i][i])
{
for (int m = i + 1; m < mat.m_iRow; ++m)
{
if (mat.m_pptMat[m][i] != 0)
{
mat.swap_line(i, m);
iResult *= -1;
}
}
// 如果还是 0
if (0 == mat.m_pptMat[i][i])
{
return 0;
}
}
// 主对角线上不是 0
else
{
for (int j = i + 1; j < mat.m_iRow; ++j)
{
if (0 == mat.m_pptMat[j][i])
continue;
// 系数
double dMult = mat.m_pptMat[j][i] / mat.m_pptMat[i][i];
for (int k = 0; k < m_iCol; ++k)
{
mat.m_pptMat[j][k] -= mat.m_pptMat[i][k] * dMult;
}
}
}
} // for (int i = 0; i < m_iRow - 1; ++i)， 从上到下的行变换结束
// 计算结果
for (int i = 0; i < mat.m_iRow; ++i)
{
iResult *= mat.m_pptMat[i][i];
}
return iResult;
}

[ttlb]

// TtlbMatrix.h
// Version: 1.2
// Part of Ttlb Library
// Version: 1.2
// Includes Matrix manipulations
// Copyright: Ttlb(Tang Tao), 2003

#ifndef _TTLB_MATRIX_H
#define _TTLB_MATRIX_H

#include <iostream>
#include <iomanip>
#include <algorithm>
#include <memory>
#include <cassert>
#include <typeinfo>
#include <cstring>
using std::ostream;
using std::cout;
using std::endl;
using std::swap;
using std::setw;
using std::setprecision;
// 不用声明 std::strcmp

namespace ttlb
{
template<typename T>
class Matrix
{
public:
// 构造函数
explicit Matrix(int iRow = 0, int iCol = 0);
// 拷贝构造函数
Matrix(const Matrix<T> &cmMat);
// 拷贝赋值
Matrix& operator= (const Matrix<T> &cmMat);
// destructor, virtual to support derivation
virtual ~Matrix() { release(); }

// 是否为空
bool is_empty() const { return 0 == m_iRow; }
// 改变（只允许扩大）容量
void resize(int iRow, int iCol);
// 交换第 i,j 两行
void swap_line(int i, int j);
// 转置
void transpose();
// 是否可以计算行列式
bool is_formulatable() const;
// 求行列式的值
double formulate();
// 是否可以求逆
bool is_contradictable();
// 求逆矩阵
void contradict();

// 下标操作符
T* operator[] (int iRowNum) const { return m_pptMat[iRowNum]; }
// 一系列算术运算操作符重载
Matrix& operator+= (const Matrix<T> &rhs);
Matrix& operator-= (const Matrix<T> &rhs);
Matrix& operator*= (const Matrix<T> &rhs);
Matrix& operator/= (const Matrix<T> &rhs);
// 输出操作符重载
template<typename T>
friend ostream& operator<< (ostream &os, const Matrix<T> &cmMat);

public:
// 全部置 0
void reset();
// 是否为单位矩阵
bool is_identity() const;
// 取得行数
int row() const { return m_iRow; }
// 取得列数
int col() const { return m_iCol; }

private:
// 分配内存（iRow 行, iCol 列）
void create(int iRow, int iCol);
// 释放内存
void release();

private:
int m_iRow; // 行
int m_iCol; // 列
T **m_pptMat; // 指向二维数组的指针

};

[minghui000]

up~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[ttlb]

// 函数定义：

template<typename T>
inline Matrix<T>::Matrix(int iRow, int iCol)
: m_iRow(iRow), m_iCol(iCol)
{
create(m_iRow, m_iCol);
}

template<typename T>
inline Matrix<T>::Matrix(const Matrix<T> &rhs)
: m_iRow(rhs.m_iRow), m_iCol(rhs.m_iCol)
{
create(m_iRow, m_iCol);
if (0 != m_iRow && 0 != m_iCol)
for (int i = 0; i < m_iRow; ++i)
memcpy(m_pptMat[i], rhs.m_pptMat[i], sizeof(T) * m_iCol);
}

template<typename T>
inline Matrix<T>& Matrix<T>::operator= (const Matrix<T> &rhs)
{
if (&rhs != this)
{
create(rhs.m_iRow, rhs.m_iCol);
if (0 != m_iRow && 0 != m_iCol)
for (int i = 0; i < m_iRow; ++i)
memcpy(m_pptMat[i], rhs.m_pptMat[i], sizeof(T) * m_iCol);
}
return *this;
}

template<typename T>
void Matrix<T>::create(int iRow, int iCol)
{
if (0 >= iRow || 0 >= iCol)
{
m_iRow = 0;
m_iCol = 0;
m_pptMat = 0;
}
else
{
m_iRow = iRow;
m_iCol = iCol;
m_pptMat = new T*[m_iRow];
for (int i = 0; i < m_iRow; ++i)
m_pptMat[i] = new T[m_iCol];
}
}

template<typename T>
inline void Matrix<T>::release()
{
if (0 != m_iRow && 0 != m_iCol)
{
for (int i = 0; i < m_iRow; ++i)
delete[] m_pptMat[i];
delete[] m_pptMat;
m_pptMat = 0;
m_iRow = m_iCol = 0;
}
}

template<typename T>
inline void Matrix<T>::resize(int iRow, int iCol)
{
if (iRow >= m_iRow && iCol >= m_iCol)
{
Matrix<T> mat(*this);
release();
create(iRow, iCol);
for (int i = 0; i < mat.m_iRow; ++i)
{
memcpy(m_pptMat[i], mat.m_pptMat[i], mat.m_iCol * sizeof(T));
}
}
}

template<typename T>
inline void Matrix<T>::swap_line(int i, int j)
{
if (i < m_iRow && i > 0 && j < m_iRow && j > 0)
{
for (int k = 0; k < m_iCol; ++k)
{
swap(m_pptMat[i][k], m_pptMat[j][k]);
}
}
}

template<typename T>
void Matrix<T>::transpose()
{
if (m_iCol != m_iRow)
{ // row col not equal, must realloc memory
Matrix<T> mat(*this);
release();
create(mat.m_iCol, mat.m_iRow);
for (int i = 0; i < m_iRow; ++i)
{
for (int j = 0; j < m_iCol; ++j)
{
m_pptMat[i][j] = mat.m_pptMat[j][i];
}
}
}
else
{ // row equals col, neednot realloc memory, swap only
for (int i = 0; i < m_iRow - 1; ++i)
{
for (int j = i + 1; j < m_iCol; ++j)
{
swap(m_pptMat[i][j], m_pptMat[j][i]);
}
}
}
}

template<typename T>
inline bool Matrix<T>::is_formulatable() const
{
return row() == col();
}

template<typename T>
double Matrix<T>::formulate()
{
assert(is_formulatable());

Matrix<T> mat(*this);

// 结果
double iResult = 1;

// 从上到下的行变换
for (int i = 0; i < m_iRow - 1; ++i)
{
// 如果主对角线上是 0， 则通过行交换使之不为零
if (0 == mat.m_pptMat[i][i])
{
for (int m = i + 1; m < mat.m_iRow; ++m)
{
if (mat.m_pptMat[m][i] != 0)
{
mat.swap_line(i, m);
iResult *= -1;
}
}
// 如果还是 0
if (0 == mat.m_pptMat[i][i])
{
return 0;
}
}
// 主对角线上不是 0
else
{
for (int j = i + 1; j < mat.m_iRow; ++j)
{
if (0 == mat.m_pptMat[j][i])
continue;
// 系数
double dMult = mat.m_pptMat[j][i] / mat.m_pptMat[i][i];
for (int k = 0; k < m_iCol; ++k)
{
mat.m_pptMat[j][k] -= mat.m_pptMat[i][k] * dMult;
}
}
}
} // for (int i = 0; i < m_iRow - 1; ++i)， 从上到下的行变换结束
// 计算结果
for (int i = 0; i < mat.m_iRow; ++i)
{
iResult *= mat.m_pptMat[i][i];
}

return iResult;
}

[ttlb]

template<typename T>
bool Matrix<T>::is_contradictable()
{
Matrix<T> mat(*this);
return 0 != mat.formulate();
}
// before contradict -> col change -> after contradict
// mat this mat this(result)
// 1 2 3 1 0 0 1 0 0 * * *
// 4 5 6 0 1 0 0 1 0 * * *
// 7 8 9 0 0 1 0 0 1 * * *
//
template<typename T>
void Matrix<T>::contradict()
{
if (is_contradictable())
{
Matrix<T> mat(*this);
// this 置为单位矩阵
reset();
for (int i = 0; i < m_iRow; ++i)
m_pptMat[i][i] = 1;
// 从上到下的行变换
for (int i = 0; i < m_iRow - 1; ++i)
{
// 如果主对角线上是 0， 则通过行交换使之不为 0
// 因为有 is_contradictable() 的保证
// 所以，行交换之后主对角线上的元素肯定不为 0
if (0 == mat.m_pptMat[i][i])
{
for (int m = i + 1; m < mat.m_iRow; ++m)
{
if (mat.m_pptMat[m][i] != 0)
{
mat.swap_line(i, m);
swap_line(i, m);
}
}
}
// 主对角线上不是 0
else
{
for (int j = i + 1; j < mat.m_iRow; ++j)
{
if (0 == mat.m_pptMat[j][i])
continue;
// 系数
double dMult = mat.m_pptMat[j][i] / mat.m_pptMat[i][i];
for (int k = 0; k < m_iCol; ++k)
{
m_pptMat[j][k] -= m_pptMat[i][k] * dMult;
mat.m_pptMat[j][k] -= mat.m_pptMat[i][k] * dMult;
}
}
}
} // for (int i = 0; i < m_iRow - 1; ++i)， 从上到下的行变换结束

// 从下到上的行变换
// 因为已经进行了从上到下的行变换，所以，主对角线上不可能有 0
for (int i = m_iRow - 1; i > 0; --i)
{
for (int j = i - 1; j >= 0; --j)
{
if (0 == mat.m_pptMat[j][i])
continue;
double dMult = mat.m_pptMat[j][i] / mat.m_pptMat[i][i];
for (int k = m_iCol - 1; k >= 0; --k)
{
m_pptMat[j][k] -= m_pptMat[i][k] * dMult;
mat.m_pptMat[j][k] -= mat.m_pptMat[i][k] * dMult;
}
}
} // for (int i = m_iRow - 1; i > 0; --i)， 从上到下的行变换结束

// mat 各位设为 1 ，使之成为单位矩阵，此时的 this 即为结果
// 因为已经进行了从上到下的行变换，所以，主对角线上不可能有 0
for (int i = 0; i < m_iRow; ++i)
{
for (int j = 0; j < m_iCol; ++j)
{
if (0 == m_pptMat[i][j])
continue;
m_pptMat[i][j] /= mat.m_pptMat[i][i];
}
}
} // if (is_contradictable()) 结束
}

template<typename T>
inline void Matrix<T>::reset()
{
for (int i = 0; i < m_iRow; ++i)
for (int j = 0; j < m_iCol; ++j)
m_pptMat[i][j] = 0;
}

template<typename T>
bool Matrix<T>::is_identity() const
{
if (m_iRow != m_iCol)
return false;
for (int i = 0; i < m_iRow; ++i)
{
if (m_pptMat[i][i] != static_cast<T>(1))
{
return false;
}
}
return true;
}

template<typename T>
inline Matrix<T>& Matrix<T>::operator+= (const Matrix<T> &rhs)
{
if (rhs.m_iRow == m_iRow && rhs.m_iCol == m_iCol)
{
for (int i = 0; i < m_iRow; ++i)
for (int j = 0; j < m_iCol; ++j)
m_pptMat[i][j] += rhs.m_pptMat[i][j];
}
return *this;
}

template<typename T>
inline Matrix<T>& Matrix<T>::operator-= (const Matrix<T> &rhs)
{
if (rhs.m_iRow == m_iRow && rhs.m_iCol == m_iCol)
{
for (int i = 0; i < m_iRow; ++i)
for (int j = 0; j < m_iCol; ++j)
m_pptMat[i][j] -= rhs.m_pptMat[i][j];
}
return *this;
}

template<typename T>
Matrix<T>& Matrix<T>::operator*= (const Matrix<T> &rhs)
{
if ((rhs.m_iRow == m_iCol && rhs.m_iCol == m_iRow)
&& !is_empty() && !rhs.is_empty())
{
Matrix<T> mat(*this);
Matrix<T> mat2(rhs);
if (m_iRow != m_iCol)
{
release();
create(mat.m_iRow, mat2.m_iCol);
}
reset();
for (int i = 0; i < m_iRow; ++i)
for (int j = 0; j < m_iCol; ++j)
for (int k = 0; k < mat.m_iCol; ++k)
m_pptMat[i][j] += mat.m_pptMat[i][k] * mat2.m_pptMat[k][j];
}
return *this;
}
// operator*= 的另一种效率较高的实现方法如下：
// template<typename T>
// inline Matrix<T>& Matrix<T>::operator+= (const Matrix<T> &rhs)
// {
// // 防止自己乘自己：
// if (&rhs == this)
// {
// if ((rhs.m_iRow == m_iCol && rhs.m_iCol == m_iRow)
// && !is_empty() && !rhs.is_empty())
// {
// Matrix<T> mat(*this);
// Matrix<T> mat2(*this);
// if (m_iRow != m_iCol)
// {
// release();
// create(mat.m_iRow, mat2.m_iCol);
// }
// reset();
// for (int i = 0; i < m_iRow; ++i)
// for (int j = 0; j < m_iCol; ++j)
// for (int k = 0; k < mat.m_iCol; ++k)
// m_pptMat[i][j] += mat.m_pptMat[i][k] * mat2.m_pptMat[k][j];
// }
// }
// else if ((rhs.m_iRow == m_iCol && rhs.m_iCol == m_iRow)
// && !is_empty() && !rhs.is_empty())
// {
// Matrix<T> mat(*this);
// if (m_iRow != m_iCol)
// {
// release();
// create(mat.m_iRow, rhs.m_iCol);
// }
// reset();
// for (int i = 0; i < m_iRow; ++i)
// for (int j = 0; j < m_iCol; ++j)
// for (int k = 0; k < mat.m_iCol; ++k)
// m_pptMat[i][j] += mat.m_pptMat[i][k] * rhs.m_pptMat[k][j];
// }
// return *this;
// }

template<typename T>
Matrix<T>& Matrix<T>::operator/= (const Matrix<T> &rhs)
{
if (rhs.m_iRow == m_iCol && rhs.m_iCol == m_iRow)
{
for (int i = 0; i < m_iRow; ++i)
for (int j = 0; j < m_iCol; ++j)
m_pptMat[i][j] += rhs.m_pptMat[i][j];
}
return *this;
}

[ttlb]

// 友元函数:

template<typename T>
ostream& operator << (ostream & os, const Matrix<T> &rhs)
{
for (int i = 0; i < rhs.m_iRow; ++i)
{
for (int j = 0; j < rhs.m_iCol; ++j)
{
if (strcmp(typeid(rhs.m_pptMat[i][j]).name(), "double") == 0)
{
os << setw(12) << setprecision(6) << rhs.m_pptMat[i][j];
}
else
{
os << setw(10) << setprecision(6) << rhs.m_pptMat[i][j];
}
}
os << endl;
}
return os;
}

// helper 函数:

template<typename T>
inline Matrix<T> operator + (const Matrix<T> &lhs, const Matrix<T> &rhs)
{
Matrix<T> mat = lhs;
mat += rhs;
return mat;
}

template<typename T>
inline Matrix<T> operator - (const Matrix<T> &lhs, const Matrix<T> &rhs)
{
Matrix<T> mat = lhs;
mat -= rhs;
return mat;
}

template<typename T>
inline Matrix<T> operator * (const Matrix<T> &lhs, const Matrix<T> &rhs)
{
Matrix<T> mat = lhs;
mat *= rhs;
return mat;
}

template<typename T>
inline Matrix<T> operator / (const Matrix<T> &lhs, const Matrix<T> &rhs)
{
Matrix<T> mat = lhs;
mat /= rhs;
return mat;
}




}


#endif // _TTLB_MATRIX_H

[minghui000]

大哥,救命啊,你的代码对我来说简直是天书！！

[lcaamtb]

/*求矩阵行列式函数 */

double surplus(double **a,int m,int n)
{
int i,j,k,p,r;
double z,temp=1,temp1=1,s=0,s1=0;

/* 如果为2*2矩阵 */
if(n==2)

{
for(i=0;i<m;i++)
for(j=0;j<n;j++)
if((i+j)%2)
temp1*=a[i][j];
else
temp*=a[i][j];
z=temp-temp1;
} /* 结束if */

/* 如果为三阶以上方阵 */
else
{

for(k=0;k<n;k++)
{
i=0;
j=k;
for(i=0,j=k;(i<m)&&(j<n);i++,j++)
temp*=a[i][j];

if(m-i)
{
for(p=m-i,r=m-1;p>0;p--,r--)
temp*=a[r][p-1];
}/* 结束if */

s+=temp;
temp=1;
}/* 结束for */


for(k=n-1;k>=0;k--)
{
for(i=0,j=k;(i<m)&&(j>=0);i++,j--)
temp1*=a[i][j];

if(m-i)
{
for(p=m-1,r=i;r<m;p--,r++)
temp1*=a[r][p];
} /* 结束if */

s1+=temp1;
temp1=1;
} /* 结束for */

z=s-s1;
} /* 结束else */

return (z);
} /* 结束求行列式函数； */

[LinuxPanther]

写的人不少了，也用不着我了！
UP了

[MPU]

我也正在写"系统工程概论"中的"单纯形法"的C实现.

就是求极值.

N 阶行列式的高斯消元法..

顺便来看看...

[BlueRockCrystal]

Up~~~~

[datalover]

这个问题应该很简单呀，我记得我在上数值分析的时候，也编写过类似的行列式计算，我是先把行列式对角化，然后就直接M=M[1][1]*M[2][2]*。。。*M[N][N]。

你自己先想一下吧，应该是比较好做的

[minghui000]

哇。好做的话我就不帖出来问呀

[datalover]

哎~~~~~~~~~`，堕落的一代！
#define N 5
main()
{
float a[N][N],b;
int i,j,k;
for(i=0;i<N;i++)
{
printf("please input the %d col data:\n",);
for(j=0;j<N;j++)
scanf("%d",&a[i][j]);
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{ if(i==j) break;
b=-a[i][i]/a[j][i];
for(k=0;k<N;k++)
a[j][k]=a[j][k]+b*a[i][k];
}}

printf("the result is as follows:\n");
for(i=0;i<N;i++)
{ printf("\n");
for(j=0;j<N;j++)
printf("%5f",a[i][j]);
}
}

你运行一下吧，我虽然是一年前做的，但是现在还能记起来该怎么做！

基本思想就是先按列逐步消。

[datalover]

不好意思，程序中的 if(i==j) break; 应该改为if(i==j) continues;

[zhouqingyuan]

up

[datalover]

不过我的方法没有检验消除第I列的时候，没考虑A[I][I]是否为0，所以还需要改善，你应该自己加一个检验函数来检验一下，然后进行合理的换行，再继续下一步的消列。