c++ - 使用带有 Eigen 的 C 风格数组来实现矩阵逆

Question

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

在我的程序中，我进行了一系列矩阵计算，然后需要在某个点找到矩阵的逆矩阵，我们将其称为矩阵L[30][30]。要使用 Eigen，我需要将其采用特殊的 Eigen 矩阵格式来调用函数 m.inverse，如下所示:

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

我的目标是找到一种方法...任何方法，将我的数据从 L 转换为 Eigen 可以接受的格式，这样我就可以运行这个东西。我花了过去两个小时研究这个问题，但一无所获。 :-( 我对 C 相当陌生，所以请尽可能彻底。我想要尽可能简单的方法。我读过有关映射的内容，但遗憾的是我对指针不是很清楚(这似乎是一个组成部分)。有没有一种方法可以循环遍历每一行和每一列并将它们复制到特征矩阵中？

当我问时，我是否需要获取生成的特征矩阵并将其转回 C 数组？这个过程将如何进行？预先感谢您的任何帮助!我花了大约 50-60 个小时在这上面，本周到期!这是我需要做的最后一件事，我将完成我的学期项目。这是一门数学课，所以编程方面的事情对我来说有点模糊，但我学到了很多东西。

可能相关的信息:-运行于 Windows 10 i7 处理器 Sony VAIO- 使用 C++ 中的 CodeBlocks 进行编译，但最初是用 C 编写的-这段代码全部在一个 while 循环中，可能会迭代 10 次左右。-这个矩阵L每次迭代都需要计算矩阵逆，并且每次的数据都会不同。

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

编辑 - 我看到了这个并尝试实现它但没有成功，但如果我能弄清楚这一点，这似乎是解决方案:

“假设您有一个具有大小为 nRows x nCols 的 double 值的数组。

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

您可以使用 map 功能创建一个 nRows x nCols 大小的双矩阵，如下所示:

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

Map 操作将现有的内存区域映射到 Eigen 的数据结构中。像这样的单行可以避免编写丑陋的矩阵创建代码、for 循环并按良好顺序复制每个元素等。”

这似乎是一个不错的解决方案，但我对如何使用“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。关键部分在最后，我需要在其中进行逆运算(滚动到底部 - 我已对其进行了评论)。

Answer 1

Have you tried Map<MatrixXd>(A[0],30,30).inverse() ?? – ggael

What you're proposing seems like it would be doing both at once or something?

对了，Map<MatrixXd>()返回 Eigen 的 MatrixXd ，其中方法 inverse()被调用。

May I ask what the [0] is after A?

[0]是数组下标运算符 []指定0 -第一个元素； A[0]是矩阵的第一行A[30][30]并转换为指向 A[0][0] 的指针对应X你看到了。