Subversion Repositories slepc-dev

Rev

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