Subversion Repositories slepc-dev

Rev

Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
6 dsic.upv.es!jroman 1 
/*
2 
      Implements the shift-and-invert technique for eigenvalue problems.
1376 slepc 3 
 
4 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1672 slepc 5 
   SLEPc - Scalable Library for Eigenvalue Problem Computations
2575 eromero 6 
   Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
1376 slepc 7 
 
1672 slepc 8 
   This file is part of SLEPc.
9 
 
10 
   SLEPc is free software: you can redistribute it and/or modify it under  the
11 
   terms of version 3 of the GNU Lesser General Public License as published by
12 
   the Free Software Foundation.
13 
 
14 
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
15 
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
16 
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
17 
   more details.
18 
 
19 
   You  should have received a copy of the GNU Lesser General  Public  License
20 
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
1376 slepc 21 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
6 dsic.upv.es!jroman 22 
*/
1376 slepc 23 
 
2283 jroman 24 
#include <private/stimpl.h>          /*I "slepcst.h" I*/
6 dsic.upv.es!jroman 25 
 
26 
#undef __FUNCT__  
27 
#define __FUNCT__ "STApply_Sinvert"
476 dsic.upv.es!antodo 28 
PetscErrorCode STApply_Sinvert(ST st,Vec x,Vec y)
6 dsic.upv.es!jroman 29 
{
476 dsic.upv.es!antodo 30 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 31 
 
32 
  PetscFunctionBegin;
33 
  if (st->B) {
34 
    /* generalized eigenproblem: y = (A - sB)^-1 B x */
344 dsic.upv.es!antodo 35 
    ierr = MatMult(st->B,x,st->w);CHKERRQ(ierr);
36 
    ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
6 dsic.upv.es!jroman 37 
  }
38 
  else {
39 
    /* standard eigenproblem: y = (A - sI)^-1 x */
18 dsic.upv.es!jroman 40 
    ierr = STAssociatedKSPSolve(st,x,y);CHKERRQ(ierr);
6 dsic.upv.es!jroman 41 
  }
42 
  PetscFunctionReturn(0);
43 
}
44 
 
45 
#undef __FUNCT__  
780 dsic.upv.es!jroman 46 
#define __FUNCT__ "STApplyTranspose_Sinvert"
47 
PetscErrorCode STApplyTranspose_Sinvert(ST st,Vec x,Vec y)
48 
{
49 
  PetscErrorCode ierr;
50 
 
51 
  PetscFunctionBegin;
52 
  if (st->B) {
53 
    /* generalized eigenproblem: y = B^T (A - sB)^-T x */
54 
    ierr = STAssociatedKSPSolveTranspose(st,x,st->w);CHKERRQ(ierr);
55 
    ierr = MatMultTranspose(st->B,st->w,y);CHKERRQ(ierr);
56 
  }
57 
  else {
58 
    /* standard eigenproblem: y = (A - sI)^-T x */
59 
    ierr = STAssociatedKSPSolveTranspose(st,x,y);CHKERRQ(ierr);
60 
  }
61 
  PetscFunctionReturn(0);
62 
}
63 
 
64 
#undef __FUNCT__  
6 dsic.upv.es!jroman 65 
#define __FUNCT__ "STBackTransform_Sinvert"
1779 antodo 66 
PetscErrorCode STBackTransform_Sinvert(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
6 dsic.upv.es!jroman 67 
{
2317 jroman 68 
  PetscInt    j;
2320 jroman 69 
#if !defined(PETSC_USE_COMPLEX)
179 dsic.upv.es!antodo 70 
  PetscScalar t;
2317 jroman 71 
#endif
72 
 
6 dsic.upv.es!jroman 73 
  PetscFunctionBegin;
2320 jroman 74 
#if !defined(PETSC_USE_COMPLEX)
1778 antodo 75 
  for (j=0;j<n;j++) {
76 
    if (eigi[j] == 0) eigr[j] = 1.0 / eigr[j] + st->sigma;
77 
    else {
78 
      t = eigr[j] * eigr[j] + eigi[j] * eigi[j];
79 
      eigr[j] = eigr[j] / t + st->sigma;
80 
      eigi[j] = - eigi[j] / t;
81 
    }
179 dsic.upv.es!antodo 82 
  }
83 
#else
1778 antodo 84 
  for (j=0;j<n;j++) {
85 
    eigr[j] = 1.0 / eigr[j] + st->sigma;
86 
  }
179 dsic.upv.es!antodo 87 
#endif
6 dsic.upv.es!jroman 88 
  PetscFunctionReturn(0);
89 
}
90 
 
91 
#undef __FUNCT__  
1029 slepc 92 
#define __FUNCT__ "STPostSolve_Sinvert"
93 
PetscErrorCode STPostSolve_Sinvert(ST st)
6 dsic.upv.es!jroman 94 
{
476 dsic.upv.es!antodo 95 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 96 
 
97 
  PetscFunctionBegin;
1940 jroman 98 
  if (st->shift_matrix == ST_MATMODE_INPLACE) {
2331 jroman 99 
    if (st->B) {
828 dsic.upv.es!antodo 100 
      ierr = MatAXPY(st->A,st->sigma,st->B,st->str);CHKERRQ(ierr);
101 
    } else {
2330 jroman 102 
      ierr = MatShift(st->A,st->sigma);CHKERRQ(ierr);
828 dsic.upv.es!antodo 103 
    }
6 dsic.upv.es!jroman 104 
    st->setupcalled = 0;
105 
  }
106 
  PetscFunctionReturn(0);
107 
}
108 
 
109 
#undef __FUNCT__  
110 
#define __FUNCT__ "STSetUp_Sinvert"
476 dsic.upv.es!antodo 111 
PetscErrorCode STSetUp_Sinvert(ST st)
6 dsic.upv.es!jroman 112 
{
476 dsic.upv.es!antodo 113 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 114 
 
115 
  PetscFunctionBegin;
2305 jroman 116 
  ierr = MatDestroy(&st->mat);CHKERRQ(ierr);
6 dsic.upv.es!jroman 117 
 
2074 jroman 118 
  /* if the user did not set the shift, use the target value */
119 
  if (!st->sigma_set) st->sigma = st->defsigma;
120 
 
2364 jroman 121 
  if (!st->ksp) { ierr = STGetKSP(st,&st->ksp);CHKERRQ(ierr); }
344 dsic.upv.es!antodo 122 
  switch (st->shift_matrix) {
1940 jroman 123 
  case ST_MATMODE_INPLACE:
360 dsic.upv.es!antodo 124 
    st->mat = PETSC_NULL;
331 dsic.upv.es!antodo 125 
    if (st->sigma != 0.0) {
326 dsic.upv.es!antodo 126 
      if (st->B) {
828 dsic.upv.es!antodo 127 
        ierr = MatAXPY(st->A,-st->sigma,st->B,st->str);CHKERRQ(ierr);
326 dsic.upv.es!antodo 128 
      } else {
828 dsic.upv.es!antodo 129 
        ierr = MatShift(st->A,-st->sigma);CHKERRQ(ierr);
326 dsic.upv.es!antodo 130 
      }
131 
    }
1241 slepc 132 
    ierr = KSPSetOperators(st->ksp,st->A,st->A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
163 dsic.upv.es!antodo 133 
    break;
1940 jroman 134 
  case ST_MATMODE_SHELL:
344 dsic.upv.es!antodo 135 
    ierr = STMatShellCreate(st,&st->mat);CHKERRQ(ierr);
136 
    ierr = KSPSetOperators(st->ksp,st->mat,st->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
163 dsic.upv.es!antodo 137 
    break;
138 
  default:
331 dsic.upv.es!antodo 139 
    if (st->sigma != 0.0) {
1418 slepc 140 
      ierr = MatDuplicate(st->A,MAT_COPY_VALUES,&st->mat);CHKERRQ(ierr);
2654 jroman 141 
      if (st->B) { ierr = MatAXPY(st->mat,-st->sigma,st->B,st->str);CHKERRQ(ierr); }
142 
      else { ierr = MatShift(st->mat,-st->sigma);CHKERRQ(ierr); }
1241 slepc 143 
      ierr = KSPSetOperators(st->ksp,st->mat,st->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
144 
    } else {
145 
      st->mat = PETSC_NULL;
146 
      ierr = KSPSetOperators(st->ksp,st->A,st->A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
326 dsic.upv.es!antodo 147 
    }
6 dsic.upv.es!jroman 148 
  }
18 dsic.upv.es!jroman 149 
  ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
6 dsic.upv.es!jroman 150 
  PetscFunctionReturn(0);
151 
}
152 
 
153 
#undef __FUNCT__  
154 
#define __FUNCT__ "STSetShift_Sinvert"
476 dsic.upv.es!antodo 155 
PetscErrorCode STSetShift_Sinvert(ST st,PetscScalar newshift)
6 dsic.upv.es!jroman 156 
{
476 dsic.upv.es!antodo 157 
  PetscErrorCode ierr;
1241 slepc 158 
  MatStructure   flg;
6 dsic.upv.es!jroman 159 
 
160 
  PetscFunctionBegin;
161 
  /* Nothing to be done if STSetUp has not been called yet */
162 
  if (!st->setupcalled) PetscFunctionReturn(0);
2654 jroman 163 
 
1241 slepc 164 
  /* Check if the new KSP matrix has the same zero structure */
165 
  if (st->B && st->str == DIFFERENT_NONZERO_PATTERN && (st->sigma == 0.0 || newshift == 0.0)) {
166 
    flg = DIFFERENT_NONZERO_PATTERN;
167 
  } else {
168 
    flg = SAME_NONZERO_PATTERN;
169 
  }
6 dsic.upv.es!jroman 170 
 
344 dsic.upv.es!antodo 171 
  switch (st->shift_matrix) {
1940 jroman 172 
  case ST_MATMODE_INPLACE:
6 dsic.upv.es!jroman 173 
    /* Undo previous operations */
331 dsic.upv.es!antodo 174 
    if (st->sigma != 0.0) {
828 dsic.upv.es!antodo 175 
      if (st->B) {
176 
        ierr = MatAXPY(st->A,st->sigma,st->B,st->str);CHKERRQ(ierr);
177 
      } else {
178 
        ierr = MatShift(st->A,st->sigma);CHKERRQ(ierr);
179 
      }
331 dsic.upv.es!antodo 180 
    }
6 dsic.upv.es!jroman 181 
    /* Apply new shift */
331 dsic.upv.es!antodo 182 
    if (newshift != 0.0) {
828 dsic.upv.es!antodo 183 
      if (st->B) {
184 
        ierr = MatAXPY(st->A,-newshift,st->B,st->str);CHKERRQ(ierr);
185 
      } else {
186 
        ierr = MatShift(st->A,-newshift);CHKERRQ(ierr);
187 
      }
331 dsic.upv.es!antodo 188 
    }
1241 slepc 189 
    ierr = KSPSetOperators(st->ksp,st->A,st->A,flg);CHKERRQ(ierr);
163 dsic.upv.es!antodo 190 
    break;
1940 jroman 191 
  case ST_MATMODE_SHELL:
2654 jroman 192 
    ierr = KSPSetOperators(st->ksp,st->mat,st->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
163 dsic.upv.es!antodo 193 
    break;
194 
  default:
1241 slepc 195 
    if (st->mat) {
2654 jroman 196 
      ierr = MatCopy(st->A,st->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
1241 slepc 197 
    } else {
198 
      ierr = MatDuplicate(st->A,MAT_COPY_VALUES,&st->mat);CHKERRQ(ierr);
199 
    }
331 dsic.upv.es!antodo 200 
    if (newshift != 0.0) {  
2654 jroman 201 
      if (st->B) { ierr = MatAXPY(st->mat,-newshift,st->B,st->str);CHKERRQ(ierr); }
202 
      else { ierr = MatShift(st->mat,-newshift);CHKERRQ(ierr); }
331 dsic.upv.es!antodo 203 
    }
1241 slepc 204 
    ierr = KSPSetOperators(st->ksp,st->mat,st->mat,flg);CHKERRQ(ierr);    
6 dsic.upv.es!jroman 205 
  }
344 dsic.upv.es!antodo 206 
  st->sigma = newshift;
18 dsic.upv.es!jroman 207 
  ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
6 dsic.upv.es!jroman 208 
  PetscFunctionReturn(0);
209 
}
210 
 
2005 eromero 211 
#undef __FUNCT__  
212 
#define __FUNCT__ "STSetFromOptions_Sinvert"
213 
PetscErrorCode STSetFromOptions_Sinvert(ST st)
214 
{
215 
  PetscErrorCode ierr;
216 
  PC             pc;
217 
  const PCType   pctype;
218 
  const KSPType  ksptype;
219 
 
220 
  PetscFunctionBegin;
2364 jroman 221 
  if (!st->ksp) { ierr = STGetKSP(st,&st->ksp);CHKERRQ(ierr); }
2005 eromero 222 
  ierr = KSPGetPC(st->ksp,&pc);CHKERRQ(ierr);
223 
  ierr = KSPGetType(st->ksp,&ksptype);CHKERRQ(ierr);
224 
  ierr = PCGetType(pc,&pctype);CHKERRQ(ierr);
225 
  if (!pctype && !ksptype) {
226 
    if (st->shift_matrix == ST_MATMODE_SHELL) {
227 
      /* in shell mode use GMRES with Jacobi as the default */
228 
      ierr = KSPSetType(st->ksp,KSPGMRES);CHKERRQ(ierr);
229 
      ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
230 
    } else {
231 
      /* use direct solver as default */
232 
      ierr = KSPSetType(st->ksp,KSPPREONLY);CHKERRQ(ierr);
233 
      ierr = PCSetType(pc,PCREDUNDANT);CHKERRQ(ierr);
234 
    }
235 
  }
236 
  PetscFunctionReturn(0);
237 
}
238 
 
6 dsic.upv.es!jroman 239 
EXTERN_C_BEGIN
240 
#undef __FUNCT__  
241 
#define __FUNCT__ "STCreate_Sinvert"
476 dsic.upv.es!antodo 242 
PetscErrorCode STCreate_Sinvert(ST st)
6 dsic.upv.es!jroman 243 
{
244 
  PetscFunctionBegin;
1358 slepc 245 
  st->ops->apply           = STApply_Sinvert;
246 
  st->ops->getbilinearform = STGetBilinearForm_Default;
247 
  st->ops->applytrans      = STApplyTranspose_Sinvert;
248 
  st->ops->postsolve       = STPostSolve_Sinvert;
249 
  st->ops->backtr          = STBackTransform_Sinvert;
250 
  st->ops->setup           = STSetUp_Sinvert;
251 
  st->ops->setshift        = STSetShift_Sinvert;
2333 jroman 252 
  st->ops->setfromoptions  = STSetFromOptions_Sinvert;
253 
  st->ops->checknullspace  = STCheckNullSpace_Default;
6 dsic.upv.es!jroman 254 
  PetscFunctionReturn(0);
255 
}
256 
EXTERN_C_END
257