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

[Zhang Wei]

Hi Jose.
In pa-> MatVecMult(Mat A, Vec x , Vec y). I do three things
1) copy the  x as  initial value for my fvm code:

const PetscScalar *px;
    PetscScalar       *py;
    PetscErrorCode    ierr;

    PetscFunctionBegin;
    ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr);
    ierr = VecGetArray(y,&py);CHKERRQ(ierr);

    if(loopID>0){
      copySLEPcDataToField(&px[0]);
    } 
2) In  fvm solver, where I implement first order eular method, fourth order RK method and several others.
3) after several timestep. I do sent the field to slepc:
   
   copyFieldToSLEPcData(&py[0]);
    VecAssemblyBegin(y);
    VecAssemblyEnd(y);
 
>From the those files I saved during the time stepping (all those methods), the  fvm solver works fine for me. They gives expected things.
Actually the numerical methods and accuracy in pa-> MatVecMult  doesn't affect the results that much. At least it the shape of eigen vectors do not change with those factors.

-----Original Message-----
From: petsc-users-bounces at mcs.anl.gov [mailto:petsc-users-bounces at mcs.anl.gov] On Behalf Of Jose E. Roman
Sent: den 15 april 2013 17:35
To: PETSc users list
Subject: Re: [petsc-users] algorithm to get pure real eigen valule for general eigenvalue problem(Non-hermitian type)


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