Subversion Repositories slepc-dev

Rev

Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 2573 Rev 2575
/*
/*
     Routines related to orthogonalization.
     Routines related to orthogonalization.
     See the SLEPc Technical Report STR-1 for a detailed explanation.
     See the SLEPc Technical Report STR-1 for a detailed explanation.
 
 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   SLEPc - Scalable Library for Eigenvalue Problem Computations
   SLEPc - Scalable Library for Eigenvalue Problem Computations
   Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
   Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
 
 
   This file is part of SLEPc.
   This file is part of SLEPc.
     
     
   SLEPc is free software: you can redistribute it and/or modify it under  the
   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
   terms of version 3 of the GNU Lesser General Public License as published by
   the Free Software Foundation.
   the Free Software Foundation.
 
 
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
   more details.
   more details.
 
 
   You  should have received a copy of the GNU Lesser General  Public  License
   You  should have received a copy of the GNU Lesser General  Public  License
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
*/
 
 
#include "private/ipimpl.h"      /*I "slepcip.h" I*/
#include "private/ipimpl.h"      /*I "slepcip.h" I*/
#include "slepcblaslapack.h"
#include "slepcblaslapack.h"
#include "../src/eps/impls/davidson/common/davidson.h"
#include "../src/eps/impls/davidson/common/davidson.h"
 
 
/*
/*
   IPOrthogonalizeCGS1 - Compute |v'| (estimated), |v| and one step of CGS with only one global synchronization
   IPOrthogonalizeCGS1 - Compute |v'| (estimated), |v| and one step of CGS with only one global synchronization
*/
*/
#undef __FUNCT__  
#undef __FUNCT__  
#define __FUNCT__ "IPBOrthogonalizeCGS1"
#define __FUNCT__ "IPBOrthogonalizeCGS1"
PetscErrorCode IPBOrthogonalizeCGS1(IP ip,PetscInt nds,Vec *DS,Vec *BDS,PetscInt n,PetscBool *which,Vec *V,Vec *BV,Vec v,Vec Bv,PetscScalar *H,PetscReal *onorm,PetscReal *norm)
PetscErrorCode IPBOrthogonalizeCGS1(IP ip,PetscInt nds,Vec *DS,Vec *BDS,PetscInt n,PetscBool *which,Vec *V,Vec *BV,Vec v,Vec Bv,PetscScalar *H,PetscReal *onorm,PetscReal *norm)
{
{
  PetscErrorCode ierr;
  PetscErrorCode ierr;
  PetscInt       i,j;
  PetscInt       i,j;
  PetscScalar    alpha;
  PetscScalar    alpha;
 
 
  PetscFunctionBegin;
  PetscFunctionBegin;
  /* h = [DS V]^* Bv ; alpha = (v , v) */
  /* h = [DS V]^* Bv ; alpha = (v , v) */
  ierr = VecsMultIa(H,0,nds,DS,0,nds,&Bv,0,1);CHKERRQ(ierr); j = nds;
  ierr = VecsMultIa(H,0,nds,DS,0,nds,&Bv,0,1);CHKERRQ(ierr); j = nds;
  if (!which) {
  if (!which) {
    ierr = VecsMultIa(H+j,0,n,V,0,n,&Bv,0,1);CHKERRQ(ierr); j+= n;
    ierr = VecsMultIa(H+j,0,n,V,0,n,&Bv,0,1);CHKERRQ(ierr); j+= n;
  } else {
  } else {
    for (i=0; i<n; i++) {
    for (i=0; i<n; i++) {
      if (which[i]) {
      if (which[i]) {
        ierr = VecsMultIa(H+j,0,1,V+i,0,1,&Bv,0,1);CHKERRQ(ierr); j++;
        ierr = VecsMultIa(H+j,0,1,V+i,0,1,&Bv,0,1);CHKERRQ(ierr); j++;
      }
      }
    }
    }
  }
  }
  if (onorm || norm) {
  if (onorm || norm) {
    ierr = VecsMultIa(H+j,0,1,&v,0,1,&Bv,0,1);CHKERRQ(ierr); j++;
    ierr = VecsMultIa(H+j,0,1,&v,0,1,&Bv,0,1);CHKERRQ(ierr); j++;
  }
  }
  ierr = VecsMultIb(H,0,j,j,1,PETSC_NULL,v);CHKERRQ(ierr);
  ierr = VecsMultIb(H,0,j,j,1,PETSC_NULL,v);CHKERRQ(ierr);
  if (onorm || norm) alpha = H[j-1];
  if (onorm || norm) alpha = H[j-1];
 
 
  /* v = v - V h */
  /* v = v - V h */
  ierr = SlepcVecMAXPBY(v,1.0,-1.0,nds,H,DS);CHKERRQ(ierr);
  ierr = SlepcVecMAXPBY(v,1.0,-1.0,nds,H,DS);CHKERRQ(ierr);
  if (which) {
  if (which) {
    for (j=0; j<n; j++)
    for (j=0; j<n; j++)
      if (which[j]) { ierr = VecAXPBY(v,-H[nds+j],1.0,V[j]);CHKERRQ(ierr); }
      if (which[j]) { ierr = VecAXPBY(v,-H[nds+j],1.0,V[j]);CHKERRQ(ierr); }
  } else {
  } else {
    ierr = SlepcVecMAXPBY(v,1.0,-1.0,n,H+nds,V);CHKERRQ(ierr);
    ierr = SlepcVecMAXPBY(v,1.0,-1.0,n,H+nds,V);CHKERRQ(ierr);
  }
  }
 
 
  /* Bv = Bv - BV h */
  /* Bv = Bv - BV h */
  ierr = SlepcVecMAXPBY(Bv,1.0,-1.0,nds,H,BDS);CHKERRQ(ierr);
  ierr = SlepcVecMAXPBY(Bv,1.0,-1.0,nds,H,BDS);CHKERRQ(ierr);
  if (which) {
  if (which) {
    for (j=0; j<n; j++)
    for (j=0; j<n; j++)
      if (which[j]) { ierr = VecAXPBY(Bv,-H[nds+j],1.0,BV[j]);CHKERRQ(ierr); }
      if (which[j]) { ierr = VecAXPBY(Bv,-H[nds+j],1.0,BV[j]);CHKERRQ(ierr); }
  } else {
  } else {
    ierr = SlepcVecMAXPBY(Bv,1.0,-1.0,n,H+nds,BV);CHKERRQ(ierr);
    ierr = SlepcVecMAXPBY(Bv,1.0,-1.0,n,H+nds,BV);CHKERRQ(ierr);
  }
  }
   
   
  /* compute |v| */
  /* compute |v| */
  if (onorm) *onorm = PetscSqrtReal(PetscRealPart(alpha));
  if (onorm) *onorm = PetscSqrtReal(PetscRealPart(alpha));
 
 
  /* compute |v'| */
  /* compute |v'| */
  if (norm) {
  if (norm) {
    ierr = VecDot(Bv, v, &alpha); CHKERRQ(ierr);
    ierr = VecDot(Bv, v, &alpha); CHKERRQ(ierr);
    *norm = PetscSqrtReal(PetscRealPart(alpha));
    *norm = PetscSqrtReal(PetscRealPart(alpha));
  }
  }
  PetscFunctionReturn(0);
  PetscFunctionReturn(0);
}
}
 
 
/*
/*
  IPOrthogonalizeCGS - Orthogonalize with classical Gram-Schmidt
  IPOrthogonalizeCGS - Orthogonalize with classical Gram-Schmidt
*/
*/
#undef __FUNCT__  
#undef __FUNCT__  
#define __FUNCT__ "IPBOrthogonalizeCGS"
#define __FUNCT__ "IPBOrthogonalizeCGS"
static PetscErrorCode IPBOrthogonalizeCGS(IP ip,PetscInt nds,Vec *DS,Vec *BDS,PetscInt n,PetscBool *which,Vec *V,Vec *BV,Vec v,Vec Bv,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
static PetscErrorCode IPBOrthogonalizeCGS(IP ip,PetscInt nds,Vec *DS,Vec *BDS,PetscInt n,PetscBool *which,Vec *V,Vec *BV,Vec v,Vec Bv,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
{
{
  PetscErrorCode ierr;
  PetscErrorCode ierr;
  PetscScalar    lh[100],*h,lc[100],*c,alpha;
  PetscScalar    lh[100],*h,lc[100],*c,alpha;
  PetscBool      allocatedh = PETSC_FALSE,allocatedc = PETSC_FALSE;
  PetscBool      allocatedh = PETSC_FALSE,allocatedc = PETSC_FALSE;
  PetscReal      onrm,nrm;
  PetscReal      onrm,nrm;
  PetscInt       j,k;
  PetscInt       j,k;
 
 
  PetscFunctionBegin;
  PetscFunctionBegin;
  /* allocate h and c if needed */
  /* allocate h and c if needed */
  if (!H || nds>0) {
  if (!H || nds>0) {
    if (nds+n+1<=100) h = lh;
    if (nds+n+1<=100) h = lh;
    else {
    else {
      ierr = PetscMalloc((nds+n+1)*sizeof(PetscScalar),&h);CHKERRQ(ierr);
      ierr = PetscMalloc((nds+n+1)*sizeof(PetscScalar),&h);CHKERRQ(ierr);
      allocatedh = PETSC_TRUE;
      allocatedh = PETSC_TRUE;
    }
    }
  } else h = H;
  } else h = H;
  if (ip->orthog_ref != IP_ORTHOG_REFINE_NEVER) {
  if (ip->orthog_ref != IP_ORTHOG_REFINE_NEVER) {
    if (nds+n+1<=100) c = lc;
    if (nds+n+1<=100) c = lc;
    else {
    else {
      ierr = PetscMalloc((nds+n+1)*sizeof(PetscScalar),&c);CHKERRQ(ierr);
      ierr = PetscMalloc((nds+n+1)*sizeof(PetscScalar),&c);CHKERRQ(ierr);
      allocatedc = PETSC_TRUE;
      allocatedc = PETSC_TRUE;
    }
    }
  }
  }
 
 
  /* orthogonalize and compute onorm */
  /* orthogonalize and compute onorm */
  switch (ip->orthog_ref) {
  switch (ip->orthog_ref) {
 
 
  case IP_ORTHOG_REFINE_NEVER:
  case IP_ORTHOG_REFINE_NEVER:
    ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,h,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
    ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,h,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
    /* compute |v| */
    /* compute |v| */
    if (norm) {
    if (norm) {
      ierr = VecDot(Bv,v,&alpha); CHKERRQ(ierr);
      ierr = VecDot(Bv,v,&alpha); CHKERRQ(ierr);
      *norm = PetscSqrtReal(PetscRealPart(alpha));
      *norm = PetscSqrtReal(PetscRealPart(alpha));
    }
    }
    /* linear dependence check does not work without refinement */
    /* linear dependence check does not work without refinement */
    if (lindep) *lindep = PETSC_FALSE;
    if (lindep) *lindep = PETSC_FALSE;
    break;
    break;
   
   
  case IP_ORTHOG_REFINE_ALWAYS:
  case IP_ORTHOG_REFINE_ALWAYS:
    ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,h,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
    ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,h,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
    if (lindep) {
    if (lindep) {
      ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,&onrm,&nrm);CHKERRQ(ierr);
      ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,&onrm,&nrm);CHKERRQ(ierr);
      if (norm) *norm = nrm;
      if (norm) *norm = nrm;
      if (nrm < ip->orthog_eta * onrm) *lindep = PETSC_TRUE;
      if (nrm < ip->orthog_eta * onrm) *lindep = PETSC_TRUE;
      else *lindep = PETSC_FALSE;
      else *lindep = PETSC_FALSE;
    } else {
    } else {
      ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,PETSC_NULL,norm);CHKERRQ(ierr);
      ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,PETSC_NULL,norm);CHKERRQ(ierr);
    }
    }
    for (j=0;j<n;j++)
    for (j=0;j<n;j++)
      if (!which || which[j]) h[nds+j] += c[nds+j];
      if (!which || which[j]) h[nds+j] += c[nds+j];
    break;
    break;
 
 
  case IP_ORTHOG_REFINE_IFNEEDED:
  case IP_ORTHOG_REFINE_IFNEEDED:
    ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,h,&onrm,&nrm);CHKERRQ(ierr);
    ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,h,&onrm,&nrm);CHKERRQ(ierr);
    /* ||q|| < eta ||h|| */
    /* ||q|| < eta ||h|| */
    k = 1;
    k = 1;
    while (k<3 && nrm < ip->orthog_eta * onrm) {
    while (k<3 && nrm < ip->orthog_eta * onrm) {
      k++;
      k++;
      if (!ip->matrix) {
      if (!ip->matrix) {
        ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,&onrm,&nrm);CHKERRQ(ierr);
        ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,&onrm,&nrm);CHKERRQ(ierr);
      } else {
      } else {
        onrm = nrm;
        onrm = nrm;
        ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,PETSC_NULL,&nrm);CHKERRQ(ierr);
        ierr = IPBOrthogonalizeCGS1(ip,nds,DS,BDS,n,which,V,BV,v,Bv,c,PETSC_NULL,&nrm);CHKERRQ(ierr);
      }
      }
      for (j=0;j<n;j++)
      for (j=0;j<n;j++)
        if (!which || which[j]) h[nds+j] += c[nds+j];
        if (!which || which[j]) h[nds+j] += c[nds+j];
    }
    }
    if (norm) *norm = nrm;
    if (norm) *norm = nrm;
    if (lindep) {
    if (lindep) {
      if (nrm < ip->orthog_eta * onrm) *lindep = PETSC_TRUE;
      if (nrm < ip->orthog_eta * onrm) *lindep = PETSC_TRUE;
      else *lindep = PETSC_FALSE;
      else *lindep = PETSC_FALSE;
    }
    }
    break;
    break;
 
 
  default:
  default:
    SETERRQ(((PetscObject)ip)->comm,PETSC_ERR_ARG_WRONG,"Unknown orthogonalization refinement");
    SETERRQ(((PetscObject)ip)->comm,PETSC_ERR_ARG_WRONG,"Unknown orthogonalization refinement");
  }
  }
 
 
  /* recover H from workspace */
  /* recover H from workspace */
  if (H && nds>0) {
  if (H && nds>0) {
    for (j=0;j<n;j++)
    for (j=0;j<n;j++)
      if (!which || which[j]) H[j] = h[nds+j];
      if (!which || which[j]) H[j] = h[nds+j];
  }
  }
 
 
  /* free work space */
  /* free work space */
  if (allocatedc) { ierr = PetscFree(c);CHKERRQ(ierr); }
  if (allocatedc) { ierr = PetscFree(c);CHKERRQ(ierr); }
  if (allocatedh) { ierr = PetscFree(h);CHKERRQ(ierr); }
  if (allocatedh) { ierr = PetscFree(h);CHKERRQ(ierr); }
  PetscFunctionReturn(0);
  PetscFunctionReturn(0);
}        
}        
 
 
#undef __FUNCT__  
#undef __FUNCT__  
#define __FUNCT__ "IPBOrthogonalize"
#define __FUNCT__ "IPBOrthogonalize"
/*@
/*@
   IPBOrthogonalize - B-Orthogonalize a vector with respect to two set of vectors.
   IPBOrthogonalize - B-Orthogonalize a vector with respect to two set of vectors.
 
 
   Collective on IP
   Collective on IP
 
 
   Input Parameters:
   Input Parameters:
+  ip     - the inner product (IP) context
+  ip     - the inner product (IP) context
.  nds    - number of columns of DS
.  nds    - number of columns of DS
.  DS     - first set of vectors
.  DS     - first set of vectors
.  BDS    - B * DS
.  BDS    - B * DS
.  n      - number of columns of V
.  n      - number of columns of V
.  which  - logical array indicating columns of V to be used
.  which  - logical array indicating columns of V to be used
.  V      - second set of vectors
.  V      - second set of vectors
-  BV     - B * V
-  BV     - B * V
 
 
   Input/Output Parameter:
   Input/Output Parameter:
+  v      - (input) vector to be orthogonalized and (output) result of
+  v      - (input) vector to be orthogonalized and (output) result of
            orthogonalization
            orthogonalization
-  Bv     - (input/output) B * v
-  Bv     - (input/output) B * v
 
 
   Output Parameter:
   Output Parameter:
+  H      - coefficients computed during orthogonalization with V, of size nds+n
+  H      - coefficients computed during orthogonalization with V, of size nds+n
            if norm == PETSC_NULL, and nds+n+1 otherwise.
            if norm == PETSC_NULL, and nds+n+1 otherwise.
.  norm   - norm of the vector after being orthogonalized
.  norm   - norm of the vector after being orthogonalized
-  lindep - flag indicating that refinement did not improve the quality
-  lindep - flag indicating that refinement did not improve the quality
            of orthogonalization
            of orthogonalization
 
 
   Notes:
   Notes:
   This function applies an orthogonal projector to project vector v onto the
   This function applies an orthogonal projector to project vector v onto the
   orthogonal complement of the span of the columns of DS and V.
   orthogonal complement of the span of the columns of DS and V.
 
 
   On exit, v0 = [V v]*H, where v0 is the original vector v.
   On exit, v0 = [V v]*H, where v0 is the original vector v.
 
 
   This routine does not normalize the resulting vector.
   This routine does not normalize the resulting vector.
 
 
   Level: developer
   Level: developer
 
 
.seealso: IPSetOrthogonalization(), IPBiOrthogonalize()
.seealso: IPSetOrthogonalization(), IPBiOrthogonalize()
@*/
@*/
PetscErrorCode IPBOrthogonalize(IP ip,PetscInt nds,Vec *DS, Vec *BDS,PetscInt n,PetscBool *which,Vec *V,Vec *BV,Vec v,Vec Bv,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
PetscErrorCode IPBOrthogonalize(IP ip,PetscInt nds,Vec *DS, Vec *BDS,PetscInt n,PetscBool *which,Vec *V,Vec *BV,Vec v,Vec Bv,PetscScalar *H,PetscReal *norm,PetscBool *lindep)
{
{
  PetscErrorCode ierr;
  PetscErrorCode ierr;
  PetscScalar    alpha;
  PetscScalar    alpha;
 
 
  PetscFunctionBegin;
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ip,IP_CLASSID,1);
  PetscValidHeaderSpecific(ip,IP_CLASSID,1);
  PetscValidLogicalCollectiveInt(ip,nds,2);
  PetscValidLogicalCollectiveInt(ip,nds,2);
  PetscValidLogicalCollectiveInt(ip,n,4);
  PetscValidLogicalCollectiveInt(ip,n,4);
  ierr = PetscLogEventBegin(IP_Orthogonalize,ip,0,0,0);CHKERRQ(ierr);
  ierr = PetscLogEventBegin(IP_Orthogonalize,ip,0,0,0);CHKERRQ(ierr);
 
 
  /* Bv <- B * v */
  /* Bv <- B * v */
  ierr = PetscLogEventBegin(IP_ApplyMatrix,ip,0,0,0);CHKERRQ(ierr);
  ierr = PetscLogEventBegin(IP_ApplyMatrix,ip,0,0,0);CHKERRQ(ierr);
  ierr = MatMult(ip->matrix, v, Bv); CHKERRQ(ierr);
  ierr = MatMult(ip->matrix, v, Bv); CHKERRQ(ierr);
  ierr = PetscLogEventEnd(IP_ApplyMatrix,ip,0,0,0);CHKERRQ(ierr);
  ierr = PetscLogEventEnd(IP_ApplyMatrix,ip,0,0,0);CHKERRQ(ierr);
   
   
  if (nds==0 && n==0) {
  if (nds==0 && n==0) {
    if (norm) {
    if (norm) {
      ierr = VecDot(Bv, v, &alpha); CHKERRQ(ierr);
      ierr = VecDot(Bv, v, &alpha); CHKERRQ(ierr);
      *norm = PetscSqrtReal(PetscRealPart(alpha));
      *norm = PetscSqrtReal(PetscRealPart(alpha));
    }
    }
    if (lindep) *lindep = PETSC_FALSE;
    if (lindep) *lindep = PETSC_FALSE;
  } else {
  } else {
    switch (ip->orthog_type) {
    switch (ip->orthog_type) {
    case IP_ORTHOG_CGS:
    case IP_ORTHOG_CGS:
      ierr = IPBOrthogonalizeCGS(ip,nds,DS,BDS,n,which,V,BV,v,Bv,H,norm,lindep);CHKERRQ(ierr);
      ierr = IPBOrthogonalizeCGS(ip,nds,DS,BDS,n,which,V,BV,v,Bv,H,norm,lindep);CHKERRQ(ierr);
      break;
      break;
    case IP_ORTHOG_MGS:
    case IP_ORTHOG_MGS:
      SETERRQ(((PetscObject)ip)->comm,PETSC_ERR_ARG_WRONG,"Unsupported orthogonalization type");
      SETERRQ(((PetscObject)ip)->comm,PETSC_ERR_ARG_WRONG,"Unsupported orthogonalization type");
      break;
      break;
    default:
    default:
      SETERRQ(((PetscObject)ip)->comm,PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
      SETERRQ(((PetscObject)ip)->comm,PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
    }
    }
  }
  }
  ierr = PetscLogEventEnd(IP_Orthogonalize,ip,0,0,0);CHKERRQ(ierr);  
  ierr = PetscLogEventEnd(IP_Orthogonalize,ip,0,0,0);CHKERRQ(ierr);  
  PetscFunctionReturn(0);
  PetscFunctionReturn(0);
}
}