Subversion Repositories slepc-dev

#include "slepc.h" /*I "slepc.h" I*/

extern PetscRandom rctx;

#undef __FUNCT__  
#define __FUNCT__ "SlepcVecSetRandom"
/*@
   SlepcVecSetRandom - Sets all components of a vector to random numbers which
   follow a uniform distribution in [-1,1].

   Collective on Vec

   Input/Output Parameter:
.  x  - the vector

   Note:
   This operation is equivalent to VecSetRandom - the difference is that the
   vector generated by SlepcVecSetRandom is the same irrespective of the size
   of the communicator.

   Level: developer

.seealso: VecSetRandom()
@*/
PetscErrorCode SlepcVecSetRandom(Vec x)
{
  PetscErrorCode ierr;
  int            i,n,low,high;
  PetscScalar    *pv,*px;
  Vec            v;

  PetscFunctionBegin;
  ierr = VecGetSize(x,&n);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&v);CHKERRQ(ierr);
  ierr = VecSetRandom(rctx,v);CHKERRQ(ierr);
  ierr = VecGetArray(v,&pv);CHKERRQ(ierr);
  ierr = VecGetArray(x,&px);CHKERRQ(ierr);
  for (i=low;i<high;i++) px[i-low]=pv[i];
  ierr = VecRestoreArray(v,&pv);CHKERRQ(ierr);
  ierr = VecRestoreArray(x,&px);CHKERRQ(ierr);
  ierr = VecDestroy(v);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "SlepcIsHermitian"
/*@
   SlepcIsHermitian - Checks if a matrix is Hermitian or not.

   Collective on Mat

   Input parameter:
.  A  - the matrix

   Output parameter:
.  is  - flag indicating if the matrix is Hermitian

   Notes:
   The result of Ax and A^Hx (with a random x) is compared, but they
   could be equal also for some non-Hermitian matrices.

   This routine will not work with BOPT=O_complex and matrix formats
   MATSEQSBAIJ or MATMPISBAIJ.
   
   Level: developer

@*/
PetscErrorCode SlepcIsHermitian(Mat A,PetscTruth *is)
{
  PetscErrorCode ierr;
  int            M,N,m,n;
  Vec            x,w1,w2;
  PetscScalar    alpha;
  MPI_Comm       comm;
  PetscReal      norm;
  PetscTruth     has;

  PetscFunctionBegin;

#if !defined(PETSC_USE_COMPLEX)
  ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,is);CHKERRQ(ierr);
  if (*is) PetscFunctionReturn(0);
  ierr = PetscTypeCompare((PetscObject)A,MATMPISBAIJ,is);CHKERRQ(ierr);
  if (*is) PetscFunctionReturn(0);
#endif

  *is = PETSC_FALSE;
  ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
  ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
  if (M!=N) PetscFunctionReturn(0);
  ierr = MatHasOperation(A,MATOP_MULT,&has);CHKERRQ(ierr);
  if (!has) PetscFunctionReturn(0);
  ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&has);CHKERRQ(ierr);
  if (!has) PetscFunctionReturn(0);

  ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
  ierr = VecCreate(comm,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,n,N);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = SlepcVecSetRandom(x);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&w1);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&w2);CHKERRQ(ierr);
  ierr = MatMult(A,x,w1);CHKERRQ(ierr);
  ierr = MatMultTranspose(A,x,w2);CHKERRQ(ierr);
  ierr = VecConjugate(w2);CHKERRQ(ierr);
  alpha = -1.0;
  ierr = VecAXPY(&alpha,w1,w2);CHKERRQ(ierr);
  ierr = VecNorm(w2,NORM_2,&norm);CHKERRQ(ierr);
  if (norm<1.0e-6) *is = PETSC_TRUE;
  ierr = VecDestroy(x);CHKERRQ(ierr);
  ierr = VecDestroy(w1);CHKERRQ(ierr);
  ierr = VecDestroy(w2);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#if !defined(PETSC_USE_COMPLEX)

#undef __FUNCT__  
#define __FUNCT__ "SlepcAbsEigenvalue"
PetscReal SlepcAbsEigenvalue(PetscScalar x,PetscScalar y)
{
  PetscReal xabs,yabs,w,z,t;
  PetscFunctionBegin;
  xabs = PetscAbsReal(x);
  yabs = PetscAbsReal(y);
  w = PetscMax(xabs,yabs);
  z = PetscMin(xabs,yabs);
  if (z == 0.0) PetscFunctionReturn(w);
  t = z/w;
  PetscFunctionReturn(w*sqrt(1.0+t*t));  
}

#endif