Subversion Repositories slepc-dev

Rev

Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 2599 Rev 2609
Line 378... Line 378...
      }
      }
      if(s){
      if(s){
        d_prev = PetscAbsReal( (sPres->value - s->value)/(sPres->inertia - s->inertia));
        d_prev = PetscAbsReal( (sPres->value - s->value)/(sPres->inertia - s->inertia));
      }else{/* First shift. Average distance obtained with values in this shift */
      }else{/* First shift. Average distance obtained with values in this shift */
        /* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
        /* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
        if( (sr->dir)*(sr->eig[0]-sPres->value)>0 && PetscAbsReal( (PetscRealPart(sr->eig[sr->indexEig-1]) - PetscRealPart(sr->eig[0]))/PetscRealPart(sr->eig[0])) > PetscSqrtReal(eps->tol) ){
        if( (sr->dir)*(PetscRealPart(sr->eig[0])-sPres->value)>0 && PetscAbsReal( (PetscRealPart(sr->eig[sr->indexEig-1]) - PetscRealPart(sr->eig[0]))/PetscRealPart(sr->eig[0])) > PetscSqrtReal(eps->tol) ){
          d_prev =  PetscAbsReal( (PetscRealPart(sr->eig[sr->indexEig-1]) - PetscRealPart(sr->eig[0])))/(sPres->neigs+0.3);
          d_prev =  PetscAbsReal( (PetscRealPart(sr->eig[sr->indexEig-1]) - PetscRealPart(sr->eig[0])))/(sPres->neigs+0.3);
        }else{
        }else{
          d_prev = PetscAbsReal( PetscRealPart(sr->eig[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);          
          d_prev = PetscAbsReal( PetscRealPart(sr->eig[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);          
        }
        }
      }
      }
Line 393... Line 393...
        *newS = sPres->value + ((sr->dir)*d_prev*(eps->nev))/2;
        *newS = sPres->value + ((sr->dir)*d_prev*(eps->nev))/2;
      }
      }
    }
    }
    /* End of interval can not be surpassed */
    /* End of interval can not be surpassed */
    if((sr->dir)*( sr->int1 - *newS) < 0) *newS = sr->int1;
    if((sr->dir)*( sr->int1 - *newS) < 0) *newS = sr->int1;
  }//of neighb[side]==null
  }/* of neighb[side]==null */
  PetscFunctionReturn(0);
  PetscFunctionReturn(0);
}
}
 
 
/*
/*
  Function for sorting an array of real values
  Function for sorting an array of real values
Line 449... Line 449...
  for( i=0; i < eps->nconv ;i++ ){
  for( i=0; i < eps->nconv ;i++ ){
    lambda = PetscRealPart(eps->eigr[eps->perm[i]]);
    lambda = PetscRealPart(eps->eigr[eps->perm[i]]);
    err = eps->errest[eps->perm[i]];
    err = eps->errest[eps->perm[i]];
    if(  (sr->dir)*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0  ){/* Valid value */
    if(  (sr->dir)*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0  ){/* Valid value */
      if(count>=sr->numEigs){/* Error found */
      if(count>=sr->numEigs){/* Error found */
         SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Unexpected error in Spectrum Slicing");
         SETERRQ(((PetscObject)eps)->comm,1,"Unexpected error in Spectrum Slicing");
      }  
      }  
      sr->eig[count] = lambda;
      sr->eig[count] = lambda;
      sr->errest[count] = err;
      sr->errest[count] = err;
      /* Purification */
      /* Purification */
      ierr = STApply(eps->OP,eps->V[eps->perm[i]],sr->V[count]);CHKERRQ(ierr);
      if (eps->isgeneralized){
      ierr = IPNorm(eps->ip,sr->V[count],&norm);CHKERRQ(ierr);
        ierr = STApply(eps->OP,eps->V[eps->perm[i]],sr->V[count]);CHKERRQ(ierr);
      ierr = VecScale(sr->V[count],1.0/norm);CHKERRQ(ierr);
        ierr = IPNorm(eps->ip,sr->V[count],&norm);CHKERRQ(ierr);
 
        ierr = VecScale(sr->V[count],1.0/norm);CHKERRQ(ierr);
 
      }else{
 
        ierr = VecCopy(eps->V[eps->perm[i]], sr->V[count]);CHKERRQ(ierr);
 
      }
      count++;
      count++;
    }
    }
  }
  }
  sPres->neigs = count - sr->indexEig;
  sPres->neigs = count - sr->indexEig;
  sr->indexEig = count;
  sr->indexEig = count;
Line 588... Line 592...
  /* Memory reservation for eig, V and perm */
  /* Memory reservation for eig, V and perm */
  ierr = PetscMalloc((sr->numEigs)*sizeof(PetscScalar),&sr->eig);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs)*sizeof(PetscScalar),&sr->eig);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs)*sizeof(PetscScalar),&sr->eigi);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs)*sizeof(PetscScalar),&sr->eigi);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs+eps->ncv) *sizeof(PetscReal),&sr->errest);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs+eps->ncv) *sizeof(PetscReal),&sr->errest);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs+eps->ncv)*sizeof(PetscReal),&errest_left);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs+eps->ncv)*sizeof(PetscReal),&errest_left);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs+eps->ncv)*sizeof(PetscReal),&sr->monit);CHKERRQ(ierr);
  ierr = PetscMalloc((sr->numEigs+eps->ncv)*sizeof(PetscScalar),&sr->monit);CHKERRQ(ierr);
  for(i=0;i<sr->numEigs;i++){sr->eigi[i]=0;sr->eig[i] = 0;}
  for(i=0;i<sr->numEigs;i++){sr->eigi[i]=0;sr->eig[i] = 0;}
  for(i=0;i<sr->numEigs+eps->ncv;i++){errest_left[i]=0;sr->errest[i]=0;sr->monit[i]=0;}
  for(i=0;i<sr->numEigs+eps->ncv;i++){errest_left[i]=0;sr->errest[i]=0;sr->monit[i]=0;}
  ierr = VecCreateMPI(((PetscObject)eps)->comm,eps->nloc,PETSC_DECIDE,&t);CHKERRQ(ierr);
  ierr = VecCreateMPI(((PetscObject)eps)->comm,eps->nloc,PETSC_DECIDE,&t);CHKERRQ(ierr);
  ierr = VecDuplicateVecs(t,sr->numEigs,&sr->V);CHKERRQ(ierr);
  ierr = VecDuplicateVecs(t,sr->numEigs,&sr->V);CHKERRQ(ierr);
  ierr = VecDestroy(&t);CHKERRQ(ierr);
  ierr = VecDestroy(&t);CHKERRQ(ierr);