Subversion Repositories slepc-dev

/*                      

       This file implements a wrapper to the LAPACK eigenvalue subroutines.

       Currently, only LAPACK routines for standard problems are used.

       Generalized problems are transformed to standard ones.

*/

#include "src/eps/impls/lapack/lapackp.h"

#include "slepcblaslapack.h"

#undef __FUNCT__  

#define __FUNCT__ "EPSSetUp_LAPACK"

static int EPSSetUp_LAPACK(EPS eps)

{

  int         ierr,i,size,rank,n,m,row,nz,*cols,dummy;

  PetscScalar *vals;

  EPS_LAPACK *la = (EPS_LAPACK *)eps->data;

  MPI_Comm    comm = eps->comm;

  PetscFunctionBegin;

  ierr = EPSComputeExplicitOperator(eps,&la->BA);CHKERRQ(ierr);

  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);

  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);

  ierr     = MatGetSize(la->BA,&n,&n);CHKERRQ(ierr);

  if (size > 1) { /* assemble matrix on first processor */

    if (!rank) {

      ierr = MatCreateMPIDense(comm,n,n,n,n,PETSC_NULL,&la->A);CHKERRQ(ierr);

    }

    else {

      ierr = MatCreateMPIDense(comm,0,n,n,n,PETSC_NULL,&la->A);CHKERRQ(ierr);

    }

    PetscLogObjectParent(la->BA,la->A);

    ierr = MatGetOwnershipRange(la->BA,&row,&dummy);CHKERRQ(ierr);

    ierr = MatGetLocalSize(la->BA,&m,&dummy);CHKERRQ(ierr);

    for (i=0; i<m; i++) {

      ierr = MatGetRow(la->BA,row,&nz,&cols,&vals);CHKERRQ(ierr);

      ierr = MatSetValues(la->A,1,&row,nz,cols,vals,INSERT_VALUES);CHKERRQ(ierr);

      ierr = MatRestoreRow(la->BA,row,&nz,&cols,&vals);CHKERRQ(ierr);

      row++;

    }

    ierr = MatAssemblyBegin(la->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    ierr = MatAssemblyEnd(la->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSSetDefaults_LAPACK"

static int EPSSetDefaults_LAPACK(EPS eps)

{

  int         ierr, N;

  PetscFunctionBegin;

  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);

  if (eps->ncv) {

    if (eps->ncv<1 || eps->ncv>N) SETERRQ(1,"Wrong value of ncv");

  }

  else eps->ncv = eps->nev;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSSolve_LAPACK"

static int  EPSSolve_LAPACK(EPS eps,int *its)

{

  int         ierr,n,size,rank;

  PetscScalar *array,*pV;

  EPS_LAPACK *la = (EPS_LAPACK *)eps->data;

  MPI_Comm    comm = eps->comm;

  PetscFunctionBegin;

  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);

  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);

  if (size>1) {

    ierr = MatGetArray(la->A,&array);CHKERRQ(ierr);

  } else {

    ierr = MatGetArray(la->BA,&array);CHKERRQ(ierr);

  }

  ierr = MatGetSize(la->BA,&n,&n);CHKERRQ(ierr);

  if (eps->dropvectors) pV = PETSC_NULL;

  else {

    ierr = VecGetArray(eps->V[0],&pV);CHKERRQ(ierr);

  }

  ierr = EPSDenseNHEPSorted(n,array,eps->eigr,eps->eigi,pV,eps->ncv,eps->which);CHKERRQ(ierr);

  if (!eps->dropvectors) {

    ierr = VecRestoreArray(eps->V[0],&pV);CHKERRQ(ierr);

  }

  if (size > 1) {

    ierr = MatRestoreArray(la->A,&array);CHKERRQ(ierr);

  } else {

    ierr = MatRestoreArray(la->BA,&array);CHKERRQ(ierr);

  }

  eps->nconv = eps->ncv;

  eps->its   = 1;

  *its       = eps->its;

  eps->reason = EPS_CONVERGED_TOL;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSDestroy_LAPACK"

/*

  EPSDestroy_LAPACK - Destroys the context variable for LAPACK.

  Input Parameter:

. eps - the iterative context

*/

int EPSDestroy_LAPACK(EPS eps)

{

  int         ierr,size;

  EPS_LAPACK *la = (EPS_LAPACK *)eps->data;

  MPI_Comm    comm = eps->comm;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE);

  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);

  if (size > 1) {

    ierr = MatDestroy(la->A);CHKERRQ(ierr);

  }

  ierr = MatDestroy(la->BA);CHKERRQ(ierr);

  ierr = EPSDefaultDestroy(eps);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

EXTERN_C_BEGIN

#undef __FUNCT__  

#define __FUNCT__ "EPSCreate_LAPACK"

int EPSCreate_LAPACK(EPS eps)

{

  EPS_LAPACK *la;

  int ierr;

  PetscFunctionBegin;

  ierr = PetscNew(EPS_LAPACK,&la);CHKERRQ(ierr);

  PetscMemzero(la,sizeof(EPS_LAPACK));

  PetscLogObjectMemory(eps,sizeof(EPS_LAPACK));

  eps->data                      = (void *) la;

  eps->ops->setup                = EPSSetUp_LAPACK;

  eps->ops->setdefaults          = EPSSetDefaults_LAPACK;

  eps->ops->solve                = EPSSolve_LAPACK;

  eps->ops->destroy              = EPSDestroy_LAPACK;

  eps->ops->view                 = 0;

  PetscFunctionReturn(0);

}

EXTERN_C_END