2015-12-12 2 views
0

У меня есть около 1000 строк кода, которые я написал в C для решения линейного программирования (алгоритм внутренней точки). Я понял, что мне нужно использовать Eigen для вычисления инверсии матрицы, поэтому теперь я запускаю свой код на C++ (работает просто отлично, похоже). Теперь у меня есть куча массивов, объявленных в формате C, например: A[30][30];Использование массивов C-стиля с Eigen для Matrix Inverse

В моей программе, я кучу матричных вычислений, а затем нужно найти обратную матрицу в какой-то момент, давайте назовем это матрица L[30][30] , Чтобы использовать Эйгена, мне нужно, чтобы иметь его в специальном формате Эйгеном матрицы для вызова функции m.inverse так:

//cout << "The inverse of L is:\n" << L.inverse() << endl; 

Моя цель состоит в том, чтобы найти путь ... В любом случае, чтобы получить мои данные L в формат, который Eigen примет, чтобы я мог запускать эту вещь. Я провел последние 2 часа, исследуя это, и ничего не придумал. :-(Я довольно новичок в C, поэтому, пожалуйста, будьте настолько тщательны, насколько сможете. Я хочу, чтобы был самый простой способ. Я читал о сопоставлениях, но я не очень четко говорю о указателях (что, кажется, неотъемлемая часть). Есть ли способ просто прокрутить каждую строку и столбец и скопировать их в собственную матрицу?

В то время как я прошу, мне нужно взять результирующую Eigen-матрицу и вернуть ее обратно в C массив? Как бы этот процесс работал? Спасибо заранее за любую помощь! Я потратил около 50-60 часов на это, и это должно произойти на этой неделе! Это последнее, что мне нужно сделать, и я покончу с моим долгосрочный проект. Математический класс, поэтому для меня немного нечеткая идея, но я многому учусь.

Возможно, релевантная информация: - Воспроизведение на Windows 10 i7-процессоре Sony VAIO -Компиляция с CodeBlocks на C++, но изначально написанная на C . Этот код находится в цикле while, который может быть повторен до 10 раз или около того. -Матричный обратный должен быть рассчитан для этой матрицы L на каждой итерации, и данные будут разными каждый раз.

Пожалуйста, помогите! Я хочу учиться, но мне нужно руководство, и этот класс в сети, поэтому у меня практически нет. Большое спасибо заранее!

Edit - я видел это и пытался осуществить это без толку, но это, кажется, как и решение, если я могу понять это:

«Предположим, у вас есть массив с двойными значениями размера NROWS х Ncols.

double *X; // non-NULL pointer to some data 

Вы можете создать NROWS х Ncols размера двойную матрицу с помощью функции карты, как это:

MatrixXd eigenX = Map<MatrixXd>(X, nRows, nCols); 

операция карты отображает существующую Regio памяти n в структуры данных Эйгена. Одна строка, подобная этой, позволяет избежать написания уродливого кода создания матрицы, для цикла с копированием каждого элемента в хорошем состоянии и т. Д. ».

Это похоже на приятное решение, но я не знаю, как это сделать. что «double * X» говорит, что «указывает на некоторые данные». Я начал искать указатели и т. д., и это не помогло прояснить - я видел всевозможные вещи, указывающие на многомерные массивы, которые, похоже, не были

Я также не совсем понимаю формат второй строки. Является ли каждый капитал X только тем же, что и матрица * X в строке раньше? Что мне нужно объявить/создать для что? Или есть более простой способ, чтобы все это?

EDIT2: Вот что я имею в своей программе, по существу - это значительно сократилось, извините, если он еще слишком длинный.

#include <iostream> 
#include <Eigen/Dense> 

using namespace Eigen; 
using namespace std; 

#include <stdio.h> 
#include <stdlib.h> 
#include <conio.h> 
#include <math.h> 

typedef Matrix<double, 30, 30> Matrix30d; 

double L[30][30] ={{0}}; 
double Ax[30][30] = {{0}};   //[A] times [x]               
double At[30][30] = {{0}};   //A transpose 
double ct[30][30] = {{0}};   //c transpose                 
double x[30][30] = {{0}};   //primal solution                 
double w[30][30] = {{0}};   //omega, dual solution                 
double s[30][30] = {{0}};   //dual slack                
double u[30][30] = {{0}};   //[c]t - [A]t x [w] - [s]                 
double Atxw[30][30] = {{0}};  //A transpose times omega                
double t[30][30] = {{0}};   //RHS - [A]x[x]                
double v[30][30] = {{0}};   //mu - xij * sij                
double p[30][30] = {{0}};   //vij/xij                
double D2[30][30] = {{0}};   //diagonal of xij/sij                
double AD2[30][30] = {{0}};  //[A] x [D2]                 
double AD2xAt[30][30] = {{0}};  //[AD2] x [At]               
double uminp[30][30] = {{0}};  //[u] - [p]                
double AD2xuminp[30][30] = {{0}}; //[AD2] x [uminp]               
double z[30][30] = {{0}};   //[AD2] x [uminp] + [t]                
double Atxdw[30][30] = {{0}};  //[At] x [dw]                
double xt[30][30] = {{0}};   //x transpose                
double bt[30][30] = {{0}};  //b transpose                 
Matrix30d Inv;     //C++ style matrix for Eigen, maybe needed? 

int main(){ 

int r1;      //rows of A 
int c1;      //columns of A              
int i;      //row and column counters 
int j;                    
int k; 
double sum = 0; 
double size;     //size of square matrix being inverted [L]               
double *pointer[30][30]; 

FILE *myLPproblem;      

myLPproblem = fopen("LPproblem.txt", "r"); //Opens file and reads in data 

float firstLine[4]; 
int Anz; 

for (i = 0; i < 4; i++) 
{ 
    fscanf(myLPproblem, "%f", &firstLine[i]); 
} 

r1 = firstLine[0]; 
c1 = firstLine[1]; 
Anz = firstLine[2]; 

double A[r1][c1]; 
double b[r1][1]; 
double c[1][c1]; 
int Ap[c1+1]; 
int Ai[Anz]; 
double Ax2[Anz]; 

for(i=0; i<r1; i++){ 
    for(j=0; j<c1; j++){ 
    A[i][j]=0; 
    } 
} 

for (i = 0; i < (c1 + 1); i++) 
{ 
    fscanf(myLPproblem, "%d", &Ap[i]); 
} 

for (i = 0; i < (Anz); i++) 
{ 
    fscanf(myLPproblem, "%d", &Ai[i]); 
} 

for (i = 0; i < (Anz); i++) 
{ 
    fscanf(myLPproblem, "%lf", &Ax2[i]); 
} 

for (i = 0; i < (r1); i++) 
{ 
    fscanf(myLPproblem, "%lf", &b[i][0]); 
} 

for (i = 0; i < (c1); i++) 
{ 
    fscanf(myLPproblem, "%lf", &c[0][i]); 
} 

fclose(myLPproblem); 

int row; 
double xval; 
int Apj; 
int Apj2; 

for(j=0; j<c1; j++){ 

Apj = Ap[j]; 
Apj2 = Ap[j+1]; 

for(i=Apj; i<Apj2; i++){ 
    row = Ai[i]; 
    xval = Ax2[i]; 
    A[row][j] = xval; 
} 
} 

size = r1; 

for(i=0; i<c1; i++)      //Create c transpose                
{ 
    ct[i][0] = c[0][i]; 
} 

for(i=0; i<r1; i++)      //Create b transpose               
{ 
    bt[i][0] = b[0][i]; 
} 

for(i=0; i<c1; i++)      //Create A transpose               
    { 
    for(j=0; j<r1; j++) 
    { 
    At[i][j] = A[j][i]; 
    } 
} 

while(1){         //Main loop for iterations 

for (i = 0; i <= r1; i++) {    //Multiply [A] times [x]           
    for (j = 0; j <= 1; j++) { 
    sum = 0; 
    for (k = 0; k <= c1; k++) { 
     sum = sum + A[i][k] * x[k][j]; 
    } 
    Ax[i][j] = sum; 
    } 
} 

sum = 0;       //Multiply [At] times [w]                 

for (i = 0; i <= c1; i++){ 
    for (j = 0; j <= 1; j++) { 
    sum = 0; 
    for (k = 0; k <= r1; k++) { 
     sum = sum + At[i][k] * w[k][j]; 
    } 
    Atxw[i][j] = sum; 
    } 
} 

for(i=0; i<c1; i++)     //Subtraction to create matrix u             
{for(j=0; j<1; j++) 
    { 
    u[i][j] = (ct[i][j]) - (Atxw[i][j]) - (s[i][j]); 
    } 
} 

for(i=0; i<r1; i++)      //Subtraction to create matrix t              
{for(j=0; j<1; j++) 
    { 
    t[i][j] = (b[i][j]) - (Ax[i][j]); 
    } 
} 

for(i=0; i<c1; i++)    //Subtract and multiply to make matrix v             
{for(j=0; j<1; j++) 
    { 
    v[i][j] = mu - x[i][j]*s[i][j]; 
    } 
} 

for(i=0; i<c1; i++)    //create matrix p                
{for(j=0; j<1; j++) 
    { 
    p[i][j] = v[i][j]/x[i][j]; 
    } 
} 

for(i=0; i<c1; i++)    //create matrix D2                
{for(j=0; j<c1; j++) 
    { 
    if(i == j){ 
    D2[i][j] = x[i][0]/s[i][0]; 
    }else{ 
    D2[i][j] = 0; 
    } 
    } 
} 

sum = 0;                   

for (i = 0; i <= r1; i++) {   //Multiply [A] times [D2] 
    for (j = 0; j <= c1; j++) { 
    sum = 0; 
    for (k = 0; k <= c1; k++) { 
     sum = sum + A[i][k] * D2[k][j]; 
    } 
    AD2[i][j] = sum; 
    } 
} 

sum = 0;                   

for (i = 0; i <= r1; i++) {  //Multiply [AD2] times [At], to be inverted! 
    for (j = 0; j <= r1; j++) { 
    sum = 0; 
    for (k = 0; k <= c1; k++) { 
     sum = sum + AD2[i][k] * At[k][j]; 
    } 
    AD2xAt[i][j] = sum; 
    } 
} 

//Here is where I need to calculate the inverse (and determinant probably)  of matrix AD2xAt. I'd like to inverse to then be stored as [L]. 
//cout << "The determinant of AD2xAt is " << AD2xAt.determinant() << endl; 
//cout << "The inverse of AD2xAt is:\n" << AD2xAt.inverse() << endl; 

printf("\n\nThe inverse of AD2xAt, L, is : \n\n"); //print matrix L        

for (i=0; i<size; i++) 
{ 
    for (j=0; j<size; j++) 
    { 
     printf("%.3f\t",AD2xAt[i][j]); 
    } 
    printf("\n"); 
} 
} 

return 0; 
} 

В двух словах, он считывает матрицу из файла, вычисляет кучу матриц, то необходимо, чтобы инвертировать AD2xAt и сохранить его как L. Критическая часть находится в конце, где мне нужно взять обратный (прокрутите вниз - я прокомментировал это).

+0

показать, что вы пробовали до сих пор – Cherubim

+0

Я отредактировал мой вопрос, чтобы добавить дополнительную информацию. – user297883

+0

Нет, он говорит о вашем коде, чтобы мы могли видеть объявления и способы заполнения и ограничения массивов. В противном случае все, что может сделать это, - это предложить общие предложения, для которых не предназначен StackOverflow. –

ответ

0

Вы пытались Map<MatrixXd>(A[0],30,30).inverse() ?? –   ggael

Что вы предлагаете кажется, что это будет делать и сразу или что-то?

Правильно, Map<MatrixXd>() возвращает Эйген-х MatrixXd, на которых метод inverse() называется.

Могу ли я спросить, что [0] после A?

[0] является оператором индекс массива [] обозначающего 0 его элемента; A[0] - это начальная строка матрицы A[30][30] и преобразуется в указатель на A[0][0], соответствующий X, который вы видели.

Смежные вопросы