2015-12-12 97 views
0

我有大约1000行代码,我用C编写了线性规划求解器(内点算法)。我意识到我需要使用Eigen来计算矩阵逆,所以现在我用C++运行我的代码(运行得很好,看起来好像)。现在我有一堆C画幅声明数组,例如:A[30][30];使用C型阵列和特征矩阵求逆

在我的节目,我做了一堆矩阵计算的,然后需要找到在某个点矩阵的逆,我们称之为矩阵L[30][30] 。要使用本征,我需要把它在一个特殊的特征矩阵格式来调用函数m.inverse这样的:

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

我的目标是找到一种方法...任何方式,从中获取我的数据L以Eigen将接受的格式,所以我可以运行这个东西。我花了最近2个小时研究这个问题,并没有提出任何建议。 :-(我对C相当陌生,所以请尽量详尽,我想要最简单的方法,我已经阅读了映射,但是我不太清楚指针的含义(这似乎是是否有一种方法可以循环遍历每一行和每一列并将它们复制到特征矩阵中?

虽然我在问,我是否需要将生成的特征矩阵转化为一个C数组?这个过程如何工作?提前感谢任何帮助!我已经花了大约50-60小时的时间,这是本周到期的!这是我需要做的最后一件事,我将与我一起完成这是一个数学课程,所以编程方面对我来说有点模糊,但是我学到了很多东西,

可能相关信息: 在Windows 10 i7处理器上运行索尼VAIO - 在C++中编译CodeBlocks,但最初编写为C - 此代码全部在while循环中,可能会迭代10次左右。 - 每次迭代需要为矩阵L计算矩阵逆,并且数据每次都会不同。

请帮忙!我愿意学习,但我需要指导,而且这门课是在线的,所以我几乎没有。非常感谢!

编辑 - 我看到了这一点,并试图将其落实无济于事,但它似乎是解决办法,如果我能想出解决办法:

“假设你有大小NROWS X NCOLS的双值的数组。

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

您可以使用这样的地图功能创建一个NROWS X NCOLS尺寸双矩阵:

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

地图操作映射现有内存REG离子进入特征的数据结构。像这样的单行允许避免写出难看的矩阵创建代码,用于以良好顺序复制每个元素的循环等。“

这似乎是一个很好的解决方案,但我对如何做任何事情都无能为力“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]是一个后?

[0]是指定0个元素的数组下标操作者[]; A[0]是矩阵A[30][30]的第一行,并转换为指向A[0][0]的指针,对应于您看到的X