renormvec.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar 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 General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: renormvec.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "util.h"

PetscReal RenormVec(Vec v, PetscReal normset, PetscReal normtol, int dim)
{
  MagparFunctionInfoBegin;

  int t_vlen;
  ierr = VecGetLocalSize(v,(PetscInt*)&t_vlen);CHKERRQ(ierr);
  assert(t_vlen%dim==0);

  PetscReal *ta_v;
  ierr = VecGetArray(v,&ta_v);CHKERRQ(ierr);

  assert(dim==ND); /* otherwise replace my_* calls by cblas calls again */

  PetscReal normmax,normmin;
  normmax=0;
  normmin=PETSC_MAX;
  for (int i=0;i<t_vlen/dim;i++) {
    PetscReal nrm;
    nrm=my_dnrm2(dim,ta_v+dim*i,1);
    normmax=PetscMax(normmax,nrm);
    normmin=PetscMin(normmin,nrm);
  }

  PetscReal gnormmax,gnormmin;
  ierr = PetscGlobalMax(&normmax,&gnormmax,PETSC_COMM_WORLD);CHKERRQ(ierr);
  ierr = PetscGlobalMin(&normmin,&gnormmin,PETSC_COMM_WORLD);CHKERRQ(ierr);

  PetscReal normdev;
  normdev = PetscMax(PetscAbsReal(gnormmax-normset),
                     PetscAbsReal(gnormmin-normset));

  /* rescale vectors only if maximum deviation is larger than tolerance */
  if (PetscAbsReal(normdev) > normtol) {
    for (int i=0;i<t_vlen/dim;i++) {
      /* recalculate norm since we did not store it anywhere */
      PetscReal nrm;
      nrm=my_dnrm2(dim,ta_v+dim*i,1);

      /* work around problems if nrm is small */
      if (nrm>D_EPS) {
        my_dscal(dim,normset/nrm,ta_v+dim*i,1);
      }
      /* if nrm is too small and normtol==0
         set vector (at will) parallel to z-axis
         TODO: good idea? problem for GMR calculation
      */
      /*
      else if (normtol==0){
      */
      else {
        ta_v[dim*i+0]=0.0;
        ta_v[dim*i+1]=0.0;
        ta_v[dim*i+2]=1.0;
      }
    }
  }

  VecRestoreArray(v,&ta_v);

  MagparFunctionInfoReturn(normdev);
}