← Back to team overview

dolfin team mailing list archive

Re: [Question #112556]: Linear solver in C++

 

On May 28 2010, Phil Marinier wrote:

New question #112556 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/112556

I am attempting to use LinearSolver in C++ to do some simple matrix and vector solution, as in 3x3 or 5x5. I get this error.


[0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Object is in wrong state! [0]PETSC ERROR: Not for unassembled matrix! [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Petsc Release Version 3.0.0, Patch 8, Fri Aug 21 14:02:12 CDT 2009


This seems to suggest that I cannot use LinearSolver in the code directly, that instead I have to make a form file. Is this in fact the case? and if so how would I make that form file?


You need to call A.apply(), X.apply(), etc before using the matrix and vectors.

Manipulating small matrices using sparse linear algebra objects will be dead slow. I would recommend using Armadillo (take a look in dolfin/adaptivity for an example).

Garth



This is my code in case it is relevant. I am doing some simple interpolation to calculate a step size for an iterative non-linear problem. This is the test program I wrote.


int main()
{
	//declare and initialise
	int degree = 4;
	double *x = new double[degree + 1];
	double *y = new double[degree + 1];
	
	std::vector<dolfin::uint> columns;
	std::vector<double> values;
	
	dolfin::LinearSolver linSolver;
	dolfin::Matrix X(degree + 1, degree + 1);
	dolfin::Vector Y(degree + 1);
	dolfin::Vector A(degree + 1);
	
	//generate test points
	x[0] = 0;
	y[0] = 0;
	
	x[1] = 1;
	y[1] = -2;
	
	x[2] = 2;
	y[2] = -1;
	
	x[3] = 3;
	y[3] = -2;
	
	x[4] = 4;
	y[4] = 0;
	
	
	//set vector
	Y.set_local(y);
	
	//set matrix
	for (int j = 0; j <= degree; j ++)
		columns.push_back(dolfin::uint(j));
	
	for (dolfin::uint row = 0; row < dolfin::uint(degree); row++){
		
		for (int j = 0; j <= degree; j++)
			values.push_back( power(x[row], (degree - j) ) );
		
		X.setrow(row, columns, values);
		
		values.clear();
	}
	
	//solve [X]{A} = {Y} where A is the coefficients of the polynomial
	linSolver.solve(X, A, Y);
	
	//output
	for (int j = 0; j <= degree; j++)
		std::cout << A[j] << std::endl;
	
	//clean up
	delete [] x;
	delete [] y;
}


Any help is appreciated. Thank you.

You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.

_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
Unsubscribe : https://launchpad.net/~dolfin
More help   : https://help.launchpad.net/ListHelp





Follow ups

References