[petsc-users] algorithm to get pure real eigen valule for general eigenvalue problem(Non-hermitian type)

[Jose E. Roman]

El 15/04/2013, a las 17:25, Zhang Wei escribió:

> Hi, thanks for your reply!
> Here is my code:
>    //  Compute the operator matrix that defines the eigensystem, Ax=kx
>    ierr = MatCreateShell(PETSC_COMM_WORLD,n,n,PETSC_DETERMINE,PETSC_DETERMINE,NULL,&A);CHKERRQ(ierr);
>    ierr = MatSetFromOptions(A);CHKERRQ(ierr);
> 
>    ierr = MatShellSetContext(A, &pa); CHKERRQ(ierr);
>    ierr = MatShellSetOperation(A, MATOP_MULT, (void(*)()) MatVecMult); CHKERRQ(ierr);
> 
>    ierr = MatGetVecs(A, PETSC_NULL,&xr); CHKERRQ(ierr);
>    ierr = MatGetVecs(A, PETSC_NULL,&xi); CHKERRQ(ierr);
>    ierr = MatGetVecs(A, PETSC_NULL,&iv); CHKERRQ(ierr);
> 
>    //  Create the eigensolver and set various options
>    ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
> 
>    //  Set operators. In this case, it is a standard eigenvalue problem
>    ierr = EPSSetOperators(eps, A, PETSC_NULL);CHKERRQ(ierr);
>    ierr = EPSSetProblemType(eps, EPS_NHEP);CHKERRQ(ierr);
>    ierr = EPSSetDimensions(eps, nev, ncv, mpd);CHKERRQ(ierr);
>    ierr = EPSSetWhichEigenpairs(eps,which);CHKERRQ(ierr);
>    ierr = EPSSetType(eps,method.c_str());CHKERRQ(ierr);
>    ierr = EPSSetTolerances(eps,s_tol,s_maxit);CHKERRQ(ierr);
>    ierr = EPSSetConvergenceTest(eps,cov);CHKERRQ(ierr);
>    ierr = EPSKrylovSchurSetRestart(eps,keep);CHKERRQ(ierr);
>    //  Set solver parameters at runtime
>    ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
> 
>    //  Set initial vector
>  ierr = pa.setInitialVector(iv);
>  ierr = EPSSetInitialSpace(eps,1,&iv);CHKERRQ(ierr);
> 
>    ierr = EPSView(eps, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
> 
>    //  Solve the eigensystem
>    ierr = EPSSolve(eps);CHKERRQ(ierr);
> 
> and pa is a class, where I define the OP,  MatVecMult. It actually looks like this:
> 
> PetscErrorCode MatVecMult(Mat A, Vec x, Vec y)
> {
>    linearizedTimeStepper *pa;
>    MatShellGetContext(A, &pa);
>    return pa->MatVecMult(A, x, y);
> }

This code is basically correct. What I was asking is whether you had tested the matrix (that is, pa->MatVecMult) in other contexts before computing eigenvalues. 

Jose