题目描述
这次,alice梦回大一上的线代课堂。她梦到线代老师在考试中出了这道压轴题:
(12分)已知一个n阶方阵,以及一个列向量A,要求编程求解:
(1)满足的列向量x(如果不存在唯一解,输出"No Solution!").(3分)
(2)矩阵的行列式det(A).(5分)
(3)矩阵的秩r(A),以及的零空间的维度dim(N(A)).(4分)
对于第(1)问,alice很清楚这就是高斯消元,因此她很快就在网上查到了解决这个问题的伪代码:
const double eps = 1e-8; int guass_elimination(double A[maxn][maxn], int n) // 高斯消元 { int i, j, k, r; for(i = 0; i < n; i++) { // 消元 r = i; for(j = i + 1; j < n; j++) if(fabs(A[j][i]) > fabs(A[r][i])) r = j; if(fabs(A[r][i]) < eps) return 0; // 无解 if(r != i) for(j = 0; j <= n; j++) swap(A[i][j], A[r][j]); for(k = i + 1; k < n; k++) { double f = A[k][i] / A[i][i]; for(j = i; j <= n; j++) A[k][j] -= f * A[i][j]; } } for(i = n - 1; i >= 0; i--) { // 回代 for(j = i + 1; j < n; j++) { A[i][n] -= A[j][n] * A[i][j]; } A[i][n] /= A[i][i]; } return 1; }
其中,fabs(a)表示浮点数a的绝对值,swap(a,b)的作用是交换两个数a,b。函数传递的两个参数分别是增广矩阵(A|b)以及行数。这个增广矩阵在函数作用之后,假设只存在唯一解,则矩阵的第(n+1)列就是符合题意的列向量x。
然而对于第(2)(3)问,则是一筹莫展。
眼看就要交卷了,这时rawslee通过心灵感应给了她指点:
"密切关注题目的第一问。"
在这个提示下,alice重新审视了上面的伪代码,只进行了一些微小的调整以及增加一些语句,便马上写出了第(2)(3)问,愉快的在梦中取得了满分。你知道她是怎么做出来的吗?
输入格式
第一行一个正整数T,表示数据组数。
对于每组数据:
第一行一个正整数n,表示矩阵A的大小为n行n列。
之后n行每行n个实数(最多有2位小数),表示矩阵A,第i行第j列表示矩阵的元素Aij。
输出格式
对于每组数据,输出三个数,分别代表det(A),(保留2位小数),r(A),dim(N(A)).
每组输出以换行符结尾。
数据范围
.
样例输入1
1 3 1 3 4 1 4 7 9 3 2
样例输出1
38.00 3 0
样例输入2
1
3
1.50 3.50 4.50
1.50 3.50 4.50
9.20 3.20 2.20
样例输出2
0.00 2 1
友情提醒
为保持精度,建议使用double类型进行计算。
# include<stdio.h>
# include<math.h>
# include<stdlib.h>
/*
(12分)已知一个阶方阵,以及一个列向量,要求编程求解:
(1)满足的列向量(如果不存在唯一解,输出"No Solution!").(3分)
(2)矩阵的行列式.(5分)
(3)矩阵的秩,以及的零空间的维度.(4分)
对于第(1)问,很清楚这就是高斯消元,因此她很快就在网上查到了解决这个问题的伪代码:
const double eps = 1e-8;
int guass_elimination(double A[maxn][maxn], int n) // 高斯消元
{
int i, j, k, r;
for(i = 0; i < n; i++) { // 消元
r = i;
for(j = i + 1; j < n; j++) if(fabs(A[j][i]) > fabs(A[r][i])) r = j;
if(fabs(A[r][i]) < eps) return 0; // 无解
if(r != i) for(j = 0; j <= n; j++) swap(A[i][j], A[r][j]);
for(k = i + 1; k < n; k++) {
double f = A[k][i] / A[i][i];
for(j = i; j <= n; j++) A[k][j] -= f * A[i][j];
}
}
for(i = n - 1; i >= 0; i--) { // 回代,求方程右边的向量
for(j = i + 1; j < n; j++) {
A[i][n] -= A[j][n] * A[i][j];
}
A[i][n] /= A[i][i];
}
return 1;
}
其中,fabs(a)表示浮点数a的绝对值,swap(a,b)的作用是交换两个数a,b。
函数传递的两个参数分别是增广矩阵以及行数。这个增广矩阵在函数作用之后,
假设只存在唯一解,则矩阵的第列就是符合题意的列向量。
然而对于第(2)(3)问,则是一筹莫展。
眼看就要交卷了,这时通过心灵感应给了她指点:
"密切关注题目的第一问。"
在这个提示下,重新审视了上面的伪代码,
只进行了一些微小的调整以及增加一些语句,便马上写出了第(2)(3)问
输入格式
第一行一个正整数,表示数据组数。
对于每组数据:
第一行一个正整数,表示矩阵的大小为行列。
之后行每行个实数(最多有2位小数),表示矩阵,第行第列表示矩阵的元素。
输出格式
对于每组数据,输出三个数,分别代表(保留2位小数),,.
每组输出以换行符结尾。
样例输入1
1 3 1 3 4 1 4 7 9 3 2
样例输出1
38.00 3
样例输入2
1
3
1.50 3.50 4.50
1.50 3.50 4.50
9.20 3.20 2.20
样例输出2
0.00 2 1
友情提醒
为保持精度,建议使用double类型进行计算。
*/
int guass_elimination(double *A,int n, double * det, int *dim) // 高斯消元
{//提示中的代码是对增广矩阵操作,而这里只要对方阵操作就行,所以
//循环上限有所修改,对增广部分的代码也删了
int i, j, k, r, flag = 1, time = 0;//time交换行的次数
double temp; //交换行,行列式取相反数
for(i = 0; i < n; i++)
{ // 消元
r = i;
for(j = i + 1; j < n; j++)//按列从左向右遍历
if(fabs(A[j*n+i]) > fabs(A[r*n+i]))//第i列中从对角线往下找最大的那个
r = j; //其所在行为r
if(fabs(A[r*n+i]) < 1e-8) //若最大的都是0
flag = 0; // 无数解,降秩
if(r != i)
{
for(j = 0; j < n; j++)
{ //将第r行与第i行交换
temp = A[i*n+j];
A[i*n+j] = A[r*n+j];
A[r*n+j] = temp;
}
time ++;
}
for(k = i + 1; k < n; k++) //将第i行以下的减f倍第i行(把最左边的
{ //即第i列)全消为0
if(fabs(A[i*n+i]) > 1e-8)
{
double f = A[k*n+i] / A[i*n+i];//f存第i列要减的倍数
for(j = i; j < n; j++)
A[k*n+j] -= f * A[i*n+j];
}
}
}//从而化为阶梯型 (注意是方阵
*det = pow(-1, time);
if(flag == 1)
{
for(int m = 0; m < n; m ++)
{
*det *= A[m*n+m];
}
}
else
{
for(int m = 0; m < n; m ++)
{
if(fabs(A[m*n+m]) > 1e-8)
{
*dim = *dim + 1;//在对角线上数维数
}
}
}
return flag;
}
int main()
{
int num;
scanf("%d", &num);//数据组数
int flag[num], dim[num], n[num];//flag是否满秩,dim为秩,n为总维数
double det[num];//行列式
for(int i = 0; i<num;i++)
{
det[i] = 1;
dim[i]=0;
}
for(int q = 0; q < num; q ++)
{
scanf("%d", &n[q]);
double *a= (double*) malloc (n[q]*n[q] * sizeof(double));
//用malloc而不用可变二维数组
//因为可变二维数组传入函数需要在函数原型中确定列数,但列数非
//常数,又函数原型在主函数获得列数前(此时列数未定义)
//故此时无法用二维数组的模式定位,不如用一维
for(int p = 0; p < n[q]*n[q]; p ++)
{
scanf("%lf", (a + p));//a为地址
}
flag[q] = guass_elimination(a, n[q], &det[q], &dim[q]);
//函数出口最多一个返回值
}
for(int q = 0;q < num; q ++)
{
if(flag[q] != 0)//满秩
{
printf("%.2lf %d %d\n", det[q], n[q], 0);
}
else
{
printf("0.00 %d %d\n", dim[q], n[q] - dim[q]);
}
}
return 0;
}