У меня есть около 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. Критическая часть находится в конце, где мне нужно взять обратный (прокрутите вниз - я прокомментировал это).
показать, что вы пробовали до сих пор – Cherubim
Я отредактировал мой вопрос, чтобы добавить дополнительную информацию. – user297883
Нет, он говорит о вашем коде, чтобы мы могли видеть объявления и способы заполнения и ограничения массивов. В противном случае все, что может сделать это, - это предложить общие предложения, для которых не предназначен StackOverflow. –