Subversion Repositories slepc-dev

/*

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain

   This file is part of SLEPc.

   SLEPc is free software: you can redistribute it and/or modify it under  the

   terms of version 3 of the GNU Lesser General Public License as published by

   the Free Software Foundation.

   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY

   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS

   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for

   more details.

   You  should have received a copy of the GNU Lesser General  Public  License

   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*/

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

#include "petscblaslapack.h"

#include <stdlib.h>

PetscLogEvent SLEPC_UpdateVectors = 0;

#undef __FUNCT__  

#define __FUNCT__ "SlepcVecSetRandom"

/*@

   SlepcVecSetRandom - Sets all components of a vector to random numbers which

   follow a uniform distribution in [0,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;

  PetscInt       i,n,low,high;

  PetscScalar    *px,t;

#if defined(PETSC_HAVE_DRAND48)

  static unsigned short seed[3] = { 1, 3, 2 };

#endif

  PetscFunctionBegin;

  ierr = VecGetSize(x,&n);CHKERRQ(ierr);

  ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr);

  ierr = VecGetArray(x,&px);CHKERRQ(ierr);

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

#if defined(PETSC_HAVE_DRAND48)

    t = erand48(seed);

#elif defined(PETSC_HAVE_RAND)

    t = rand()/(PetscReal)((unsigned int)RAND_MAX+1);

#else

    t = 0.5;

#endif

    if (i>=low && i<high) px[i-low] = t;

  }

  ierr = VecRestoreArray(x,&px);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 matrix formats MATSEQSBAIJ or MATMPISBAIJ,

   or when PETSc is configured with complex scalars.

   Level: developer

@*/

PetscErrorCode SlepcIsHermitian(Mat A,PetscTruth *is)

{

  PetscErrorCode ierr;

  PetscInt       M,N,m,n;

  Vec            x,w1,w2;

  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);

  ierr = VecAXPY(w2,-1.0,w1);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"

/*@C

   SlepcAbsEigenvalue - Returns the absolute value of a complex number given

   its real and imaginary parts.

   Not collective

   Input parameters:

+  x  - the real part of the complex number

-  y  - the imaginary part of the complex number

   Notes:

   This function computes sqrt(x**2+y**2), taking care not to cause unnecessary

   overflow. It is based on LAPACK's DLAPY2.

   Level: developer

@*/

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

#undef __FUNCT__  

#define __FUNCT__ "SlepcMatConvertSeqDense"

/*@C

   SlepcMatConvertSeqDense - Converts a parallel matrix to another one in sequential

   dense format replicating the values in every processor.

   Collective

   Input parameters:

+  A  - the source matrix

-  B  - the target matrix

   Level: developer

@*/

PetscErrorCode SlepcMatConvertSeqDense(Mat mat,Mat *newmat)

{

  PetscErrorCode ierr;

  PetscInt       m,n;

  PetscMPIInt    size;

  MPI_Comm       comm;

  Mat            *M;

  IS             isrow, iscol;

  PetscTruth     flg;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(mat,MAT_COOKIE,1);

  PetscValidPointer(newmat,2);

  ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr);

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

  if (size > 1) {

    /* assemble full matrix on every processor */

    ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr);

    ierr = ISCreateStride(PETSC_COMM_SELF,m,0,1,&isrow);CHKERRQ(ierr);

    ierr = ISCreateStride(PETSC_COMM_SELF,n,0,1,&iscol);CHKERRQ(ierr);

    ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&M);CHKERRQ(ierr);

    ierr = ISDestroy(isrow);CHKERRQ(ierr);

    ierr = ISDestroy(iscol);CHKERRQ(ierr);

    /* Fake support for "inplace" convert */

    if (*newmat == mat) {

      ierr = MatDestroy(mat);CHKERRQ(ierr);

    }

    *newmat = *M;

    ierr = PetscFree(M);CHKERRQ(ierr);    

    /* convert matrix to MatSeqDense */

    ierr = PetscTypeCompare((PetscObject)*newmat,MATSEQDENSE,&flg); CHKERRQ(ierr);

    if (!flg) {

      ierr = MatConvert(*newmat,MATSEQDENSE,MAT_INITIAL_MATRIX,newmat);CHKERRQ(ierr);

    }

  } else {

    /* convert matrix to MatSeqDense */

    ierr = MatConvert(mat,MATSEQDENSE,MAT_INITIAL_MATRIX,newmat);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);  

}

#undef __FUNCT__  

#define __FUNCT__ "SlepcCheckOrthogonality"

/*@

   SlepcCheckOrthogonality - Checks (or prints) the level of orthogonality

   of a set of vectors.

   Collective on Vec

   Input parameters:

+  V  - a set of vectors

.  nv - number of V vectors

.  W  - an alternative set of vectors (optional)

.  nw - number of W vectors

-  B  - matrix defining the inner product (optional)

   Output parameter:

.  lev - level of orthogonality (optional)

   Notes:

   This function computes W'*V and prints the result. It is intended to check

   the level of bi-orthogonality of the vectors in the two sets. If W is equal

   to PETSC_NULL then V is used, thus checking the orthogonality of the V vectors.

   If matrix B is provided then the check uses the B-inner product, W'*B*V.

   If lev is not PETSC_NULL, it will contain the level of orthogonality

   computed as ||W'*V - I|| in the Frobenius norm. Otherwise, the matrix W'*V

   is printed.

   Level: developer

@*/

PetscErrorCode SlepcCheckOrthogonality(Vec *V,PetscInt nv,Vec *W,PetscInt nw,Mat B,PetscScalar *lev)

{

  PetscErrorCode ierr;

  PetscInt       i,j;

  PetscScalar    *vals;

  Vec            w;

  MPI_Comm       comm;

  PetscFunctionBegin;

  if (nv<=0 || nw<=0) PetscFunctionReturn(0);

  ierr = PetscObjectGetComm((PetscObject)V[0],&comm);CHKERRQ(ierr);

  ierr = PetscMalloc(nv*sizeof(PetscScalar),&vals);CHKERRQ(ierr);

  if (B) { ierr = VecDuplicate(V[0],&w);CHKERRQ(ierr); }

  if (lev) *lev = 0.0;

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

    if (B) {

      if (W) { ierr = MatMultTranspose(B,W[i],w);CHKERRQ(ierr); }

      else { ierr = MatMultTranspose(B,V[i],w);CHKERRQ(ierr); }

    }

    else {

      if (W) w = W[i];

      else w = V[i];

    }

    ierr = VecMDot(w,nv,V,vals);CHKERRQ(ierr);

    for (j=0;j<nv;j++) {

      if (lev) *lev += (j==i)? (vals[j]-1.0)*(vals[j]-1.0): vals[j]*vals[j];

      else {

#ifndef PETSC_USE_COMPLEX

        ierr = PetscPrintf(comm," %12g  ",vals[j]);CHKERRQ(ierr);

#else

        ierr = PetscPrintf(comm," %12g%+12gi ",PetscRealPart(vals[j]),PetscImaginaryPart(vals[j]));CHKERRQ(ierr);    

#endif

      }

    }

    if (!lev) { ierr = PetscPrintf(comm,"\n");CHKERRQ(ierr); }

  }

  ierr = PetscFree(vals);CHKERRQ(ierr);

  if (B) { ierr = VecDestroy(w);CHKERRQ(ierr); }

  if (lev) *lev = PetscSqrtScalar(*lev);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "SlepcUpdateVectors"

/*@

   SlepcUpdateVectors - Update a set of vectors V as V(:,s:e-1) = V*Q(:,s:e-1).

   Collective on Vec

   Input parameters:

+  n      - number of vectors in V

.  s      - first column of V to be overwritten

.  e      - first column of V not to be overwritten

.  Q      - matrix containing the coefficients of the update

.  ldq    - leading dimension of Q

-  qtrans - flag indicating if Q is to be transposed

   Input/Output parameter:

.  V      - set of vectors

   Notes:

   This function computes V(:,s:e-1) = V*Q(:,s:e-1), that is, given a set of

   vectors V, columns from s to e-1 are overwritten with columns from s to

   e-1 of the matrix-matrix product V*Q.

   Matrix V is represented as an array of Vec, whereas Q is represented as

   a column-major dense array of leading dimension ldq. Only columns s to e-1

   of Q are referenced.

   If qtrans=PETSC_TRUE, the operation is V*Q'.

   Level: developer

@*/

PetscErrorCode SlepcUpdateVectors(PetscInt n_,Vec *V,PetscInt s,PetscInt e,const PetscScalar *Q,PetscInt ldq_,PetscTruth qtrans)

{

  PetscErrorCode ierr;

  PetscInt       l;

  PetscBLASInt   i,j,k,bs=64,m,n,ldq,ls;

  PetscScalar    *pv,*pw,*pq,*work,*pwork,one=1.0,zero=0.0;

  const char     *qt;

  PetscFunctionBegin;

  n = PetscBLASIntCast(n_);

  ldq = PetscBLASIntCast(ldq_);

  m = e-s;

  if (m==0) PetscFunctionReturn(0);

  PetscValidIntPointer(Q,5);

  if (m<0 || n<0 || s<0 || e>n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Index argument out of range");

  ierr = PetscLogEventBegin(SLEPC_UpdateVectors,0,0,0,0);CHKERRQ(ierr);

  ierr = VecGetLocalSize(V[0],&l);CHKERRQ(ierr);

  ls = PetscBLASIntCast(l);

  ierr = VecGetArray(V[0],&pv);CHKERRQ(ierr);

  if (qtrans) {

    pq = (PetscScalar*)Q+s;

    qt = "T";

  } else {

    pq = (PetscScalar*)Q+s*ldq;

    qt = "N";

  }

  ierr = PetscMalloc(sizeof(PetscScalar)*bs*m,&work);CHKERRQ(ierr);

  k = ls % bs;

  if (k) {

    BLASgemm_("N",qt,&k,&m,&n,&one,pv,&ls,pq,&ldq,&zero,work,&k);

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

      pw = pv+(s+j)*ls;

      pwork = work+j*k;

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

        *pw++ = *pwork++;

      }

    }        

  }

  for (;k<ls;k+=bs) {

    BLASgemm_("N",qt,&bs,&m,&n,&one,pv+k,&ls,pq,&ldq,&zero,work,&bs);

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

      pw = pv+(s+j)*ls+k;

      pwork = work+j*bs;

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

        *pw++ = *pwork++;

      }

    }

  }

  ierr = VecRestoreArray(V[0],&pv);CHKERRQ(ierr);

  ierr = PetscFree(work);CHKERRQ(ierr);

  ierr = PetscLogFlops(m*n*2*ls);CHKERRQ(ierr);

  ierr = PetscLogEventEnd(SLEPC_UpdateVectors,0,0,0,0);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}