标题:一个关于矩阵求逆的问题
只看楼主
黑龙风月
Rank: 1
等 级:新手上路
帖 子:17
专家分:0
注 册:2014-6-6
结帖率:100%
已结贴  问题点数:10 回复次数:6 
一个关于矩阵求逆的问题
遇到的问题是矩阵的阶数最大只能到200,再大就会出错,求问如何解决内存不足的问题。一下是源代码:
#include<math.h>
#include<malloc.h>
#include<iomanip.h>
#include<time.h>
#include<stdio.h>
#include<stdlib.h>
#define N 200 //定义方阵的最大阶数为10000
#define random(x)((double)rand()/x)
//函数声明部分
int Gauss(double A[][N],double B[][N],int n);//采用部分主元的高斯消去法求方阵A的逆矩阵B

int main()
{
    int i,j;

    double a[N][N],b[N][N];
    int n;
    double x;
    double start,end;
    int test;

    x = pow(2,20);

    cout<<"采用部分主元的高斯消去法求方阵的逆矩阵!\n";
    cout<<"请输入方阵的阶数:";
    cin>>n;
    srand((int)time(0));
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            a[i][j] = random(x);
        }
    }

    start = (double)clock()/CLOCKS_PER_SEC;

    test=Gauss(a,b,n);
    end = (double)clock()/CLOCKS_PER_SEC;

    //运用高斯消去法求该矩阵的逆矩阵并输出
    if(test)
    {
        cout<<"该方阵的逆矩阵为:\n";

        for(i=0;i<n;i++)
        {
            cout<<setw(4);
            for(j=0;j<n;j++)
            {
                cout<<b[i][j]<<setw(10);
            }
            cout<<endl;
        }
    }
//end = (double)clock()/CLOCKS_PER_SEC;


    cout<<start<<endl;
    cout<<end<<endl;
    cout<<end-start<<"s"<<endl;

    return 0;
}

//----------------------------------------------
//功能: 采用部分主元的高斯消去法求方阵A的逆矩阵B
//入口参数: 输入方阵,输出方阵,方阵阶数
//返回值: 0 or 1
//----------------------------------------------
int Gauss(double A[][N], double B[][N], int n)
{

    int i, j, k;
    double max, temp;
    double t[N][N]; //临时矩阵


    //将A矩阵存放在临时矩阵t[n][n]中
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            t[i][j] = A[i][j];
        }
    }
    //初始化B矩阵为单位阵
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            B[i][j] = (i == j) ? (float)1 : 0;
        }
    }

    for (i = 0; i < n; i++)
    {
        //寻找主元
        max = t[i][i];
        k = i;
        for (j = i+1; j < n; j++)
        {
            if (fabs(t[j][i]) > fabs(max))
            {
                max = t[j][i];
                k = j;
            }
        }
        //如果主元所在行不是第i行,进行行交换
        if (k != i)
        {
            for (j = 0; j < n; j++)
            {
                temp = t[i][j];
                t[i][j] = t[k][j];
                t[k][j] = temp;
                //B伴随交换
                temp = B[i][j];
                B[i][j] = B[k][j];
                B[k][j] = temp;
            }
        }
        //判断主元是否为0, 若是, 则矩阵A不是满秩矩阵,不存在逆矩阵
        if (t[i][i] == 0)
        {
            cout << "There is no inverse matrix!";
            return 0;
        }
        //消去A的第i列除去i行以外的各行元素
        temp = t[i][i];
        for (j = 0; j < n; j++)
        {
            t[i][j] = t[i][j] / temp; //主对角线上的元素变为1
            B[i][j] = B[i][j] / temp; //伴随计算
        }
        for (j = 0; j < n; j++) //第0行->第n行
        {
            if (j != i) //不是第i行
            {
                temp = t[j][i];
                for (k = 0; k < n; k++) //第j行元素 - i行元素*j列i行元素
                {
                    t[j][k] = t[j][k] - t[i][k]*temp;
                    B[j][k] = B[j][k] - B[i][k]*temp;
                }
            }
        }
    }

    return 1;
}
搜索更多相关主题的帖子: include 源代码 double 高斯 star 
2015-03-09 22:38
天使梦魔
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:37
帖 子:564
专家分:2754
注 册:2007-8-29
得分:0 
#define N 200
不是已经固定住了吗,再调试要自己设置吧
2015-03-10 12:08
黑龙风月
Rank: 1
等 级:新手上路
帖 子:17
专家分:0
注 册:2014-6-6
得分:0 
当#define N 200设置成更大的数的时候就会出错,会提示出现一个问题导致程序关闭。我需要10000阶的矩阵,可是最大只能到200.
2015-03-10 12:58
天使梦魔
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:37
帖 子:564
专家分:2754
注 册:2007-8-29
得分:10 
你跟进一下呀,反汇编部分:
test    dword ptr [eax],eax     ; probe page.
静态变量超出,堆栈溢出。

改成new方式,但还是不要设置太大,否则会抛出std::bad_alloc,内存无法分配。
double占用8字节,8x10000x10000/1024/1024=762MB。new肯定会失败的,模板也不要用,理论上模板比new的上限还要小。
我改成new方式,你自己看着办吧

程序代码:
#include <iostream>
#include<math.h>

 #include<malloc.h>

 #include<iomanip>

 #include<time.h>

 #include<stdio.h>

 #include<stdlib.h>


 #define N 5000 //注意大小

 #define random(x)((double)rand()/x)

 //函数声明部分
int Gauss(double** A,double** B,int n);//采用部分主元的高斯消去法求方阵A的逆矩阵B
using namespace std;

 int main()

 {
     int i,j;

     double** a=new double*[N];
     for(i=0;i<N;++i)a[i]=new double[N]; 

     double** b=new double*[N];
     for(i=0;i<N;++i)b[i]=new double[N];

     int n;
     double x;
     double start,end;
     int test;

     x = pow(2.0,20);

     cout<<"采用部分主元的高斯消去法求方阵的逆矩阵!\n";
     cout<<"请输入方阵的阶数:";
     cin>>n;
     srand((int)time(0));
     for(i=0;i<n;i++)
     {
         for(j=0;j<n;j++)
         {
             a[i][j] = random(x);
         }
     }

     start = (double)clock()/CLOCKS_PER_SEC;

     test=Gauss(a,b,n);
     end = (double)clock()/CLOCKS_PER_SEC;

     //运用高斯消去法求该矩阵的逆矩阵并输出
    if(test)
     {
         cout<<"该方阵的逆矩阵为:\n";

         for(i=0;i<n;i++)
         {
             cout<<setw(4);
             for(j=0;j<n;j++)
             {
                 cout<<b[i][j]<<setw(10);
             }
             cout<<endl;
         }
     }

 //end = (double)clock()/CLOCKS_PER_SEC;


     cout<<start<<endl;
     cout<<end<<endl;
     cout<<end-start<<"s"<<endl;

     for(i=0;i<N;++i){delete[] a[i]; delete[] b[i];}
     delete[] a; delete[] b;
     return 0;

 }


 //----------------------------------------------

 //功能: 采用部分主元的高斯消去法求方阵A的逆矩阵B

 //入口参数: 输入方阵,输出方阵,方阵阶数
//返回值: 0 or 1

 //----------------------------------------------
 int Gauss(double** A, double** B, int n)

 {
     int i, j, k;
     double max, temp;

     A=new double*[N];
     for(i=0;i<N;++i)A[i]=new double[N]; 

     B=new double*[N];
     for(i=0;i<N;++i)B[i]=new double[N];

     double** t=new double*[N];
     for(i=0;i<N;++i)t[i]=new double[N];


    //将A矩阵存放在临时矩阵t[n][n]中
    for (i = 0; i < n; i++)
     {
         for (j = 0; j < n; j++)
         {
             t[i][j] = A[i][j];
         }
     }
     //初始化B矩阵为单位阵
    for (i = 0; i < n; i++)
     {
         for (j = 0; j < n; j++)
         {
             B[i][j] = (i == j) ? (float)1 : 0;
         }
     }

     for (i = 0; i < n; i++)
     {
         //寻找主元
        max = t[i][i];
         k = i;
         for (j = i+1; j < n; j++)
         {
             if (fabs(t[j][i]) > fabs(max))
             {
                 max = t[j][i];
                 k = j;
             }
         }
         //如果主元所在行不是第i行,进行行交换
        if (k != i)
         {
             for (j = 0; j < n; j++)
             {
                 temp = t[i][j];
                 t[i][j] = t[k][j];
                 t[k][j] = temp;
                 //B伴随交换
                temp = B[i][j];
                 B[i][j] = B[k][j];
                 B[k][j] = temp;
             }
         }
         //判断主元是否为0, 若是, 则矩阵A不是满秩矩阵,不存在逆矩阵
        if (t[i][i] == 0)
         {
             cout << "There is no inverse matrix!";
             return 0;
         }
         //消去A的第i列除去i行以外的各行元素
        temp = t[i][i];
         for (j = 0; j < n; j++)
         {
             t[i][j] = t[i][j] / temp; //主对角线上的元素变为1
             B[i][j] = B[i][j] / temp; //伴随计算
        }
         for (j = 0; j < n; j++) //第0行->第n行
        {
             if (j != i) //不是第i行
            {
                 temp = t[j][i];
                 for (k = 0; k < n; k++) //第j行元素 - i行元素*j列i行元素
                {
                     t[j][k] = t[j][k] - t[i][k]*temp;
                     B[j][k] = B[j][k] - B[i][k]*temp;
                 }
             }
         }
     }
     for(i=0;i<N;++i){delete[] A[i]; delete[] B[i]; delete[] t[i]; }
     delete[] A; delete[] B; delete[] t; 

     return 1;

 } 
2015-03-11 10:42
天使梦魔
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:37
帖 子:564
专家分:2754
注 册:2007-8-29
得分:0 
为毛老是一回二

还有一个,在32位机器里内存寻址上限是4G,你的代码里有5个地方new过大概有5x762MB约3个多G
虽然超过内存大小的会放在虚拟内存,但它有保护机制的,一般程序在占用1G的时候就开始启用硬盘读写机制(程序员自己做)
c++里除了动态内存的数据不是连续的,模板好象都是连续的。我之前在做修改器时发现这个问题的。所以模板就放弃吧。
2015-03-11 10:54
黑龙风月
Rank: 1
等 级:新手上路
帖 子:17
专家分:0
注 册:2014-6-6
得分:0 
非常感谢!
我不是学c++的所以不是太懂。

[ 本帖最后由 黑龙风月 于 2015-3-11 11:50 编辑 ]
2015-03-11 11:44
天使梦魔
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:37
帖 子:564
专家分:2754
注 册:2007-8-29
得分:0 
回复 6楼 黑龙风月
多大都行,但是数据处理需要自己在硬盘换算。(临时文件)
你的程序还可以改进,比如函数内的new部分完全没有必要,用指针就可以节约空间。
2015-03-11 11:52



参与讨论请移步原网站贴子:https://bbs.bccn.net/thread-442489-1-1.html




关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.161809 second(s), 8 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved