Lineares Gleichungssystem Lösen mit lapack oder Eigen



  • Hi,

    ich habe ein Problem beim Lösen von Gleichungssystemen. Ich habe mir ein kleines Beispiel zum Lösen eines linearen Gleichungssystems mit lapack und der Eigen-Bibliothek geschrieben. Leider kommen beide Verfahren auf unterschiedliche Ergebnisse wobei Eigen das richtige zu bestimmen scheint. Kann mir jemand sagen was ich bei der Verwendung von lapack falsch mache? Lapack schreibt die Lösung des Gleichungssystems in das Array, welches voher die RHS gebildet hat, mit ist aber zudem aufgefallen, dass die Koeffizientenmatrix nach dem Lösen mit Lapack auch verändert wurde.
    Hier mein Beispiel. Eigentlich sollte da doch das gleich raus kommen!?!

    int size = 2;
    
    	VectorXd eigenB = VectorXd::Random(size);
    	ublas::matrix<double,  ublas::column_major>  boostB(eigenB.rows(),1);
    
    	for(unsigned i = 0; i<eigenB.rows(); i++)
    	{
    		boostB(i,0) = eigenB(i);
    	}
    
    //	cout<<"eigen RHS:"<<endl;
    //	cout<<eigenB<<endl;
    //
    //	cout<<"Boost RHS:"<<endl;
    //	cout<<boostB<<endl;
    
    	MatrixXd eigenMatrix = MatrixXd::Random(size,size);
    	ublas::matrix<double> boostMatrix(size,size);
    
    	for(unsigned i = 0; i < eigenMatrix.rows(); i++)
    	{
    		for(unsigned j = 0; j < eigenMatrix.cols(); j++)
    		{
    			boostMatrix(i,j) = eigenMatrix(i,j);
    		}
    	}
    
    //		cout<<"eigen Matrix:"<<endl;
    //		cout<<eigenMatrix<<endl;
    //
    //
    //			cout<<"Boost Matrix:"<<endl;
    //			for(unsigned i = 0; i<boostMatrix.size1(); i++)
    //			{
    //				for(unsigned j = 0; j<boostMatrix.size2(); j++)
    //				{
    //					cout<<boostMatrix(i,j)<<"  ";
    //				}
    //				cout<<endl;
    //			}
    
    	VectorXd x = eigenMatrix.fullPivLu().solve(eigenB);
    	int info = lapack::gesv(boostMatrix, boostB);
    
    // LU factorisation
    //	std::vector<int> ipiv(boostMatrix.size1());
    //	int test = lapack::getrf(boostMatrix, ipiv);
    //
    //	if(test != 0) {
    //		std::stringstream str;
    //		str << "singular system matrix! (info == " << test << ")";
    //		throw std::runtime_error(str.str());
    //	}
    //
    //	//	// solve using LU
    //		lapack::getrs('N', boostMatrix, ipiv, boostB);
    
    	cout<<"eigen result"<<endl;
    	cout<<x<<endl;
    	cout<<"------ Boost Results"<<endl;
    	cout<<boostB<<endl;
    
    //	cout<<"eigen Matrix:"<<endl;
    //	cout<<eigenMatrix<<endl;
    //
    //
    //		cout<<"Boost Matrix:"<<endl;
    //		for(unsigned i = 0; i<boostMatrix.size1(); i++)
    //		{
    //			for(unsigned j = 0; j<boostMatrix.size2(); j++)
    //			{
    //				cout<<boostMatrix(i,j)<<"  ";
    //			}
    //			cout<<endl;
    //		}
    


  • okay, hätte ich auch schneller drauf kommen können. Row-major von "Eigen" und Column-major von lapack muss beachtet werden. Ne, also so ganz verstehe ich es doch nicht. Wenn ich die BoostMatrix transponiere, bevor ich das Gleichungssystem löse, dann erhalte ich die richtige Lösung. Row major oder column major ist doch nur die Reihenfolge, in der die Einträge der Matrizen gespeichert werden. Egel welche Reihenfolge verwendet wird, es sollten doch beide auf das selbe Ergebnis kommen.


Log in to reply