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