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 
      Implements the shift-and-invert technique for eigenvalue problems.
3 
*/
4 
#include "src/st/stimpl.h"          /*I "slepcst.h" I*/
5 
 
6 
#undef __FUNCT__  
7 
#define __FUNCT__ "STApply_Sinvert"
476 dsic.upv.es!antodo 8 
PetscErrorCode STApply_Sinvert(ST st,Vec x,Vec y)
6 dsic.upv.es!jroman 9 
{
476 dsic.upv.es!antodo 10 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 11 
 
12 
  PetscFunctionBegin;
13 
  if (st->B) {
14 
    /* generalized eigenproblem: y = (A - sB)^-1 B x */
344 dsic.upv.es!antodo 15 
    ierr = MatMult(st->B,x,st->w);CHKERRQ(ierr);
16 
    ierr = STAssociatedKSPSolve(st,st->w,y);CHKERRQ(ierr);
6 dsic.upv.es!jroman 17 
  }
18 
  else {
19 
    /* standard eigenproblem: y = (A - sI)^-1 x */
18 dsic.upv.es!jroman 20 
    ierr = STAssociatedKSPSolve(st,x,y);CHKERRQ(ierr);
6 dsic.upv.es!jroman 21 
  }
22 
  PetscFunctionReturn(0);
23 
}
24 
 
25 
#undef __FUNCT__  
26 
#define __FUNCT__ "STApplyNoB_Sinvert"
476 dsic.upv.es!antodo 27 
PetscErrorCode STApplyNoB_Sinvert(ST st,Vec x,Vec y)
6 dsic.upv.es!jroman 28 
{
476 dsic.upv.es!antodo 29 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 30 
 
31 
  PetscFunctionBegin;
18 dsic.upv.es!jroman 32 
  ierr = STAssociatedKSPSolve(st,x,y);CHKERRQ(ierr);
6 dsic.upv.es!jroman 33 
  PetscFunctionReturn(0);
34 
}
35 
 
36 
#undef __FUNCT__  
413 dsic.upv.es!antodo 37 
#define __FUNCT__ "STApplyB_Sinvert"
476 dsic.upv.es!antodo 38 
PetscErrorCode STApplyB_Sinvert(ST st,Vec x,Vec y)
413 dsic.upv.es!antodo 39 
{
476 dsic.upv.es!antodo 40 
  PetscErrorCode ierr;
413 dsic.upv.es!antodo 41 
 
42 
  PetscFunctionBegin;
43 
  if( st->B ) {
44 
    ierr = MatMult( st->B, x, y ); CHKERRQ(ierr);
45 
  }
46 
  else {
47 
    ierr = VecCopy( x, y ); CHKERRQ(ierr);
48 
  }
49 
  PetscFunctionReturn(0);
50 
}
51 
 
52 
#undef __FUNCT__  
6 dsic.upv.es!jroman 53 
#define __FUNCT__ "STBackTransform_Sinvert"
476 dsic.upv.es!antodo 54 
PetscErrorCode STBackTransform_Sinvert(ST st,PetscScalar *eigr,PetscScalar *eigi)
6 dsic.upv.es!jroman 55 
{
267 dsic.upv.es!antodo 56 
#ifndef PETSC_USE_COMPLEX
179 dsic.upv.es!antodo 57 
  PetscScalar t;
6 dsic.upv.es!jroman 58 
  PetscFunctionBegin;
179 dsic.upv.es!antodo 59 
  PetscValidPointer(eigr,2);
60 
  PetscValidPointer(eigi,3);
61 
  if (*eigi == 0) *eigr = 1.0 / *eigr + st->sigma;
62 
  else {
63 
    t = *eigr * *eigr + *eigi * *eigi;
64 
    *eigr = *eigr / t + st->sigma;
65 
    *eigi = - *eigi / t;
66 
  }
67 
#else
267 dsic.upv.es!antodo 68 
  PetscFunctionBegin;
179 dsic.upv.es!antodo 69 
  PetscValidPointer(eigr,2);
70 
  *eigr = 1.0 / *eigr + st->sigma;
71 
#endif
6 dsic.upv.es!jroman 72 
  PetscFunctionReturn(0);
73 
}
74 
 
75 
#undef __FUNCT__  
76 
#define __FUNCT__ "STPost_Sinvert"
476 dsic.upv.es!antodo 77 
PetscErrorCode STPost_Sinvert(ST st)
6 dsic.upv.es!jroman 78 
{
476 dsic.upv.es!antodo 79 
  PetscErrorCode ierr;
80 
  PetscScalar    alpha;
6 dsic.upv.es!jroman 81 
 
82 
  PetscFunctionBegin;
344 dsic.upv.es!antodo 83 
  if (st->shift_matrix == STMATMODE_INPLACE) {
6 dsic.upv.es!jroman 84 
    alpha = st->sigma;
344 dsic.upv.es!antodo 85 
    if( st->B ) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
6 dsic.upv.es!jroman 86 
    else { ierr = MatShift( &alpha, st->A ); CHKERRQ(ierr); }
87 
    st->setupcalled = 0;
88 
  }
89 
  PetscFunctionReturn(0);
90 
}
91 
 
92 
#undef __FUNCT__  
93 
#define __FUNCT__ "STSetUp_Sinvert"
476 dsic.upv.es!antodo 94 
PetscErrorCode STSetUp_Sinvert(ST st)
6 dsic.upv.es!jroman 95 
{
476 dsic.upv.es!antodo 96 
  PetscErrorCode ierr;
97 
  PetscScalar    alpha;
6 dsic.upv.es!jroman 98 
 
99 
  PetscFunctionBegin;
360 dsic.upv.es!antodo 100 
 
101 
  if (st->mat) { ierr = MatDestroy(st->mat);CHKERRQ(ierr); }
6 dsic.upv.es!jroman 102 
 
344 dsic.upv.es!antodo 103 
  switch (st->shift_matrix) {
104 
  case STMATMODE_INPLACE:
360 dsic.upv.es!antodo 105 
    st->mat = PETSC_NULL;
331 dsic.upv.es!antodo 106 
    if (st->sigma != 0.0) {
107 
      alpha = -st->sigma;
326 dsic.upv.es!antodo 108 
      if (st->B) {
344 dsic.upv.es!antodo 109 
        ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr);
326 dsic.upv.es!antodo 110 
      } else {
111 
        ierr = MatShift(&alpha,st->A);CHKERRQ(ierr);
112 
      }
113 
    }
15 dsic.upv.es!jroman 114 
    /* In the following line, the SAME_NONZERO_PATTERN flag has been used to
115 
     * improve performance when solving a number of related eigenproblems */
18 dsic.upv.es!jroman 116 
    ierr = KSPSetOperators(st->ksp,st->A,st->A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
163 dsic.upv.es!antodo 117 
    break;
344 dsic.upv.es!antodo 118 
  case STMATMODE_SHELL:
119 
    ierr = STMatShellCreate(st,&st->mat);CHKERRQ(ierr);
120 
    ierr = KSPSetOperators(st->ksp,st->mat,st->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
163 dsic.upv.es!antodo 121 
    break;
122 
  default:
344 dsic.upv.es!antodo 123 
    ierr = MatDuplicate(st->A,MAT_COPY_VALUES,&st->mat);CHKERRQ(ierr);
331 dsic.upv.es!antodo 124 
    if (st->sigma != 0.0) {
125 
      alpha = -st->sigma;
326 dsic.upv.es!antodo 126 
      if (st->B) {
344 dsic.upv.es!antodo 127 
        ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr);
326 dsic.upv.es!antodo 128 
      } else {
344 dsic.upv.es!antodo 129 
        ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr);
326 dsic.upv.es!antodo 130 
      }
131 
    }
163 dsic.upv.es!antodo 132 
    /* In the following line, the SAME_NONZERO_PATTERN flag has been used to
133 
     * improve performance when solving a number of related eigenproblems */
344 dsic.upv.es!antodo 134 
    ierr = KSPSetOperators(st->ksp,st->mat,st->mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
6 dsic.upv.es!jroman 135 
  }
18 dsic.upv.es!jroman 136 
  ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
6 dsic.upv.es!jroman 137 
  PetscFunctionReturn(0);
138 
}
139 
 
140 
#undef __FUNCT__  
141 
#define __FUNCT__ "STSetShift_Sinvert"
476 dsic.upv.es!antodo 142 
PetscErrorCode STSetShift_Sinvert(ST st,PetscScalar newshift)
6 dsic.upv.es!jroman 143 
{
476 dsic.upv.es!antodo 144 
  PetscErrorCode ierr;
145 
  PetscScalar    alpha;
6 dsic.upv.es!jroman 146 
 
147 
  PetscFunctionBegin;
148 
 
149 
  /* Nothing to be done if STSetUp has not been called yet */
150 
  if (!st->setupcalled) PetscFunctionReturn(0);
151 
 
344 dsic.upv.es!antodo 152 
  switch (st->shift_matrix) {
153 
  case STMATMODE_INPLACE:
6 dsic.upv.es!jroman 154 
    /* Undo previous operations */
331 dsic.upv.es!antodo 155 
    if (st->sigma != 0.0) {
156 
      alpha = st->sigma;
344 dsic.upv.es!antodo 157 
      if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
331 dsic.upv.es!antodo 158 
      else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
159 
    }
6 dsic.upv.es!jroman 160 
    /* Apply new shift */
331 dsic.upv.es!antodo 161 
    if (newshift != 0.0) {
162 
      alpha = -newshift;
344 dsic.upv.es!antodo 163 
      if (st->B) { ierr = MatAXPY(&alpha,st->B,st->A,st->str);CHKERRQ(ierr); }
331 dsic.upv.es!antodo 164 
      else { ierr = MatShift(&alpha,st->A);CHKERRQ(ierr); }
165 
    }
18 dsic.upv.es!jroman 166 
    ierr = KSPSetOperators(st->ksp,st->A,st->A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
163 dsic.upv.es!antodo 167 
    break;
344 dsic.upv.es!antodo 168 
  case STMATMODE_SHELL:
354 dsic.upv.es!jroman 169 
    ierr = KSPSetOperators(st->ksp,st->mat,st->mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);    
163 dsic.upv.es!antodo 170 
    break;
171 
  default:
354 dsic.upv.es!jroman 172 
    ierr = MatCopy(st->A, st->mat,SUBSET_NONZERO_PATTERN); CHKERRQ(ierr);
331 dsic.upv.es!antodo 173 
    if (newshift != 0.0) {  
174 
      alpha = -newshift;
344 dsic.upv.es!antodo 175 
      if (st->B) { ierr = MatAXPY(&alpha,st->B,st->mat,st->str);CHKERRQ(ierr); }
176 
      else { ierr = MatShift(&alpha,st->mat);CHKERRQ(ierr); }
331 dsic.upv.es!antodo 177 
    }
163 dsic.upv.es!antodo 178 
    /* In the following line, the SAME_NONZERO_PATTERN flag has been used to
179 
     * improve performance when solving a number of related eigenproblems */
344 dsic.upv.es!antodo 180 
    ierr = KSPSetOperators(st->ksp,st->mat,st->mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);    
6 dsic.upv.es!jroman 181 
  }
344 dsic.upv.es!antodo 182 
  st->sigma = newshift;
18 dsic.upv.es!jroman 183 
  ierr = KSPSetUp(st->ksp);CHKERRQ(ierr);
6 dsic.upv.es!jroman 184 
  PetscFunctionReturn(0);
185 
}
186 
 
187 
EXTERN_C_BEGIN
188 
#undef __FUNCT__  
189 
#define __FUNCT__ "STCreate_Sinvert"
476 dsic.upv.es!antodo 190 
PetscErrorCode STCreate_Sinvert(ST st)
6 dsic.upv.es!jroman 191 
{
192 
  PetscFunctionBegin;
344 dsic.upv.es!antodo 193 
  st->data                = 0;
6 dsic.upv.es!jroman 194 
 
195 
  st->ops->apply          = STApply_Sinvert;
413 dsic.upv.es!antodo 196 
  st->ops->applyB         = STApplyB_Sinvert;
6 dsic.upv.es!jroman 197 
  st->ops->applynoB       = STApplyNoB_Sinvert;
198 
  st->ops->postsolve      = STPost_Sinvert;
199 
  st->ops->backtr         = STBackTransform_Sinvert;
200 
  st->ops->setup          = STSetUp_Sinvert;
201 
  st->ops->setshift       = STSetShift_Sinvert;
438 dsic.upv.es!antodo 202 
  st->ops->view           = STView_Default;
310 dsic.upv.es!antodo 203 
 
204 
  st->checknullspace      = STCheckNullSpace_Default;
6 dsic.upv.es!jroman 205 
 
206 
  PetscFunctionReturn(0);
207 
}
208 
EXTERN_C_END
209