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
1887 jroman 1 
/*                      
2 
 
3 
   Straightforward linearization for quadratic eigenproblems.
4 
 
5 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
6 
   SLEPc - Scalable Library for Eigenvalue Problem Computations
2116 eromero 7 
   Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
1887 jroman 8 
 
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/>.
22 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
23 
*/
24 
 
2284 jroman 25 
#include <private/qepimpl.h>         /*I "slepcqep.h" I*/
26 
#include <private/epsimpl.h>         /*I "slepceps.h" I*/
1908 jroman 27 
#include "linearp.h"
1887 jroman 28 
 
29 
#undef __FUNCT__  
2324 jroman 30 
#define __FUNCT__ "QEPSetUp_Linear"
31 
PetscErrorCode QEPSetUp_Linear(QEP qep)
1887 jroman 32 
{
2317 jroman 33 
  PetscErrorCode ierr;
34 
  QEP_LINEAR     *ctx = (QEP_LINEAR *)qep->data;
35 
  PetscInt       i=0;
36 
  EPSWhich       which;
37 
  PetscBool      trackall;
1908 jroman 38 
  /* function tables */
39 
  PetscErrorCode (*fcreate[][2])(MPI_Comm,QEP_LINEAR*,Mat*) = {
2324 jroman 40 
    { MatCreateExplicit_Linear_N1A, MatCreateExplicit_Linear_N1B },   /* N1 */
41 
    { MatCreateExplicit_Linear_N2A, MatCreateExplicit_Linear_N2B },   /* N2 */
42 
    { MatCreateExplicit_Linear_S1A, MatCreateExplicit_Linear_S1B },   /* S1 */
43 
    { MatCreateExplicit_Linear_S2A, MatCreateExplicit_Linear_S2B },   /* S2 */
44 
    { MatCreateExplicit_Linear_H1A, MatCreateExplicit_Linear_H1B },   /* H1 */
45 
    { MatCreateExplicit_Linear_H2A, MatCreateExplicit_Linear_H2B }    /* H2 */
1908 jroman 46 
  };
47 
  PetscErrorCode (*fmult[][2])(Mat,Vec,Vec) = {
2324 jroman 48 
    { MatMult_Linear_N1A, MatMult_Linear_N1B },
49 
    { MatMult_Linear_N2A, MatMult_Linear_N2B },
50 
    { MatMult_Linear_S1A, MatMult_Linear_S1B },
51 
    { MatMult_Linear_S2A, MatMult_Linear_S2B },
52 
    { MatMult_Linear_H1A, MatMult_Linear_H1B },
53 
    { MatMult_Linear_H2A, MatMult_Linear_H2B }
1908 jroman 54 
  };
55 
  PetscErrorCode (*fgetdiagonal[][2])(Mat,Vec) = {
2324 jroman 56 
    { MatGetDiagonal_Linear_N1A, MatGetDiagonal_Linear_N1B },
57 
    { MatGetDiagonal_Linear_N2A, MatGetDiagonal_Linear_N2B },
58 
    { MatGetDiagonal_Linear_S1A, MatGetDiagonal_Linear_S1B },
59 
    { MatGetDiagonal_Linear_S2A, MatGetDiagonal_Linear_S2B },
60 
    { MatGetDiagonal_Linear_H1A, MatGetDiagonal_Linear_H1B },
61 
    { MatGetDiagonal_Linear_H2A, MatGetDiagonal_Linear_H2B }
1908 jroman 62 
  };
1891 jroman 63 
 
64 
  PetscFunctionBegin;
1942 jroman 65 
  if (!qep->which) qep->which = QEP_LARGEST_MAGNITUDE;
1919 jroman 66 
  ctx->M = qep->M;
67 
  ctx->C = qep->C;
68 
  ctx->K = qep->K;
69 
  ctx->sfactor = qep->sfactor;
70 
 
2305 jroman 71 
  ierr = MatDestroy(&ctx->A);CHKERRQ(ierr);
72 
  ierr = MatDestroy(&ctx->B);CHKERRQ(ierr);
73 
  ierr = VecDestroy(&ctx->x1);CHKERRQ(ierr);
74 
  ierr = VecDestroy(&ctx->x2);CHKERRQ(ierr);
75 
  ierr = VecDestroy(&ctx->y1);CHKERRQ(ierr);
76 
  ierr = VecDestroy(&ctx->y2);CHKERRQ(ierr);
1891 jroman 77 
 
1910 jroman 78 
  switch (qep->problem_type) {
79 
    case QEP_GENERAL:    i = 0; break;
1911 jroman 80 
    case QEP_HERMITIAN:  i = 2; break;
1913 jroman 81 
    case QEP_GYROSCOPIC: i = 4; break;
2214 jroman 82 
    default: SETERRQ(((PetscObject)qep)->comm,1,"Wrong value of qep->problem_type");
1910 jroman 83 
  }
84 
  i += ctx->cform-1;
1908 jroman 85 
 
1891 jroman 86 
  if (ctx->explicitmatrix) {
87 
    ctx->x1 = ctx->x2 = ctx->y1 = ctx->y2 = PETSC_NULL;
1910 jroman 88 
    ierr = (*fcreate[i][0])(((PetscObject)qep)->comm,ctx,&ctx->A);CHKERRQ(ierr);
89 
    ierr = (*fcreate[i][1])(((PetscObject)qep)->comm,ctx,&ctx->B);CHKERRQ(ierr);
1891 jroman 90 
  } else {
1930 jroman 91 
    ierr = VecCreateMPIWithArray(((PetscObject)qep)->comm,qep->nloc,qep->n,PETSC_NULL,&ctx->x1);CHKERRQ(ierr);
92 
    ierr = VecCreateMPIWithArray(((PetscObject)qep)->comm,qep->nloc,qep->n,PETSC_NULL,&ctx->x2);CHKERRQ(ierr);
93 
    ierr = VecCreateMPIWithArray(((PetscObject)qep)->comm,qep->nloc,qep->n,PETSC_NULL,&ctx->y1);CHKERRQ(ierr);
94 
    ierr = VecCreateMPIWithArray(((PetscObject)qep)->comm,qep->nloc,qep->n,PETSC_NULL,&ctx->y2);CHKERRQ(ierr);
95 
    ierr = MatCreateShell(((PetscObject)qep)->comm,2*qep->nloc,2*qep->nloc,2*qep->n,2*qep->n,ctx,&ctx->A);CHKERRQ(ierr);
1910 jroman 96 
    ierr = MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))fmult[i][0]);CHKERRQ(ierr);
97 
    ierr = MatShellSetOperation(ctx->A,MATOP_GET_DIAGONAL,(void(*)(void))fgetdiagonal[i][0]);CHKERRQ(ierr);
1930 jroman 98 
    ierr = MatCreateShell(((PetscObject)qep)->comm,2*qep->nloc,2*qep->nloc,2*qep->n,2*qep->n,ctx,&ctx->B);CHKERRQ(ierr);
1910 jroman 99 
    ierr = MatShellSetOperation(ctx->B,MATOP_MULT,(void(*)(void))fmult[i][1]);CHKERRQ(ierr);
100 
    ierr = MatShellSetOperation(ctx->B,MATOP_GET_DIAGONAL,(void(*)(void))fgetdiagonal[i][1]);CHKERRQ(ierr);
1891 jroman 101 
  }
102 
 
103 
  ierr = EPSSetOperators(ctx->eps,ctx->A,ctx->B);CHKERRQ(ierr);
104 
  ierr = EPSSetProblemType(ctx->eps,EPS_GNHEP);CHKERRQ(ierr);
105 
  switch (qep->which) {
106 
      case QEP_LARGEST_MAGNITUDE:  which = EPS_LARGEST_MAGNITUDE; break;
107 
      case QEP_SMALLEST_MAGNITUDE: which = EPS_SMALLEST_MAGNITUDE; break;
108 
      case QEP_LARGEST_REAL:       which = EPS_LARGEST_REAL; break;
109 
      case QEP_SMALLEST_REAL:      which = EPS_SMALLEST_REAL; break;
110 
      case QEP_LARGEST_IMAGINARY:  which = EPS_LARGEST_IMAGINARY; break;
111 
      case QEP_SMALLEST_IMAGINARY: which = EPS_SMALLEST_IMAGINARY; break;
2214 jroman 112 
      default: SETERRQ(((PetscObject)qep)->comm,1,"Wrong value of which");
1891 jroman 113 
  }
114 
  ierr = EPSSetWhichEigenpairs(ctx->eps,which);CHKERRQ(ierr);
1944 jroman 115 
  ierr = EPSSetLeftVectorsWanted(ctx->eps,qep->leftvecs);CHKERRQ(ierr);
1891 jroman 116 
  ierr = EPSSetDimensions(ctx->eps,qep->nev,qep->ncv,qep->mpd);CHKERRQ(ierr);
117 
  ierr = EPSSetTolerances(ctx->eps,qep->tol,qep->max_it);CHKERRQ(ierr);
2055 eromero 118 
  /* Transfer the trackall option from qep to eps */
119 
  ierr = QEPGetTrackAll(qep,&trackall);CHKERRQ(ierr);
120 
  ierr = EPSSetTrackAll(ctx->eps,trackall);CHKERRQ(ierr);
2177 jroman 121 
  if (ctx->setfromoptionscalled) {
2043 eromero 122 
    ierr = EPSSetFromOptions(ctx->eps);CHKERRQ(ierr);
123 
    ctx->setfromoptionscalled = PETSC_FALSE;
124 
  }
1891 jroman 125 
  ierr = EPSSetUp(ctx->eps);CHKERRQ(ierr);
126 
  ierr = EPSGetDimensions(ctx->eps,PETSC_NULL,&qep->ncv,&qep->mpd);CHKERRQ(ierr);
127 
  ierr = EPSGetTolerances(ctx->eps,&qep->tol,&qep->max_it);CHKERRQ(ierr);
1952 jroman 128 
  if (qep->nini>0 || qep->ninil>0) PetscInfo(qep,"Ignoring initial vectors\n");
1887 jroman 129 
  PetscFunctionReturn(0);
130 
}
131 
 
132 
#undef __FUNCT__  
2118 jroman 133 
#define __FUNCT__ "QEPLinearSelect_Norm"
1902 jroman 134 
/*
2118 jroman 135 
   QEPLinearSelect_Norm - Auxiliary routine that copies the solution of the
1902 jroman 136 
   linear eigenproblem to the QEP object. The eigenvector of the generalized
137 
   problem is supposed to be
138 
                               z = [  x  ]
139 
                                   [ l*x ]
2118 jroman 140 
   The eigenvector is taken from z(1:n) or z(n+1:2*n) depending on the explicitly
141 
   computed residual norm.
142 
   Finally, x is normalized so that ||x||_2 = 1.
143 
*/
2295 jroman 144 
PetscErrorCode QEPLinearSelect_Norm(QEP qep,EPS eps)
2118 jroman 145 
{
146 
  PetscErrorCode ierr;
2295 jroman 147 
  PetscInt       i;
2152 eromero 148 
  PetscScalar    *px;
2118 jroman 149 
  PetscReal      rn1,rn2;
2295 jroman 150 
  Vec            xr,xi,wr,wi;
2118 jroman 151 
  Mat            A;
2152 eromero 152 
#if !defined(PETSC_USE_COMPLEX)
153 
  PetscScalar    *py;
154 
#endif
2118 jroman 155 
 
156 
  PetscFunctionBegin;
157 
  ierr = EPSGetOperators(eps,&A,PETSC_NULL);CHKERRQ(ierr);
2295 jroman 158 
  ierr = MatGetVecs(A,&xr,PETSC_NULL);CHKERRQ(ierr);
159 
  ierr = VecDuplicate(xr,&xi);CHKERRQ(ierr);
160 
  ierr = VecCreateMPIWithArray(((PetscObject)qep)->comm,qep->nloc,qep->n,PETSC_NULL,&wr);CHKERRQ(ierr);
161 
  ierr = VecCreateMPIWithArray(((PetscObject)qep)->comm,qep->nloc,qep->n,PETSC_NULL,&wi);CHKERRQ(ierr);
162 
  for (i=0;i<qep->nconv;i++) {
163 
    ierr = EPSGetEigenpair(eps,i,&qep->eigr[i],&qep->eigi[i],xr,xi);CHKERRQ(ierr);
164 
    qep->eigr[i] *= qep->sfactor;
165 
    qep->eigi[i] *= qep->sfactor;
2118 jroman 166 
#if !defined(PETSC_USE_COMPLEX)
2295 jroman 167 
    if (qep->eigi[i]>0.0) {   /* first eigenvalue of a complex conjugate pair */
168 
      ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
169 
      ierr = VecGetArray(xi,&py);CHKERRQ(ierr);
170 
      ierr = VecPlaceArray(wr,px);CHKERRQ(ierr);
171 
      ierr = VecPlaceArray(wi,py);CHKERRQ(ierr);
172 
      ierr = SlepcVecNormalize(wr,wi,PETSC_TRUE,PETSC_NULL);CHKERRQ(ierr);
173 
      ierr = QEPComputeResidualNorm_Private(qep,qep->eigr[i],qep->eigi[i],wr,wi,&rn1);CHKERRQ(ierr);
174 
      ierr = VecCopy(wr,qep->V[i]);CHKERRQ(ierr);
175 
      ierr = VecCopy(wi,qep->V[i+1]);CHKERRQ(ierr);
176 
      ierr = VecResetArray(wr);CHKERRQ(ierr);
177 
      ierr = VecResetArray(wi);CHKERRQ(ierr);
178 
      ierr = VecPlaceArray(wr,px+qep->nloc);CHKERRQ(ierr);
179 
      ierr = VecPlaceArray(wi,py+qep->nloc);CHKERRQ(ierr);
180 
      ierr = SlepcVecNormalize(wr,wi,PETSC_TRUE,PETSC_NULL);CHKERRQ(ierr);
181 
      ierr = QEPComputeResidualNorm_Private(qep,qep->eigr[i],qep->eigi[i],wr,wi,&rn2);CHKERRQ(ierr);
182 
      if (rn1>rn2) {
2118 jroman 183 
        ierr = VecCopy(wr,qep->V[i]);CHKERRQ(ierr);
184 
        ierr = VecCopy(wi,qep->V[i+1]);CHKERRQ(ierr);
185 
      }
2295 jroman 186 
      ierr = VecResetArray(wr);CHKERRQ(ierr);
187 
      ierr = VecResetArray(wi);CHKERRQ(ierr);
188 
      ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
189 
      ierr = VecRestoreArray(xi,&py);CHKERRQ(ierr);
2118 jroman 190 
    }
2295 jroman 191 
    else if (qep->eigi[i]==0.0)   /* real eigenvalue */
2118 jroman 192 
#endif
2295 jroman 193 
    {
194 
      ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
195 
      ierr = VecPlaceArray(wr,px);CHKERRQ(ierr);
196 
      ierr = SlepcVecNormalize(wr,PETSC_NULL,PETSC_FALSE,PETSC_NULL);CHKERRQ(ierr);
197 
      ierr = QEPComputeResidualNorm_Private(qep,qep->eigr[i],qep->eigi[i],wr,PETSC_NULL,&rn1);CHKERRQ(ierr);
198 
      ierr = VecCopy(wr,qep->V[i]);CHKERRQ(ierr);
199 
      ierr = VecResetArray(wr);CHKERRQ(ierr);
200 
      ierr = VecPlaceArray(wr,px+qep->nloc);CHKERRQ(ierr);
201 
      ierr = SlepcVecNormalize(wr,PETSC_NULL,PETSC_FALSE,PETSC_NULL);CHKERRQ(ierr);
202 
      ierr = QEPComputeResidualNorm_Private(qep,qep->eigr[i],qep->eigi[i],wr,PETSC_NULL,&rn2);CHKERRQ(ierr);
203 
      if (rn1>rn2) {
2118 jroman 204 
        ierr = VecCopy(wr,qep->V[i]);CHKERRQ(ierr);
205 
      }
2295 jroman 206 
      ierr = VecResetArray(wr);CHKERRQ(ierr);
207 
      ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
2118 jroman 208 
    }
209 
  }
2305 jroman 210 
  ierr = VecDestroy(&wr);CHKERRQ(ierr);
211 
  ierr = VecDestroy(&wi);CHKERRQ(ierr);
212 
  ierr = VecDestroy(&xr);CHKERRQ(ierr);
213 
  ierr = VecDestroy(&xi);CHKERRQ(ierr);
2118 jroman 214 
  PetscFunctionReturn(0);
215 
}
216 
 
217 
#undef __FUNCT__  
218 
#define __FUNCT__ "QEPLinearSelect_Simple"
219 
/*
220 
   QEPLinearSelect_Simple - Auxiliary routine that copies the solution of the
221 
   linear eigenproblem to the QEP object. The eigenvector of the generalized
222 
   problem is supposed to be
223 
                               z = [  x  ]
224 
                                   [ l*x ]
1902 jroman 225 
   If |l|<1.0, the eigenvector is taken from z(1:n), otherwise from z(n+1:2*n).
226 
   Finally, x is normalized so that ||x||_2 = 1.
227 
*/
2295 jroman 228 
PetscErrorCode QEPLinearSelect_Simple(QEP qep,EPS eps)
1887 jroman 229 
{
1891 jroman 230 
  PetscErrorCode ierr;
2295 jroman 231 
  PetscInt       i,offset;
1902 jroman 232 
  PetscScalar    *px;
2295 jroman 233 
  Vec            xr,xi,w;
1933 jroman 234 
  Mat            A;
1902 jroman 235 
 
236 
  PetscFunctionBegin;
1933 jroman 237 
  ierr = EPSGetOperators(eps,&A,PETSC_NULL);CHKERRQ(ierr);
2295 jroman 238 
  ierr = MatGetVecs(A,&xr,PETSC_NULL);CHKERRQ(ierr);
239 
  ierr = VecDuplicate(xr,&xi);CHKERRQ(ierr);
240 
  ierr = VecCreateMPIWithArray(((PetscObject)qep)->comm,qep->nloc,qep->n,PETSC_NULL,&w);CHKERRQ(ierr);
241 
  for (i=0;i<qep->nconv;i++) {
242 
    ierr = EPSGetEigenpair(eps,i,&qep->eigr[i],&qep->eigi[i],xr,xi);CHKERRQ(ierr);
243 
    qep->eigr[i] *= qep->sfactor;
244 
    qep->eigi[i] *= qep->sfactor;
245 
    if (SlepcAbsEigenvalue(qep->eigr[i],qep->eigi[i])>1.0) offset = qep->nloc;
246 
    else offset = 0;
1894 jroman 247 
#if !defined(PETSC_USE_COMPLEX)
2295 jroman 248 
    if (qep->eigi[i]>0.0) {   /* first eigenvalue of a complex conjugate pair */
249 
      ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
250 
      ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
251 
      ierr = VecCopy(w,qep->V[i]);CHKERRQ(ierr);
252 
      ierr = VecResetArray(w);CHKERRQ(ierr);
253 
      ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
254 
      ierr = VecGetArray(xi,&px);CHKERRQ(ierr);
255 
      ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
256 
      ierr = VecCopy(w,qep->V[i+1]);CHKERRQ(ierr);
257 
      ierr = VecResetArray(w);CHKERRQ(ierr);
258 
      ierr = VecRestoreArray(xi,&px);CHKERRQ(ierr);
259 
      ierr = SlepcVecNormalize(qep->V[i],qep->V[i+1],PETSC_TRUE,PETSC_NULL);CHKERRQ(ierr);
1891 jroman 260 
    }
2295 jroman 261 
    else if (qep->eigi[i]==0.0)   /* real eigenvalue */
1891 jroman 262 
#endif
2295 jroman 263 
    {
264 
      ierr = VecGetArray(xr,&px);CHKERRQ(ierr);
265 
      ierr = VecPlaceArray(w,px+offset);CHKERRQ(ierr);
266 
      ierr = VecCopy(w,qep->V[i]);CHKERRQ(ierr);
267 
      ierr = VecResetArray(w);CHKERRQ(ierr);
268 
      ierr = VecRestoreArray(xr,&px);CHKERRQ(ierr);
269 
      ierr = SlepcVecNormalize(qep->V[i],PETSC_NULL,PETSC_FALSE,PETSC_NULL);CHKERRQ(ierr);
1891 jroman 270 
    }
271 
  }
2305 jroman 272 
  ierr = VecDestroy(&w);CHKERRQ(ierr);
273 
  ierr = VecDestroy(&xr);CHKERRQ(ierr);
274 
  ierr = VecDestroy(&xi);CHKERRQ(ierr);
1887 jroman 275 
  PetscFunctionReturn(0);
276 
}
277 
 
1891 jroman 278 
#undef __FUNCT__  
2324 jroman 279 
#define __FUNCT__ "QEPSolve_Linear"
280 
PetscErrorCode QEPSolve_Linear(QEP qep)
1902 jroman 281 
{
282 
  PetscErrorCode ierr;
283 
  QEP_LINEAR     *ctx = (QEP_LINEAR *)qep->data;
2216 jroman 284 
  PetscBool      flg=PETSC_FALSE;
1902 jroman 285 
 
286 
  PetscFunctionBegin;
287 
  ierr = EPSSolve(ctx->eps);CHKERRQ(ierr);
288 
  ierr = EPSGetConverged(ctx->eps,&qep->nconv);CHKERRQ(ierr);
289 
  ierr = EPSGetIterationNumber(ctx->eps,&qep->its);CHKERRQ(ierr);
290 
  ierr = EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&qep->reason);CHKERRQ(ierr);
291 
  ierr = EPSGetOperationCounters(ctx->eps,&qep->matvecs,PETSC_NULL,&qep->linits);CHKERRQ(ierr);
292 
  qep->matvecs *= 2;  /* convention: count one matvec for each non-trivial block in A */
2216 jroman 293 
  ierr = PetscOptionsGetBool(((PetscObject)qep)->prefix,"-qep_linear_select_simple",&flg,PETSC_NULL);CHKERRQ(ierr);
2118 jroman 294 
  if (flg) {
2295 jroman 295 
    ierr = QEPLinearSelect_Simple(qep,ctx->eps);CHKERRQ(ierr);
2118 jroman 296 
  } else {
2295 jroman 297 
    ierr = QEPLinearSelect_Norm(qep,ctx->eps);CHKERRQ(ierr);
2118 jroman 298 
  }
1902 jroman 299 
  PetscFunctionReturn(0);
300 
}
301 
 
302 
#undef __FUNCT__  
2324 jroman 303 
#define __FUNCT__ "EPSMonitor_Linear"
304 
PetscErrorCode EPSMonitor_Linear(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
1891 jroman 305 
{
2054 eromero 306 
  PetscInt       i;
307 
  QEP            qep = (QEP)ctx;
308 
  PetscErrorCode ierr;
1891 jroman 309 
 
310 
  PetscFunctionBegin;
311 
  nconv = 0;
312 
  for (i=0;i<nest;i++) {
313 
    qep->eigr[i] = eigr[i];
314 
    qep->eigi[i] = eigi[i];
315 
    qep->errest[i] = errest[i];
2038 eromero 316 
    if (0.0 < errest[i] && errest[i] < qep->tol) nconv++;
1891 jroman 317 
  }
2054 eromero 318 
  ierr = STBackTransform(eps->OP,nest,qep->eigr,qep->eigi);CHKERRQ(ierr);
2313 jroman 319 
  ierr = QEPMonitor(qep,its,nconv,qep->eigr,qep->eigi,qep->errest,nest);CHKERRQ(ierr);
1891 jroman 320 
  PetscFunctionReturn(0);
321 
}
322 
 
323 
#undef __FUNCT__  
2324 jroman 324 
#define __FUNCT__ "QEPSetFromOptions_Linear"
325 
PetscErrorCode QEPSetFromOptions_Linear(QEP qep)
1891 jroman 326 
{
327 
  PetscErrorCode ierr;
2216 jroman 328 
  PetscBool      set,val;
1909 jroman 329 
  PetscInt       i;
1891 jroman 330 
  QEP_LINEAR     *ctx = (QEP_LINEAR *)qep->data;
331 
  ST             st;
332 
 
333 
  PetscFunctionBegin;
2324 jroman 334 
  ierr = PetscOptionsBegin(((PetscObject)qep)->comm,((PetscObject)qep)->prefix,"Linear Quadratic Eigenvalue Problem solver Options","QEP");CHKERRQ(ierr);
1909 jroman 335 
  ierr = PetscOptionsInt("-qep_linear_cform","Number of the companion form","QEPLinearSetCompanionForm",ctx->cform,&i,&set);CHKERRQ(ierr);
336 
  if (set) {
337 
    ierr = QEPLinearSetCompanionForm(qep,i);CHKERRQ(ierr);
338 
  }
2216 jroman 339 
  ierr = PetscOptionsBool("-qep_linear_explicitmatrix","Use explicit matrix in linearization","QEPLinearSetExplicitMatrix",ctx->explicitmatrix,&val,&set);CHKERRQ(ierr);
1923 jroman 340 
  if (set) {
341 
    ierr = QEPLinearSetExplicitMatrix(qep,val);CHKERRQ(ierr);
342 
  }
1891 jroman 343 
  if (!ctx->explicitmatrix) {
344 
    /* use as default an ST with shell matrix and Jacobi */
345 
    ierr = EPSGetST(ctx->eps,&st);CHKERRQ(ierr);
1940 jroman 346 
    ierr = STSetMatMode(st,ST_MATMODE_SHELL);CHKERRQ(ierr);
1891 jroman 347 
  }
348 
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
2043 eromero 349 
  ctx->setfromoptionscalled = PETSC_TRUE;
1891 jroman 350 
  PetscFunctionReturn(0);
351 
}
352 
 
1887 jroman 353 
EXTERN_C_BEGIN
354 
#undef __FUNCT__  
2324 jroman 355 
#define __FUNCT__ "QEPLinearSetCompanionForm_Linear"
356 
PetscErrorCode QEPLinearSetCompanionForm_Linear(QEP qep,PetscInt cform)
1909 jroman 357 
{
358 
  QEP_LINEAR *ctx = (QEP_LINEAR *)qep->data;
359 
 
360 
  PetscFunctionBegin;
361 
  if (cform==PETSC_IGNORE) PetscFunctionReturn(0);
362 
  if (cform==PETSC_DECIDE || cform==PETSC_DEFAULT) ctx->cform = 1;
363 
  else {
2214 jroman 364 
    if (cform!=1 && cform!=2) SETERRQ(((PetscObject)qep)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid value of argument 'cform'");
1909 jroman 365 
    ctx->cform = cform;
366 
  }
367 
  PetscFunctionReturn(0);
368 
}
369 
EXTERN_C_END
370 
 
371 
#undef __FUNCT__
372 
#define __FUNCT__ "QEPLinearSetCompanionForm"
373 
/*@
374 
   QEPLinearSetCompanionForm - Choose between the two companion forms available
375 
   for the linearization of the quadratic problem.
376 
 
377 
   Collective on QEP
378 
 
379 
   Input Parameters:
380 
+  qep   - quadratic eigenvalue solver
381 
-  cform - 1 or 2 (first or second companion form)
382 
 
383 
   Options Database Key:
384 
.  -qep_linear_cform <int> - Choose the companion form
385 
 
386 
   Level: advanced
387 
 
388 
.seealso: QEPLinearGetCompanionForm()
389 
@*/
390 
PetscErrorCode QEPLinearSetCompanionForm(QEP qep,PetscInt cform)
391 
{
2221 jroman 392 
  PetscErrorCode ierr;
1909 jroman 393 
 
394 
  PetscFunctionBegin;
2213 jroman 395 
  PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
2326 jroman 396 
  PetscValidLogicalCollectiveInt(qep,cform,2);
2221 jroman 397 
  ierr = PetscTryMethod(qep,"QEPLinearSetCompanionForm_C",(QEP,PetscInt),(qep,cform));CHKERRQ(ierr);
1909 jroman 398 
  PetscFunctionReturn(0);
399 
}
400 
 
401 
EXTERN_C_BEGIN
402 
#undef __FUNCT__  
2324 jroman 403 
#define __FUNCT__ "QEPLinearGetCompanionForm_Linear"
404 
PetscErrorCode QEPLinearGetCompanionForm_Linear(QEP qep,PetscInt *cform)
1909 jroman 405 
{
406 
  QEP_LINEAR *ctx = (QEP_LINEAR *)qep->data;
407 
 
408 
  PetscFunctionBegin;
409 
  *cform = ctx->cform;
410 
  PetscFunctionReturn(0);
411 
}
412 
EXTERN_C_END
413 
 
414 
#undef __FUNCT__
415 
#define __FUNCT__ "QEPLinearGetCompanionForm"
2077 eromero 416 
/*@
1909 jroman 417 
   QEPLinearGetCompanionForm - Returns the number of the companion form that
418 
   will be used for the linearization of the quadratic problem.
419 
 
420 
   Not collective
421 
 
422 
   Input Parameter:
423 
.  qep  - quadratic eigenvalue solver
424 
 
425 
   Output Parameter:
426 
.  cform - the companion form number (1 or 2)
427 
 
428 
   Level: advanced
429 
 
430 
.seealso: QEPLinearSetCompanionForm()
431 
@*/
432 
PetscErrorCode QEPLinearGetCompanionForm(QEP qep,PetscInt *cform)
433 
{
2221 jroman 434 
  PetscErrorCode ierr;
1909 jroman 435 
 
436 
  PetscFunctionBegin;
2213 jroman 437 
  PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
2326 jroman 438 
  PetscValidIntPointer(cform,2);
2221 jroman 439 
  ierr = PetscTryMethod(qep,"QEPLinearGetCompanionForm_C",(QEP,PetscInt*),(qep,cform));CHKERRQ(ierr);
1909 jroman 440 
  PetscFunctionReturn(0);
441 
}
442 
 
443 
EXTERN_C_BEGIN
444 
#undef __FUNCT__  
2324 jroman 445 
#define __FUNCT__ "QEPLinearSetExplicitMatrix_Linear"
446 
PetscErrorCode QEPLinearSetExplicitMatrix_Linear(QEP qep,PetscBool explicitmatrix)
1891 jroman 447 
{
448 
  QEP_LINEAR *ctx = (QEP_LINEAR *)qep->data;
449 
 
450 
  PetscFunctionBegin;
451 
  ctx->explicitmatrix = explicitmatrix;
452 
  PetscFunctionReturn(0);
453 
}
454 
EXTERN_C_END
455 
 
456 
#undef __FUNCT__
457 
#define __FUNCT__ "QEPLinearSetExplicitMatrix"
458 
/*@
459 
   QEPLinearSetExplicitMatrix - Indicate if the matrices A and B for the linearization
460 
   of the quadratic problem must be built explicitly.
461 
 
462 
   Collective on QEP
463 
 
464 
   Input Parameters:
465 
+  qep      - quadratic eigenvalue solver
466 
-  explicit - boolean flag indicating if the matrices are built explicitly
467 
 
468 
   Options Database Key:
469 
.  -qep_linear_explicitmatrix <boolean> - Indicates the boolean flag
470 
 
471 
   Level: advanced
472 
 
473 
.seealso: QEPLinearGetExplicitMatrix()
474 
@*/
2216 jroman 475 
PetscErrorCode QEPLinearSetExplicitMatrix(QEP qep,PetscBool explicitmatrix)
1891 jroman 476 
{
2221 jroman 477 
  PetscErrorCode ierr;
1891 jroman 478 
 
479 
  PetscFunctionBegin;
2213 jroman 480 
  PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
2326 jroman 481 
  PetscValidLogicalCollectiveBool(qep,explicitmatrix,2);
2221 jroman 482 
  ierr = PetscTryMethod(qep,"QEPLinearSetExplicitMatrix_C",(QEP,PetscBool),(qep,explicitmatrix));CHKERRQ(ierr);
1891 jroman 483 
  PetscFunctionReturn(0);
484 
}
485 
 
486 
EXTERN_C_BEGIN
487 
#undef __FUNCT__  
2324 jroman 488 
#define __FUNCT__ "QEPLinearGetExplicitMatrix_Linear"
489 
PetscErrorCode QEPLinearGetExplicitMatrix_Linear(QEP qep,PetscBool *explicitmatrix)
1891 jroman 490 
{
491 
  QEP_LINEAR *ctx = (QEP_LINEAR *)qep->data;
492 
 
493 
  PetscFunctionBegin;
494 
  *explicitmatrix = ctx->explicitmatrix;
495 
  PetscFunctionReturn(0);
496 
}
497 
EXTERN_C_END
498 
 
499 
#undef __FUNCT__
500 
#define __FUNCT__ "QEPLinearGetExplicitMatrix"
2077 eromero 501 
/*@
1891 jroman 502 
   QEPLinearGetExplicitMatrix - Returns the flag indicating if the matrices A and B
1909 jroman 503 
   for the linearization of the quadratic problem are built explicitly.
1891 jroman 504 
 
505 
   Not collective
506 
 
507 
   Input Parameter:
508 
.  qep  - quadratic eigenvalue solver
509 
 
510 
   Output Parameter:
1904 jroman 511 
.  explicitmatrix - the mode flag
1891 jroman 512 
 
513 
   Level: advanced
514 
 
515 
.seealso: QEPLinearSetExplicitMatrix()
516 
@*/
2216 jroman 517 
PetscErrorCode QEPLinearGetExplicitMatrix(QEP qep,PetscBool *explicitmatrix)
1891 jroman 518 
{
2221 jroman 519 
  PetscErrorCode ierr;
1891 jroman 520 
 
521 
  PetscFunctionBegin;
2213 jroman 522 
  PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
2326 jroman 523 
  PetscValidPointer(explicitmatrix,2);
2221 jroman 524 
  ierr = PetscTryMethod(qep,"QEPLinearGetExplicitMatrix_C",(QEP,PetscBool*),(qep,explicitmatrix));CHKERRQ(ierr);
1891 jroman 525 
  PetscFunctionReturn(0);
526 
}
527 
 
528 
EXTERN_C_BEGIN
529 
#undef __FUNCT__  
2324 jroman 530 
#define __FUNCT__ "QEPLinearSetEPS_Linear"
531 
PetscErrorCode QEPLinearSetEPS_Linear(QEP qep,EPS eps)
1891 jroman 532 
{
2317 jroman 533 
  PetscErrorCode ierr;
534 
  QEP_LINEAR     *ctx = (QEP_LINEAR *)qep->data;
1891 jroman 535 
 
536 
  PetscFunctionBegin;
537 
  ierr = PetscObjectReference((PetscObject)eps);CHKERRQ(ierr);
2308 jroman 538 
  ierr = EPSDestroy(&ctx->eps);CHKERRQ(ierr);  
1891 jroman 539 
  ctx->eps = eps;
540 
  qep->setupcalled = 0;
541 
  PetscFunctionReturn(0);
542 
}
543 
EXTERN_C_END
544 
 
545 
#undef __FUNCT__  
546 
#define __FUNCT__ "QEPLinearSetEPS"
547 
/*@
548 
   QEPLinearSetEPS - Associate an eigensolver object (EPS) to the
549 
   quadratic eigenvalue solver.
550 
 
551 
   Collective on QEP
552 
 
553 
   Input Parameters:
554 
+  qep - quadratic eigenvalue solver
555 
-  eps - the eigensolver object
556 
 
557 
   Level: advanced
558 
 
559 
.seealso: QEPLinearGetEPS()
560 
@*/
561 
PetscErrorCode QEPLinearSetEPS(QEP qep,EPS eps)
562 
{
2221 jroman 563 
  PetscErrorCode ierr;
1891 jroman 564 
 
565 
  PetscFunctionBegin;
2213 jroman 566 
  PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
2221 jroman 567 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
2326 jroman 568 
  PetscCheckSameComm(qep,1,eps,2);
2221 jroman 569 
  ierr = PetscTryMethod(qep,"QEPLinearSetEPS_C",(QEP,EPS),(qep,eps));CHKERRQ(ierr);
1891 jroman 570 
  PetscFunctionReturn(0);
571 
}
572 
 
573 
EXTERN_C_BEGIN
574 
#undef __FUNCT__  
2324 jroman 575 
#define __FUNCT__ "QEPLinearGetEPS_Linear"
576 
PetscErrorCode QEPLinearGetEPS_Linear(QEP qep,EPS *eps)
1891 jroman 577 
{
578 
  QEP_LINEAR *ctx = (QEP_LINEAR *)qep->data;
579 
 
580 
  PetscFunctionBegin;
581 
  *eps = ctx->eps;
582 
  PetscFunctionReturn(0);
583 
}
584 
EXTERN_C_END
585 
 
586 
#undef __FUNCT__  
587 
#define __FUNCT__ "QEPLinearGetEPS"
2077 eromero 588 
/*@
1891 jroman 589 
   QEPLinearGetEPS - Retrieve the eigensolver object (EPS) associated
590 
   to the quadratic eigenvalue solver.
591 
 
592 
   Not Collective
593 
 
594 
   Input Parameter:
595 
.  qep - quadratic eigenvalue solver
596 
 
597 
   Output Parameter:
598 
.  eps - the eigensolver object
599 
 
600 
   Level: advanced
601 
 
602 
.seealso: QEPLinearSetEPS()
603 
@*/
604 
PetscErrorCode QEPLinearGetEPS(QEP qep,EPS *eps)
605 
{
2221 jroman 606 
  PetscErrorCode ierr;
1891 jroman 607 
 
608 
  PetscFunctionBegin;
2213 jroman 609 
  PetscValidHeaderSpecific(qep,QEP_CLASSID,1);
2326 jroman 610 
  PetscValidPointer(eps,2);
2221 jroman 611 
  ierr = PetscTryMethod(qep,"QEPLinearGetEPS_C",(QEP,EPS*),(qep,eps));CHKERRQ(ierr);
1891 jroman 612 
  PetscFunctionReturn(0);
613 
}
614 
 
615 
#undef __FUNCT__  
2324 jroman 616 
#define __FUNCT__ "QEPView_Linear"
617 
PetscErrorCode QEPView_Linear(QEP qep,PetscViewer viewer)
1891 jroman 618 
{
2317 jroman 619 
  PetscErrorCode ierr;
620 
  QEP_LINEAR     *ctx = (QEP_LINEAR *)qep->data;
1891 jroman 621 
 
622 
  PetscFunctionBegin;
623 
  if (ctx->explicitmatrix) {
624 
    ierr = PetscViewerASCIIPrintf(viewer,"linearized matrices: explicit\n");CHKERRQ(ierr);
625 
  } else {
626 
    ierr = PetscViewerASCIIPrintf(viewer,"linearized matrices: implicit\n");CHKERRQ(ierr);
627 
  }
1909 jroman 628 
  ierr = PetscViewerASCIIPrintf(viewer,"companion form: %d\n",ctx->cform);CHKERRQ(ierr);
1891 jroman 629 
  ierr = EPSView(ctx->eps,viewer);CHKERRQ(ierr);
630 
  PetscFunctionReturn(0);
631 
}
632 
 
633 
#undef __FUNCT__  
2324 jroman 634 
#define __FUNCT__ "QEPDestroy_Linear"
635 
PetscErrorCode QEPDestroy_Linear(QEP qep)
1891 jroman 636 
{
2317 jroman 637 
  PetscErrorCode ierr;
638 
  QEP_LINEAR     *ctx = (QEP_LINEAR *)qep->data;
1891 jroman 639 
 
640 
  PetscFunctionBegin;
2308 jroman 641 
  ierr = EPSDestroy(&ctx->eps);CHKERRQ(ierr);
2305 jroman 642 
  ierr = MatDestroy(&ctx->A);CHKERRQ(ierr);
643 
  ierr = MatDestroy(&ctx->B);CHKERRQ(ierr);
644 
  ierr = VecDestroy(&ctx->x1);CHKERRQ(ierr);
645 
  ierr = VecDestroy(&ctx->x2);CHKERRQ(ierr);
646 
  ierr = VecDestroy(&ctx->y1);CHKERRQ(ierr);
647 
  ierr = VecDestroy(&ctx->y2);CHKERRQ(ierr);
1925 jroman 648 
 
649 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearSetCompanionForm_C","",PETSC_NULL);CHKERRQ(ierr);
650 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearGetCompanionForm_C","",PETSC_NULL);CHKERRQ(ierr);
651 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearSetEPS_C","",PETSC_NULL);CHKERRQ(ierr);
652 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearGetEPS_C","",PETSC_NULL);CHKERRQ(ierr);
653 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearSetExplicitMatrix_C","",PETSC_NULL);CHKERRQ(ierr);
654 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearGetExplicitMatrix_C","",PETSC_NULL);CHKERRQ(ierr);
655 
 
1903 jroman 656 
  ierr = QEPDestroy_Default(qep);CHKERRQ(ierr);
1891 jroman 657 
  PetscFunctionReturn(0);
658 
}
659 
 
660 
EXTERN_C_BEGIN
661 
#undef __FUNCT__  
2324 jroman 662 
#define __FUNCT__ "QEPCreate_Linear"
663 
PetscErrorCode QEPCreate_Linear(QEP qep)
1887 jroman 664 
{
665 
  PetscErrorCode ierr;
1891 jroman 666 
  QEP_LINEAR     *ctx;
1887 jroman 667 
 
668 
  PetscFunctionBegin;
1891 jroman 669 
  ierr = PetscNew(QEP_LINEAR,&ctx);CHKERRQ(ierr);
1887 jroman 670 
  PetscLogObjectMemory(qep,sizeof(QEP_LINEAR));
1891 jroman 671 
  qep->data                      = (void *)ctx;
2324 jroman 672 
  qep->ops->solve                = QEPSolve_Linear;
673 
  qep->ops->setup                = QEPSetUp_Linear;
674 
  qep->ops->setfromoptions       = QEPSetFromOptions_Linear;
675 
  qep->ops->destroy              = QEPDestroy_Linear;
676 
  qep->ops->view                 = QEPView_Linear;
677 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearSetCompanionForm_C","QEPLinearSetCompanionForm_Linear",QEPLinearSetCompanionForm_Linear);CHKERRQ(ierr);
678 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearGetCompanionForm_C","QEPLinearGetCompanionForm_Linear",QEPLinearGetCompanionForm_Linear);CHKERRQ(ierr);
679 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearSetEPS_C","QEPLinearSetEPS_Linear",QEPLinearSetEPS_Linear);CHKERRQ(ierr);
680 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearGetEPS_C","QEPLinearGetEPS_Linear",QEPLinearGetEPS_Linear);CHKERRQ(ierr);
681 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearSetExplicitMatrix_C","QEPLinearSetExplicitMatrix_Linear",QEPLinearSetExplicitMatrix_Linear);CHKERRQ(ierr);
682 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)qep,"QEPLinearGetExplicitMatrix_C","QEPLinearGetExplicitMatrix_Linear",QEPLinearGetExplicitMatrix_Linear);CHKERRQ(ierr);
1891 jroman 683 
 
684 
  ierr = EPSCreate(((PetscObject)qep)->comm,&ctx->eps);CHKERRQ(ierr);
685 
  ierr = EPSSetOptionsPrefix(ctx->eps,((PetscObject)qep)->prefix);CHKERRQ(ierr);
686 
  ierr = EPSAppendOptionsPrefix(ctx->eps,"qep_");CHKERRQ(ierr);
687 
  ierr = PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)qep,1);CHKERRQ(ierr);  
688 
  PetscLogObjectParent(qep,ctx->eps);
689 
  ierr = EPSSetIP(ctx->eps,qep->ip);CHKERRQ(ierr);
2324 jroman 690 
  ierr = EPSMonitorSet(ctx->eps,EPSMonitor_Linear,qep,PETSC_NULL);CHKERRQ(ierr);
1891 jroman 691 
  ctx->explicitmatrix = PETSC_FALSE;
1909 jroman 692 
  ctx->cform = 1;
1891 jroman 693 
  ctx->A = PETSC_NULL;
694 
  ctx->B = PETSC_NULL;
695 
  ctx->x1 = PETSC_NULL;
696 
  ctx->x2 = PETSC_NULL;
697 
  ctx->y1 = PETSC_NULL;
698 
  ctx->y2 = PETSC_NULL;
2043 eromero 699 
  ctx->setfromoptionscalled = PETSC_FALSE;
1887 jroman 700 
  PetscFunctionReturn(0);
701 
}
702 
EXTERN_C_END
703