Subversion Repositories slepc-dev

/*

     Dot product routines

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

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2010, 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 <private/ipimpl.h>      /*I "slepcip.h" I*/

#undef __FUNCT__  

#define __FUNCT__ "IPNorm"

/*@

   IPNorm - Computes the norm of a vector as the square root of the inner

   product (x,x) as defined by IPInnerProduct().

   Collective on IP and Vec

   Input Parameters:

+  ip - the inner product context

-  x  - input vector

   Output Parameter:

.  norm - the computed norm

   Notes:

   This function will usually compute the 2-norm of a vector, ||x||_2. But

   this behaviour may be different if using a non-standard inner product changed

   via IPSetBilinearForm(). For example, if using the B-inner product for

   positive definite B, (x,y)_B=y^H Bx, then the computed norm is ||x||_B =

   sqrt( x^H Bx ).

   Level: developer

.seealso: IPInnerProduct()

@*/

PetscErrorCode IPNorm(IP ip,Vec x,PetscReal *norm)

{

  PetscErrorCode ierr;

  PetscScalar    p;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,2);

  PetscValidPointer(norm,3);

  if (!ip->matrix && ip->bilinear_form == IP_INNER_HERMITIAN) {

    ierr = VecNorm(x,NORM_2,norm);CHKERRQ(ierr);

  } else {

    ierr = IPInnerProduct(ip,x,x,&p);CHKERRQ(ierr);

    if (PetscAbsScalar(p)<PETSC_MACHINE_EPSILON)

      PetscInfo(ip,"Zero norm, either the vector is zero or a semi-inner product is being used\n");

#if defined(PETSC_USE_COMPLEX)

    if (PetscRealPart(p)<0.0 || PetscAbsReal(PetscImaginaryPart(p))>PETSC_MACHINE_EPSILON)

       SETERRQ(((PetscObject)ip)->comm,1,"IPNorm: The inner product is not well defined");

    *norm = PetscSqrtScalar(PetscRealPart(p));

#else

    if (p<0.0) SETERRQ(((PetscObject)ip)->comm,1,"IPNorm: The inner product is not well defined");

    *norm = PetscSqrtScalar(p);

#endif

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPNormBegin"

/*@

   IPNormBegin - Starts a split phase norm computation.

   Input Parameters:

+  ip   - the inner product context

.  x    - input vector

-  norm - where the result will go

   Level: developer

   Notes:

   Each call to IPNormBegin() should be paired with a call to IPNormEnd().

.seealso: IPNormEnd(), IPNorm(), IPInnerProduct(), IPMInnerProduct(),

          IPInnerProductBegin(), IPInnerProductEnd()

@*/

PetscErrorCode IPNormBegin(IP ip,Vec x,PetscReal *norm)

{

  PetscErrorCode ierr;

  PetscScalar    p;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,2);

  PetscValidPointer(norm,3);

  if (!ip->matrix && ip->bilinear_form == IP_INNER_HERMITIAN) {

    ierr = VecNormBegin(x,NORM_2,norm);CHKERRQ(ierr);

  } else {

    ierr = IPInnerProductBegin(ip,x,x,&p);CHKERRQ(ierr);

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPNormEnd"

/*@

   IPNormEnd - Ends a split phase norm computation.

   Input Parameters:

+  ip   - the inner product context

-  x    - input vector

   Output Parameter:

.  norm - the computed norm

   Level: developer

   Notes:

   Each call to IPNormBegin() should be paired with a call to IPNormEnd().

.seealso: IPNormBegin(), IPNorm(), IPInnerProduct(), IPMInnerProduct(),

          IPInnerProductBegin(), IPInnerProductEnd()

@*/

PetscErrorCode IPNormEnd(IP ip,Vec x,PetscReal *norm)

{

  PetscErrorCode ierr;

  PetscScalar    p;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,2);

  PetscValidPointer(norm,3);

  if (!ip->matrix && ip->bilinear_form == IP_INNER_HERMITIAN) {

    ierr = VecNormEnd(x,NORM_2,norm);CHKERRQ(ierr);

  } else {

    ierr = IPInnerProductEnd(ip,x,x,&p);CHKERRQ(ierr);

    if (PetscAbsScalar(p)<PETSC_MACHINE_EPSILON)

      PetscInfo(ip,"Zero norm, either the vector is zero or a semi-inner product is being used\n");

#if defined(PETSC_USE_COMPLEX)

    if (PetscRealPart(p)<0.0 || PetscAbsReal(PetscImaginaryPart(p))>PETSC_MACHINE_EPSILON)

       SETERRQ(((PetscObject)ip)->comm,1,"IPNorm: The inner product is not well defined");

    *norm = PetscSqrtScalar(PetscRealPart(p));

#else

    if (p<0.0) SETERRQ(((PetscObject)ip)->comm,1,"IPNorm: The inner product is not well defined");

    *norm = PetscSqrtScalar(p);

#endif

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPInnerProduct"

/*@

   IPInnerProduct - Computes the inner product of two vectors.

   Collective on IP and Vec

   Input Parameters:

+  ip - the inner product context

.  x  - input vector

-  y  - input vector

   Output Parameter:

.  p - result of the inner product

   Notes:

   This function will usually compute the standard dot product of vectors

   x and y, (x,y)=y^H x. However this behaviour may be different if changed

   via IPSetBilinearForm(). This allows use of other inner products such as

   the indefinite product y^T x for complex symmetric problems or the

   B-inner product for positive definite B, (x,y)_B=y^H Bx.

   Level: developer

.seealso: IPSetBilinearForm(), VecDot(), IPMInnerProduct()

@*/

PetscErrorCode IPInnerProduct(IP ip,Vec x,Vec y,PetscScalar *p)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,2);

  PetscValidHeaderSpecific(y,VEC_CLASSID,3);

  PetscValidScalarPointer(p,4);

  ierr = PetscLogEventBegin(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  ip->innerproducts++;

  if (ip->matrix) {

    ierr = IPApplyMatrix_Private(ip,x);CHKERRQ(ierr);

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecDot(ip->Bx,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecTDot(ip->Bx,y,p);CHKERRQ(ierr);

    }

  } else {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecDot(x,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecTDot(x,y,p);CHKERRQ(ierr);

    }

  }

  ierr = PetscLogEventEnd(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPInnerProductBegin"

/*@

   IPInnerProductBegin - Starts a split phase inner product computation.

   Input Parameters:

+  ip - the inner product context

.  x  - the first vector

.  y  - the second vector

-  p  - where the result will go

   Level: developer

   Notes:

   Each call to IPInnerProductBegin() should be paired with a call to IPInnerProductEnd().

.seealso: IPInnerProductEnd(), IPInnerProduct(), IPNorm(), IPNormBegin(),

          IPNormEnd(), IPMInnerProduct()

@*/

PetscErrorCode IPInnerProductBegin(IP ip,Vec x,Vec y,PetscScalar *p)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,2);

  PetscValidHeaderSpecific(y,VEC_CLASSID,3);

  PetscValidScalarPointer(p,4);

  ierr = PetscLogEventBegin(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  ip->innerproducts++;

  if (ip->matrix) {

    ierr = IPApplyMatrix_Private(ip,x);CHKERRQ(ierr);

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecDotBegin(ip->Bx,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecTDotBegin(ip->Bx,y,p);CHKERRQ(ierr);

    }

  } else {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecDotBegin(x,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecTDotBegin(x,y,p);CHKERRQ(ierr);

    }

  }

  ierr = PetscLogEventEnd(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPInnerProductEnd"

/*@

   IPInnerProductEnd - Ends a split phase inner product computation.

   Input Parameters:

+  ip - the inner product context

.  x  - the first vector

-  y  - the second vector

   Output Parameter:

.  p  - result of the inner product

   Level: developer

   Notes:

   Each call to IPInnerProductBegin() should be paired with a call to IPInnerProductEnd().

.seealso: IPInnerProductBegin(), IPInnerProduct(), IPNorm(), IPNormBegin(),

          IPNormEnd(), IPMInnerProduct()

@*/

PetscErrorCode IPInnerProductEnd(IP ip,Vec x,Vec y,PetscScalar *p)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,2);

  PetscValidHeaderSpecific(y,VEC_CLASSID,3);

  PetscValidScalarPointer(p,4);

  ierr = PetscLogEventBegin(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  if (ip->matrix) {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecDotEnd(ip->Bx,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecTDotEnd(ip->Bx,y,p);CHKERRQ(ierr);

    }

  } else {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecDotEnd(x,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecTDotEnd(x,y,p);CHKERRQ(ierr);

    }

  }

  ierr = PetscLogEventEnd(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPMInnerProduct"

/*@

   IPMInnerProduct - Computes the inner products a vector x with a set of

   vectors (columns of Y).

   Collective on IP and Vec

   Input Parameters:

+  ip - the inner product context

.  x  - the first input vector

.  n  - number of vectors in y

-  y  - array of vectors

   Output Parameter:

.  p - result of the inner products

   Notes:

   This function will usually compute the standard dot product of x and y_i,

   (x,y_i)=y_i^H x, for each column of Y. However this behaviour may be different

   if changed via IPSetBilinearForm(). This allows use of other inner products

   such as the indefinite product y_i^T x for complex symmetric problems or the

   B-inner product for positive definite B, (x,y_i)_B=y_i^H Bx.

   Level: developer

.seealso: IPSetBilinearForm(), VecMDot(), IPInnerProduct()

@*/

PetscErrorCode IPMInnerProduct(IP ip,Vec x,PetscInt n,const Vec y[],PetscScalar *p)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,3);

  PetscValidPointer(y,4);

  PetscValidHeaderSpecific(*y,VEC_CLASSID,4);

  PetscValidScalarPointer(p,5);

  ierr = PetscLogEventBegin(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  ip->innerproducts += n;

  if (ip->matrix) {

    ierr = IPApplyMatrix_Private(ip,x);CHKERRQ(ierr);

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecMDot(ip->Bx,n,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecMTDot(ip->Bx,n,y,p);CHKERRQ(ierr);

    }

  } else {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecMDot(x,n,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecMTDot(x,n,y,p);CHKERRQ(ierr);

    }

  }

  ierr = PetscLogEventEnd(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPMInnerProductBegin"

/*@

   IPMInnerProductBegin - Starts a split phase multiple inner product computation.

   Input Parameters:

+  ip - the inner product context

.  x  - the first input vector

.  n  - number of vectors in y

.  y  - array of vectors

-  p  - where the result will go

   Level: developer

   Notes:

   Each call to IPMInnerProductBegin() should be paired with a call to IPMInnerProductEnd().

.seealso: IPMInnerProductEnd(), IPMInnerProduct(), IPNorm(), IPNormBegin(),

          IPNormEnd(), IPInnerProduct()

@*/

PetscErrorCode IPMInnerProductBegin(IP ip,Vec x,PetscInt n,const Vec y[],PetscScalar *p)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,3);

  if (n == 0) PetscFunctionReturn(0);

  PetscValidPointer(y,4);

  PetscValidHeaderSpecific(*y,VEC_CLASSID,4);

  PetscValidScalarPointer(p,5);

  ierr = PetscLogEventBegin(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  ip->innerproducts += n;

  if (ip->matrix) {

    ierr = IPApplyMatrix_Private(ip,x);CHKERRQ(ierr);

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecMDotBegin(ip->Bx,n,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecMTDotBegin(ip->Bx,n,y,p);CHKERRQ(ierr);

    }

  } else {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecMDotBegin(x,n,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecMTDotBegin(x,n,y,p);CHKERRQ(ierr);

    }

  }

  ierr = PetscLogEventEnd(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "IPMInnerProductEnd"

/*@

   IPMInnerProductEnd - Ends a split phase multiple inner product computation.

   Input Parameters:

+  ip - the inner product context

.  x  - the first input vector

.  n  - number of vectors in y

-  y  - array of vectors

   Output Parameter:

.  p - result of the inner products

   Level: developer

   Notes:

   Each call to IPMInnerProductBegin() should be paired with a call to IPMInnerProductEnd().

.seealso: IPMInnerProductBegin(), IPMInnerProduct(), IPNorm(), IPNormBegin(),

          IPNormEnd(), IPInnerProduct()

@*/

PetscErrorCode IPMInnerProductEnd(IP ip,Vec x,PetscInt n,const Vec y[],PetscScalar *p)

{

  PetscErrorCode ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(ip,IP_CLASSID,1);

  PetscValidHeaderSpecific(x,VEC_CLASSID,3);

  if (n == 0) PetscFunctionReturn(0);

  PetscValidPointer(y,4);

  PetscValidHeaderSpecific(*y,VEC_CLASSID,4);

  PetscValidScalarPointer(p,5);

  ierr = PetscLogEventBegin(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  if (ip->matrix) {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecMDotEnd(ip->Bx,n,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecMTDotEnd(ip->Bx,n,y,p);CHKERRQ(ierr);

    }

  } else {

    if (ip->bilinear_form == IP_INNER_HERMITIAN) {

      ierr = VecMDotEnd(x,n,y,p);CHKERRQ(ierr);

    } else {

      ierr = VecMTDotEnd(x,n,y,p);CHKERRQ(ierr);

    }

  }

  ierr = PetscLogEventEnd(IP_InnerProduct,ip,x,0,0);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}