使用带有特征的 C 样式数组进行矩阵逆

Using C-style Arrays with Eigen for Matrix Inverse

本文关键字:数组 样式 特征      更新时间:2023-10-16

我有大约 1000 行代码是用 C 语言编写的,用于线性规划求解器(内部点算法(。 我意识到我需要使用 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 将接受的格式,以便我可以运行这个东西。 我花了过去 2 个小时研究这个,但一无所获。 :-( 我对 C 语言相当陌生,所以请尽可能彻底。 我想要最简单的方法。 我已经读过有关映射的信息,但遗憾的是,我对指针不是很清楚(这似乎是一个不可或缺的一部分(。 有没有办法遍历每一行和每一列并将它们复制到特征矩阵中?

当我问时,我是否需要获取生成的特征矩阵并将其重新转换为 C 数组? 这个过程将如何运作? 提前感谢任何帮助! 我花了大约 50-60 个小时,本周到期! 这是我需要做的最后一件事,我将完成我的学期项目。 这是一堂数学课,所以编程方面的东西对我来说有点模糊,但我学到了很多东西。

可能相关信息:- 运行在视窗 10 i7 处理器索尼 VAIO- 用 C++ 种代码块编译,但最初是用 C 语言编写的-这段代码都在一个while循环中,可以迭代10次左右。-每次迭代都需要为这个矩阵L计算矩阵逆,每次数据都会不同。

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

编辑 - 我看到了这个并试图实现它无济于事,但如果我能弄清楚这一点,这似乎是解决方案:

"假设你有一个数组,其大小为 nRow x nCols 的双精度值。

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

您可以使用映射功能创建 nRows x nCols 大小的双精度矩阵,如下所示:

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

映射操作将现有内存区域映射到特征的数据结构中。像这样的一行可以避免编写矩阵创建的丑陋代码,用于循环,以良好的顺序复制每个元素等。

这似乎是一个很好的解决方案,但我对如何处理"指向一些数据"的"双 * 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("nnThe inverse of AD2xAt, L, is : nn");   //print matrix L                               
 for (i=0; i<size; i++)
 {
     for (j=0; j<size; j++)
     {
         printf("%.3ft",AD2xAt[i][j]);
     }
     printf("n");
 }
}
return 0;
}

简而言之,它从文件中读取矩阵,计算一堆矩阵,然后需要反转AD2xAt并将其存储为L。 关键部分在最后,我需要采取相反的方法(滚动到底部 - 我有评论(。

你试过吗 Map<MatrixXd>(A[0],30,30).inverse()?? - 加尔

您的建议似乎可以同时执行或 东西?

对,Map<MatrixXd>()返回特征MatrixXd,该方法inverse()被调用。

请问 A 后面的 [0] 是什么?

[0] 是指定0 -th 元素[]数组下标运算符; A[0] 是矩阵A[30][30]的初始行,并转换为指向与您看到的X对应的A[0][0]的指针。