Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате
Форум программистов > C/C++: Для новичков > Не работает обращение больших матриц


Автор: hello19 21.7.2011, 18:43
Есть матрица порядка 3638. Мне нужно найти обратную к ней. ( В последствии порядок матрицы будет около 100 000)
Есть код, вот только он не совсем рабочий. На маленьких матрицах порядка 3-4 он работает отлично, а вот на матрице порядка 3638 - не работает. В качестве результата выдает исходную матрицу.
Помогите найти косяк. И по возможности сделать код "по-шустрее"
Код

#include <iostream>
#include <fstream>
#include <stdio.h>
#include <sstream>
#include <vector>
#include <string>
#include <ctime>
#include <cmath>
#include <math.h>
#include <conio.h>

using namespace std;

template <typename T> void FreeMem(T **matr, int n);
template <typename T> void PrintMtx(T **matr, int n);
template <typename T> void SetMtx(T **matr, int n);
template <typename T> void TransponMtx(T **matr, T **tMatr, int n);
void Get_matr(double **matr, int n, double **temp_matr, int indRow, int indCol);
double Det(double **matr, int n);
 
void main()
{
    srand((unsigned)time(NULL));
    setlocale(0, "");
    int n; 
    double det;
    cout << "Введите размер матрицы: ";
    cin >> n;
    double **matr = new double * [n];
    double **obr_matr = new double * [n];
    double **tobr_matr = new double * [n];
    for(int i = 0; i < n; i++)
        {
            matr[i] = new double[n];
            obr_matr[i] = new double[n];
            tobr_matr[i] = new double[n];
        }
    
        SetMtx(matr, n);
        PrintMtx(matr, n);
        det = Det(matr, n);
        cout << "Определитель матрицы = " << det << endl;
        if(det){
                for(int i = 0; i < n; i++){
                        for(int j = 0; j < n; j++){
                                int m = n - 1;
                                double **temp_matr = new double * [m];
                                for(int k = 0; k < m; k++)
                                        temp_matr[k] = new double[m];
                                Get_matr(matr, n, temp_matr, i, j);
                                obr_matr[i][j] = pow(-1.0, i + j + 2) * Det(temp_matr, m) / det;
                                FreeMem(temp_matr, m);
                        }
                }       
        }
        else
  cout << "Т.к. определитель матрицы = 0,\nто матрица вырожденная и обратной не имеет!!!" << endl;
        //Транспонирование матрицы
        TransponMtx(obr_matr, tobr_matr, n);
        //Печать обратной матрицы после транспонирования
        PrintMtx(tobr_matr, n);
        // Освобождение памяти
        FreeMem(tobr_matr, n);
        FreeMem(matr, n);
        FreeMem(obr_matr, n);
}
//Функция транспонирования матрицы
template <typename T> void TransponMtx(T **matr, T **tMatr, int n){
    for (int i = 0; i < n; i++) 
        for (int j = 0; j < n; j++)
            tMatr[j][i] = matr[i][j];
}
//Функция освобождения памяти
template <typename T> void FreeMem(T **matr, int n)
{
        for(int i = 0; i < n; i++)
                delete [] matr[i];
        delete [] matr;
}
 
//функция заполнения матрицы
template <typename T> void SetMtx(T **matr, int n)
{
    ifstream ifs1("3.txt");
        for (int i = 0; i < n; i++)
                for (int j = 0; j < n; j++)
                {
                    ifs1 >> matr[i][j];
                }
}
 
//функция печати матрицы
template <typename T> void PrintMtx(T **matr, int n)
{
    ofstream ofs4("2.txt");
        for (int i = 0; i < n; i++){
                for (int j = 0; j < n; j++)
                    ofs4 << matr[i][j] << " ";
                ofs4 << endl;
        }                
}

//функция вычеркивания строки и столбца
void Get_matr(double **matr, int n, double **temp_matr, int indRow, int indCol)       
{
        int ki = 0; 
        for (int i = 0; i < n; i++){
                if(i != indRow){
                        for (int j = 0, kj = 0; j < n; j++){
                                if (j != indCol){
                                        temp_matr[ki][kj] = matr[i][j];
                                        kj++;
                                }
                        }
                        ki++;           
                }
        }
}

 //функция вычисления определителя матрицы
double Det(double **matr, int n)      
{
        double temp = 0;   //временная переменная для хранения определителя
        int k = 1;              //степень
        if(n < 1){
                cout<<"Неверный размер матрицы!!!" << endl;
        return 0;
    }
        else if (n == 1)
                temp = matr[0][0];
        else if (n == 2)
                temp = matr[0][0] * matr[1][1] - matr[1][0] * matr[0][1];
        else{
                for(int i = 0; i < n; i++){
                        int m = n - 1; 
                        double **temp_matr = new double * [m];
                        for(int j = 0; j < m; j++)
                                temp_matr[j] = new double [m];
                        Get_matr(matr, n, temp_matr, 0, i);
                        temp = temp + k * matr[0][i] * Det(temp_matr, m);
                        k = -k;
                        FreeMem(temp_matr, m);
                }
        }
        return temp;
}

Сами матрицы находятся в текстовых файлах. Прикрепил и их:
1.txt - матрица порядка 4 (пробная)
2.txt - сюда записывается результат работы программы
3.txt - матрица порядка 3638

Автор: Сыроежка 21.7.2011, 20:07
Я вам уже на другом форуме ответил, что мне не понятно следующее (цитирую по другому форуму):

template <typename T> void SetMtx(T **matr, int n);

double **matr = new double * [n];

SetMtx(matr, n);

Я уже в этих строчках вижу "косяк". Вопрос к вам, чему у вас равен тип T, когда вы вызываете SetMtx(matr, n);?! Точнее говоря, какой тип у параметра matr в результате инстанциации вашей функции для аргумента matr?!


Автор: borisbn 21.7.2011, 22:30
Цитата(Сыроежка @  21.7.2011,  20:07 Найти цитируемый пост)
Вопрос к вам, чему у вас равен тип T

так легко ж проверить
Код

#include <iostream>
#include <typeinfo>

using namespace std;

template <typename T>
void foo( T **matr ) {
    cout << "T : " << typeid( T ).name() << endl;
    cout << "matr : " << typeid( matr ).name() << endl;
}

int main()
{
    double ** x;
    foo( x );
}

http://liveworkspace.org/code/9195f30a7f23cd50a48c9092323daecf
меня смущает другое:
Цитата(hello19 @  21.7.2011,  18:43 Найти цитируемый пост)
В последствии порядок матрицы будет около 100 000

т.е. у тебя будет матрица 100 000 x 100 000 дублей ???
посчитай, сколько это будет занимать памяти...
я уже подсчитал - 10 000 000 000 * sizeof( double ) = 80 000 000 000 = 80 ГБт.
а ты выделяешь 4 таких массива
Цитата

matr[i] = new double[n];
obr_matr[i] = new double[n];
tobr_matr[i] = new double[n];
...
temp_matr[k] = new double[m];

это 320 ГБт памяти....
я не знаю, какое у тебя "железо", но сдаётся мне, что 320 ГБт у тебя не будет...
мож я, конечно, и ошибаюсь, но IMHO сама постановка (да и реализация) задачи - неправильная.
Think 'bout it  smile 

Автор: bsa 21.7.2011, 23:13
В свое время мне один математик (кандидат наук) рассказал, что чем больше размерность матриц, тем выше требования к точности операций (читай к количеству значащих разрядов в числе). Другими словами, чтобы результаты умножения и обращения матриц таких размеров были адекватны, необходимо оперировать с числами ЗНАЧИТЕЛЬНО большей точности, чем даже long double.

Автор: hello19 22.7.2011, 09:44
Как же тогда найти обратную матрицу для матрицы такого порядка?
Для меня главное скорость, а что же касается ресурсо, то их хватит и на работу  вышеуказанного кода. 

Автор: bsa 22.7.2011, 10:08
hello19, смотри в сторону рациональных чисел с бесконечной точностью. Но тут уже 320 ГБ оперативной памяти не отделаешься. Кстати, у меня есть большие сомнения, что твой код в том виде, в котором он нам представлен, вообще будет адекватно работать на машинах с достаточным количеством ресурсов (если не секрет, как называется супер компьютер, к которому у тебя есть доступ?) - там используются совсем другие принципы построения программ (так как фактически они имеют не один мощный процессор с гигантским объемом ОЗУ, а тысячи обычных процессоров каждый со своей памятью обычного объема).

Автор: Silent 22.7.2011, 11:58
Сидел, обедал, читал любимый форум...
Цитата

Как же тогда найти обратную матрицу для матрицы такого порядка?
Для меня главное скорость, а что же касается ресурсо, то их хватит и на работу  вышеуказанного кода. 

а потом чуть не подавился О_О
иметь возможность запихнуть в ОЗУ 320Гб данных и писать такой код?

уважаемый hello19, задачи матричной алгебры появились не вчера и не позавчера, их решали еще в бородатые годы, и есть же наработки, что мешает их посмотреть?
Чтобы не быть голословным, пример: середина 60х, вычислительная машина М-20 (аж 20 000 операций в секунду!, 300 слов оперативной памяти!, аще круть) - и решали СЛАУ практически любого порядка (200го - за 9 минут).

Автор: bsa 22.7.2011, 13:02

 ! 
bsa
hello19, прекращаем плодить одинаковые темы, ответ от этого быстрее не получишь

Автор: newbieone 22.7.2011, 13:10
Цитата(Silent @  22.7.2011,  11:58 Найти цитируемый пост)
и решали СЛАУ практически любого порядка

так программисты, сдается мне, хорошо знали курс численного решения математических задач, разделы по СЛАУ, основные теоремы (LU-разложение, формы представления)... smile 

Автор: xvr 24.7.2011, 11:06
Цитата(newbieone @  22.7.2011,  13:10 Найти цитируемый пост)
так программисты, сдается мне, хорошо знали курс численного решения математических задач, разделы по СЛАУ, основные теоремы (LU-разложение, формы представления)...

Угу, у них же не было лишних 320ГБ памяти  smile 

2 ТС - обратись в раздел по алгоритмам, там быстрее скажут в какую сторону смотреть

Автор: Фантом 24.7.2011, 13:55
Цитата(hello19 @  21.7.2011,  18:43 Найти цитируемый пост)

Помогите найти косяк. И по возможности сделать код "по-шустрее"

Это один сплошной косяк, который надо выкинуть и начать все сначала. Почитайте хотя бы Википедию, даже там можно найти описания методов обращения матриц с указанием сложности соответствующих алгоритмов. 


P.S. И "шедевры" вроде if(det) уберите.

Автор: alexvs11 24.7.2011, 14:12
надо сказать вычисление одного детерминанта зависнит надолго при таких размерностях

Автор: bsa 24.7.2011, 16:59
Цитата(alexvs11 @  24.7.2011,  14:12 Найти цитируемый пост)
надо сказать вычисление одного детерминанта зависнит надолго при таких размерностях 

А без рекурсии его реализация проблематична (без использования LU-разложения).

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)