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