Subversion Repositories slepc-dev

Rev

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

Rev Author Line No. Line
521 dsic.upv.es!antodo 1 
/*
545 dsic.upv.es!jroman 2 
     This file contains some simple default routines for common operations.  
3 
*/
521 dsic.upv.es!antodo 4 
#include "src/eps/epsimpl.h"   /*I "slepceps.h" I*/
5 
#include "slepcblaslapack.h"
6 
 
7 
#undef __FUNCT__  
8 
#define __FUNCT__ "EPSDestroy_Default"
9 
PetscErrorCode EPSDestroy_Default(EPS eps)
10 
{
11 
  PetscErrorCode ierr;
12 
 
13 
  PetscFunctionBegin;
14 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
1040 slepc 15 
  ierr = PetscFree(eps->data);CHKERRQ(ierr);
521 dsic.upv.es!antodo 16 
 
17 
  /* free work vectors */
18 
  ierr = EPSDefaultFreeWork(eps);CHKERRQ(ierr);
19 
  ierr = EPSFreeSolution(eps);CHKERRQ(ierr);
20 
  PetscFunctionReturn(0);
21 
}
22 
 
23 
#undef __FUNCT__  
24 
#define __FUNCT__ "EPSBackTransform_Default"
25 
PetscErrorCode EPSBackTransform_Default(EPS eps)
26 
{
27 
  PetscErrorCode ierr;
28 
  int            i;
29 
 
30 
  PetscFunctionBegin;
31 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
32 
  for (i=0;i<eps->nconv;i++) {
33 
    ierr = STBackTransform(eps->OP,&eps->eigr[i],&eps->eigi[i]);CHKERRQ(ierr);
34 
  }
35 
  PetscFunctionReturn(0);
36 
}
37 
 
38 
#undef __FUNCT__  
39 
#define __FUNCT__ "EPSComputeVectors_Default"
545 dsic.upv.es!jroman 40 
/*
41 
  EPSComputeVectors_Default - Compute eigenvectors from the vectors
42 
  provided by the eigensolver. This version just copies the vectors
43 
  and is intended for solvers such as power that provide the eigenvector.
44 
 */
521 dsic.upv.es!antodo 45 
PetscErrorCode EPSComputeVectors_Default(EPS eps)
46 
{
47 
  PetscErrorCode ierr;
48 
  int            i;
49 
 
50 
  PetscFunctionBegin;
51 
  for (i=0;i<eps->nconv;i++) {
52 
    ierr = VecCopy(eps->V[i],eps->AV[i]);CHKERRQ(ierr);
780 dsic.upv.es!jroman 53 
    if (eps->solverclass == EPS_TWO_SIDE) {
54 
      ierr = VecCopy(eps->W[i],eps->AW[i]);CHKERRQ(ierr);
55 
    }
521 dsic.upv.es!antodo 56 
  }
57 
  eps->evecsavailable = PETSC_TRUE;
58 
  PetscFunctionReturn(0);
59 
}
60 
 
61 
#undef __FUNCT__  
1243 slepc 62 
#define __FUNCT__ "EPSComputeVectors_Hermitian"
63 
/*
64 
  EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
65 
  using purification for generalized eigenproblems.
66 
 */
67 
PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
68 
{
69 
  PetscErrorCode ierr;
70 
  int            i;
71 
  PetscReal      norm;
72 
 
73 
  PetscFunctionBegin;
74 
  for (i=0;i<eps->nconv;i++) {
75 
    if (eps->isgeneralized) {
76 
      /* Purify eigenvectors */
77 
      ierr = STApply(eps->OP,eps->V[i],eps->AV[i]);CHKERRQ(ierr);
78 
      ierr = VecNormalize(eps->AV[i],&norm);CHKERRQ(ierr);
79 
    } else {
80 
      ierr = VecCopy(eps->V[i],eps->AV[i]);CHKERRQ(ierr);  
81 
    }
82 
    if (eps->solverclass == EPS_TWO_SIDE) {
83 
      ierr = VecCopy(eps->W[i],eps->AW[i]);CHKERRQ(ierr);
84 
    }
85 
  }
86 
  eps->evecsavailable = PETSC_TRUE;
87 
  PetscFunctionReturn(0);
88 
}
89 
 
90 
#undef __FUNCT__  
521 dsic.upv.es!antodo 91 
#define __FUNCT__ "EPSComputeVectors_Schur"
545 dsic.upv.es!jroman 92 
/*
93 
  EPSComputeVectors_Schur - Compute eigenvectors from the vectors
94 
  provided by the eigensolver. This version is intended for solvers
95 
  that provide Schur vectors. Given the partial Schur decomposition
96 
  OP*V=V*T, the following steps are performed:
780 dsic.upv.es!jroman 97 
      1) compute eigenvectors of T: T*Z=Z*D
98 
      2) compute eigenvectors of OP: X=V*Z
99 
  If left eigenvectors are required then also do Z'*Tl=D*Z', Y=W*Z
545 dsic.upv.es!jroman 100 
 */
521 dsic.upv.es!antodo 101 
PetscErrorCode EPSComputeVectors_Schur(EPS eps)
102 
{
806 dsic.upv.es!antodo 103 
#if defined(SLEPC_MISSING_LAPACK_TREVC)
104 
  SETERRQ(PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable.");
105 
#else
521 dsic.upv.es!antodo 106 
  PetscErrorCode ierr;
1057 slepc 107 
  int            i,mout,info,nv=eps->nv;
780 dsic.upv.es!jroman 108 
  PetscScalar    *Z,*work;
521 dsic.upv.es!antodo 109 
#if defined(PETSC_USE_COMPLEX)
110 
  PetscReal      *rwork;
111 
#endif
1345 slepc 112 
  IPBilinearForm form;
1243 slepc 113 
  PetscReal      norm;
114 
  Vec            w;
1351 slepc 115 
  Mat            B;
521 dsic.upv.es!antodo 116 
 
117 
  PetscFunctionBegin;
118 
  if (eps->ishermitian) {
1243 slepc 119 
    ierr = EPSComputeVectors_Hermitian(eps);CHKERRQ(ierr);
521 dsic.upv.es!antodo 120 
    PetscFunctionReturn(0);
121 
  }
1243 slepc 122 
 
1351 slepc 123 
  ierr = IPGetBilinearForm(eps->ip,&B,&form);CHKERRQ(ierr);
124 
  if (B && form == IPINNER_HERMITIAN) {
1243 slepc 125 
    ierr = VecDuplicate(eps->V[0],&w);CHKERRQ(ierr);
126 
  }
521 dsic.upv.es!antodo 127 
 
1057 slepc 128 
  ierr = PetscMalloc(nv*nv*sizeof(PetscScalar),&Z);CHKERRQ(ierr);
129 
  ierr = PetscMalloc(3*nv*sizeof(PetscScalar),&work);CHKERRQ(ierr);
521 dsic.upv.es!antodo 130 
#if defined(PETSC_USE_COMPLEX)
1057 slepc 131 
  ierr = PetscMalloc(nv*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
521 dsic.upv.es!antodo 132 
#endif
133 
 
780 dsic.upv.es!jroman 134 
  /* right eigenvectors */
521 dsic.upv.es!antodo 135 
#if !defined(PETSC_USE_COMPLEX)
1057 slepc 136 
  LAPACKtrevc_("R","A",PETSC_NULL,&nv,eps->T,&eps->ncv,PETSC_NULL,&nv,Z,&nv,&nv,&mout,work,&info,1,1);
521 dsic.upv.es!antodo 137 
#else
1057 slepc 138 
  LAPACKtrevc_("R","A",PETSC_NULL,&nv,eps->T,&eps->ncv,PETSC_NULL,&nv,Z,&nv,&nv,&mout,work,rwork,&info,1,1);
521 dsic.upv.es!antodo 139 
#endif
140 
  if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xTREVC %i",info);
141 
 
870 dsic.upv.es!antodo 142 
  /* AV = V * Z */
521 dsic.upv.es!antodo 143 
  for (i=0;i<eps->nconv;i++) {
1351 slepc 144 
    if (B && form == IPINNER_HERMITIAN) {
1243 slepc 145 
      /* Purify eigenvectors */
146 
      ierr = VecSet(w,0.0);CHKERRQ(ierr);
147 
      ierr = VecMAXPY(w,nv,Z+nv*i,eps->V);CHKERRQ(ierr);
148 
      ierr = STApply(eps->OP,w,eps->AV[i]);CHKERRQ(ierr);
149 
      ierr = VecNormalize(eps->AV[i],&norm);CHKERRQ(ierr);
150 
    } else {
151 
      ierr = VecSet(eps->AV[i],0.0);CHKERRQ(ierr);
152 
      ierr = VecMAXPY(eps->AV[i],nv,Z+nv*i,eps->V);CHKERRQ(ierr);
153 
    }
870 dsic.upv.es!antodo 154 
  }    
521 dsic.upv.es!antodo 155 
 
780 dsic.upv.es!jroman 156 
  /* left eigenvectors */
157 
  if (eps->solverclass == EPS_TWO_SIDE) {
158 
#if !defined(PETSC_USE_COMPLEX)
1057 slepc 159 
    LAPACKtrevc_("R","A",PETSC_NULL,&nv,eps->Tl,&eps->nv,PETSC_NULL,&nv,Z,&nv,&nv,&mout,work,&info,1,1);
780 dsic.upv.es!jroman 160 
#else
1057 slepc 161 
    LAPACKtrevc_("R","A",PETSC_NULL,&nv,eps->Tl,&eps->nv,PETSC_NULL,&nv,Z,&nv,&nv,&mout,work,rwork,&info,1,1);
780 dsic.upv.es!jroman 162 
#endif
163 
    if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xTREVC %i",info);
164 
 
870 dsic.upv.es!antodo 165 
    /* AW = W * Z */
780 dsic.upv.es!jroman 166 
    for (i=0;i<eps->nconv;i++) {
870 dsic.upv.es!antodo 167 
      ierr = VecSet(eps->AW[i],0.0);CHKERRQ(ierr);
1057 slepc 168 
      ierr = VecMAXPY(eps->AW[i],nv,Z+nv*i,eps->W);CHKERRQ(ierr);
870 dsic.upv.es!antodo 169 
    }    
780 dsic.upv.es!jroman 170 
  }
171 
 
172 
  ierr = PetscFree(Z);CHKERRQ(ierr);
521 dsic.upv.es!antodo 173 
  ierr = PetscFree(work);CHKERRQ(ierr);
174 
#if defined(PETSC_USE_COMPLEX)
175 
  ierr = PetscFree(rwork);CHKERRQ(ierr);
176 
#endif
1351 slepc 177 
  if (B && form == IPINNER_HERMITIAN) {
1243 slepc 178 
    ierr = VecDestroy(w);CHKERRQ(ierr);
179 
  }
521 dsic.upv.es!antodo 180 
  eps->evecsavailable = PETSC_TRUE;
181 
  PetscFunctionReturn(0);
806 dsic.upv.es!antodo 182 
#endif
521 dsic.upv.es!antodo 183 
}
533 dsic.upv.es!antodo 184 
 
185 
#undef __FUNCT__  
186 
#define __FUNCT__ "EPSDefaultGetWork"
187 
/*
188 
  EPSDefaultGetWork - Gets a number of work vectors.
189 
 
190 
  Input Parameters:
191 
+ eps  - eigensolver context
192 
- nw   - number of work vectors to allocate
193 
 
194 
  Notes:
195 
  Call this only if no work vectors have been allocated.
196 
 
197 
 */
198 
PetscErrorCode EPSDefaultGetWork(EPS eps, int nw)
199 
{
200 
  PetscErrorCode ierr;
201 
 
202 
  PetscFunctionBegin;
203 
 
204 
  if (eps->nwork != nw) {
205 
    if (eps->nwork > 0) {
206 
      ierr = VecDestroyVecs(eps->work,eps->nwork); CHKERRQ(ierr);
207 
    }
208 
    eps->nwork = nw;
1343 slepc 209 
    ierr = VecDuplicateVecs(eps->IV[0],nw,&eps->work); CHKERRQ(ierr);
790 dsic.upv.es!antodo 210 
    ierr = PetscLogObjectParents(eps,nw,eps->work);
533 dsic.upv.es!antodo 211 
  }
212 
 
213 
  PetscFunctionReturn(0);
214 
}
215 
 
216 
#undef __FUNCT__  
217 
#define __FUNCT__ "EPSDefaultFreeWork"
218 
/*
219 
  EPSDefaultFreeWork - Free work vectors.
220 
 
221 
  Input Parameters:
222 
. eps  - eigensolver context
223 
 
224 
 */
225 
PetscErrorCode EPSDefaultFreeWork(EPS eps)
226 
{
227 
  PetscErrorCode ierr;
228 
 
229 
  PetscFunctionBegin;
230 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
231 
  if (eps->work)  {
232 
    ierr = VecDestroyVecs(eps->work,eps->nwork); CHKERRQ(ierr);
233 
  }
234 
  PetscFunctionReturn(0);
235 
}