Subversion Repositories slepc-dev

Rev

Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
6 dsic.upv.es!jroman 1 
/*
2 
    Shift spectral transformation, applies (A + sigma I) as operator, or
3 
    inv(B)(A + sigma B) for generalized problems
1376 slepc 4 
 
5 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1672 slepc 6 
   SLEPc - Scalable Library for Eigenvalue Problem Computations
2116 eromero 7 
   Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
1376 slepc 8 
 
1672 slepc 9 
   This file is part of SLEPc.
10 
 
11 
   SLEPc is free software: you can redistribute it and/or modify it under  the
12 
   terms of version 3 of the GNU Lesser General Public License as published by
13 
   the Free Software Foundation.
14 
 
15 
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
16 
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
17 
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
18 
   more details.
19 
 
20 
   You  should have received a copy of the GNU Lesser General  Public  License
21 
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
1376 slepc 22 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
6 dsic.upv.es!jroman 23 
*/
1376 slepc 24 
 
1521 slepc 25 
#include "private/stimpl.h"          /*I "slepcst.h" I*/
6 dsic.upv.es!jroman 26 
 
27 
#undef __FUNCT__  
28 
#define __FUNCT__ "STApply_Shift"
476 dsic.upv.es!antodo 29 
PetscErrorCode STApply_Shift(ST st,Vec x,Vec y)
6 dsic.upv.es!jroman 30 
{
476 dsic.upv.es!antodo 31 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 32 
 
33 
  PetscFunctionBegin;
34 
  if (st->B) {
35 
    /* generalized eigenproblem: y = (B^-1 A + sI) x */
344 dsic.upv.es!antodo 36 
    ierr = MatMult(st->A,x,st->w);CHKERRQ(ierr);
37 
    ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
6 dsic.upv.es!jroman 38 
  }
39 
  else {
40 
    /* standard eigenproblem: y = (A + sI) x */
41 
    ierr = MatMult(st->A,x,y);CHKERRQ(ierr);
109 dsic.upv.es!antodo 42 
  }
133 dsic.upv.es!antodo 43 
  if (st->sigma != 0.0) {
828 dsic.upv.es!antodo 44 
    ierr = VecAXPY(y,st->sigma,x);CHKERRQ(ierr);
6 dsic.upv.es!jroman 45 
  }
46 
  PetscFunctionReturn(0);
47 
}
48 
 
49 
#undef __FUNCT__  
780 dsic.upv.es!jroman 50 
#define __FUNCT__ "STApplyTranspose_Shift"
51 
PetscErrorCode STApplyTranspose_Shift(ST st,Vec x,Vec y)
52 
{
53 
  PetscErrorCode ierr;
54 
 
55 
  PetscFunctionBegin;
56 
  if (st->B) {
57 
    /* generalized eigenproblem: y = (A^T B^-T + sI) x */
58 
    ierr = STAssociatedKSPSolveTranspose(st,x,st->w);CHKERRQ(ierr);
59 
    ierr = MatMultTranspose(st->A,st->w,y);CHKERRQ(ierr);
60 
  }
61 
  else {
62 
    /* standard eigenproblem: y = (A^T + sI) x */
63 
    ierr = MatMultTranspose(st->A,x,y);CHKERRQ(ierr);
64 
  }
65 
  if (st->sigma != 0.0) {
828 dsic.upv.es!antodo 66 
    ierr = VecAXPY(y,st->sigma,x);CHKERRQ(ierr);
780 dsic.upv.es!jroman 67 
  }
68 
  PetscFunctionReturn(0);
69 
}
70 
 
71 
#undef __FUNCT__  
6 dsic.upv.es!jroman 72 
#define __FUNCT__ "STBackTransform_Shift"
1779 antodo 73 
PetscErrorCode STBackTransform_Shift(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
6 dsic.upv.es!jroman 74 
{
1778 antodo 75 
  PetscInt j;
6 dsic.upv.es!jroman 76 
  PetscFunctionBegin;
1778 antodo 77 
  PetscValidPointer(eigr,3);
78 
  for (j=0;j<n;j++) {
79 
    eigr[j] -= st->sigma;
80 
  }
6 dsic.upv.es!jroman 81 
  PetscFunctionReturn(0);
82 
}
83 
 
84 
#undef __FUNCT__  
85 
#define __FUNCT__ "STSetUp_Shift"
476 dsic.upv.es!antodo 86 
PetscErrorCode STSetUp_Shift(ST st)
6 dsic.upv.es!jroman 87 
{
476 dsic.upv.es!antodo 88 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 89 
 
90 
  PetscFunctionBegin;
246 dsic.upv.es!antodo 91 
  if (st->B) {
92 
    ierr = KSPSetOperators(st->ksp,st->B,st->B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
18 dsic.upv.es!jroman 93 
    ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
6 dsic.upv.es!jroman 94 
  }
95 
  PetscFunctionReturn(0);
96 
}
97 
 
438 dsic.upv.es!antodo 98 
#undef __FUNCT__  
99 
#define __FUNCT__ "STView_Shift"
476 dsic.upv.es!antodo 100 
PetscErrorCode STView_Shift(ST st,PetscViewer viewer)
438 dsic.upv.es!antodo 101 
{
476 dsic.upv.es!antodo 102 
  PetscErrorCode ierr;
438 dsic.upv.es!antodo 103 
 
104 
  PetscFunctionBegin;
105 
  if (st->B) {
106 
    ierr = STView_Default(st,viewer);CHKERRQ(ierr);
107 
  }
108 
  PetscFunctionReturn(0);
109 
}
110 
 
2005 eromero 111 
#undef __FUNCT__  
112 
#define __FUNCT__ "STSetFromOptions_Shift"
113 
PetscErrorCode STSetFromOptions_Shift(ST st)
114 
{
115 
  PetscErrorCode ierr;
116 
  PC             pc;
117 
  const PCType   pctype;
118 
  const KSPType  ksptype;
119 
 
120 
  PetscFunctionBegin;
121 
 
122 
  ierr = KSPGetPC(st->ksp,&pc);CHKERRQ(ierr);
123 
  ierr = KSPGetType(st->ksp,&ksptype);CHKERRQ(ierr);
124 
  ierr = PCGetType(pc,&pctype);CHKERRQ(ierr);
125 
  if (!pctype && !ksptype) {
126 
    if (st->shift_matrix == ST_MATMODE_SHELL) {
127 
      /* in shell mode use GMRES with Jacobi as the default */
128 
      ierr = KSPSetType(st->ksp,KSPGMRES);CHKERRQ(ierr);
129 
      ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
130 
    } else {
131 
      /* use direct solver as default */
132 
      ierr = KSPSetType(st->ksp,KSPPREONLY);CHKERRQ(ierr);
133 
      ierr = PCSetType(pc,PCREDUNDANT);CHKERRQ(ierr);
134 
    }
135 
  }
136 
 
137 
  PetscFunctionReturn(0);
138 
}
139 
 
6 dsic.upv.es!jroman 140 
EXTERN_C_BEGIN
141 
#undef __FUNCT__  
142 
#define __FUNCT__ "STCreate_Shift"
476 dsic.upv.es!antodo 143 
PetscErrorCode STCreate_Shift(ST st)
6 dsic.upv.es!jroman 144 
{
145 
  PetscFunctionBegin;
1358 slepc 146 
  st->ops->apply           = STApply_Shift;
147 
  st->ops->getbilinearform = STGetBilinearForm_Default;
148 
  st->ops->applytrans      = STApplyTranspose_Shift;
149 
  st->ops->backtr          = STBackTransform_Shift;
2005 eromero 150 
  st->ops->setfromoptions  = STSetFromOptions_Shift;
1358 slepc 151 
  st->ops->setup           = STSetUp_Shift;
152 
  st->ops->view            = STView_Shift;
153 
  st->checknullspace       = 0;
6 dsic.upv.es!jroman 154 
  PetscFunctionReturn(0);
155 
}
156 
EXTERN_C_END
157