使用带有特征的 C 样式数组进行矩阵逆
Using C-style Arrays with Eigen for Matrix Inverse
我有大约 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]
的指针。