Actual source code: test2.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2013, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.

  8:    SLEPc is free software: you can redistribute it and/or modify it under  the
  9:    terms of version 3 of the GNU Lesser General Public License as published by
 10:    the Free Software Foundation.

 12:    SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
 13:    WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
 14:    FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
 15:    more details.

 17:    You  should have received a copy of the GNU Lesser General  Public  License
 18:    along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
 19:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 21:    Example based on spring problem in NLEVP collection [1]. See the parameters
 22:    meaning at Example 2 in [2].

 24:    [1] T. Betcke, N. J. Higham, V. Mehrmann, C. Schroder, and F. Tisseur,
 25:        NLEVP: A Collection of Nonlinear Eigenvalue Problems, MIMS EPrint
 26:        2010.98, November 2010.
 27:    [2] F. Tisseur, Backward error and condition of polynomial eigenvalue
 28:        problems, Linear Algebra and its Applications, 309 (2000), pp. 339--361,
 29:        April 2000.
 30: */

 32: static char help[] = "Test the solution of a QEP from a finite element model of "
 33:   "damped mass-spring system (problem from NLEVP collection).\n\n"
 34:   "The command line options are:\n"
 35:   "  -n <n>, where <n> = number of grid subdivisions.\n"
 36:   "  -tau <tau>, where <tau> = tau parameter.\n"
 37:   "  -kappa <kappa>, where <kappa> = kappa paramter.\n\n";

 39: #include <slepcqep.h>

 43: int main(int argc,char **argv)
 44: {
 45:   Mat            M,C,K;           /* problem matrices */
 46:   QEP            qep;             /* quadratic eigenproblem solver context */
 47:   QEPType        type;
 49:   PetscInt       n=30,Istart,Iend,i,maxit,nev;
 50:   PetscScalar    mu=1.0,tau=10.0,kappa=5.0;

 52:   SlepcInitialize(&argc,&argv,(char*)0,help);

 54:   PetscOptionsGetInt(NULL,"-n",&n,NULL);
 55:   PetscOptionsGetScalar(NULL,"-mu",&mu,NULL);
 56:   PetscOptionsGetScalar(NULL,"-tau",&tau,NULL);
 57:   PetscOptionsGetScalar(NULL,"-kappa",&kappa,NULL);

 59:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 60:      Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
 61:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 63:   /* K is a tridiagonal */
 64:   MatCreate(PETSC_COMM_WORLD,&K);
 65:   MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
 66:   MatSetFromOptions(K);
 67:   MatSetUp(K);

 69:   MatGetOwnershipRange(K,&Istart,&Iend);
 70:   for (i=Istart; i<Iend; i++) {
 71:     if (i>0) {
 72:       MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
 73:     }
 74:     MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
 75:     if (i<n-1) {
 76:       MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
 77:     }
 78:   }

 80:   MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
 81:   MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);

 83:   /* C is a tridiagonal */
 84:   MatCreate(PETSC_COMM_WORLD,&C);
 85:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
 86:   MatSetFromOptions(C);
 87:   MatSetUp(C);

 89:   MatGetOwnershipRange(C,&Istart,&Iend);
 90:   for (i=Istart; i<Iend; i++) {
 91:     if (i>0) {
 92:       MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
 93:     }
 94:     MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
 95:     if (i<n-1) {
 96:       MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
 97:     }
 98:   }

100:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
101:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

103:   /* M is a diagonal matrix */
104:   MatCreate(PETSC_COMM_WORLD,&M);
105:   MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
106:   MatSetFromOptions(M);
107:   MatSetUp(M);
108:   MatGetOwnershipRange(M,&Istart,&Iend);
109:   for (i=Istart; i<Iend; i++) {
110:     MatSetValue(M,i,i,mu,INSERT_VALUES);
111:   }
112:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
113:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:                 Create the eigensolver and set various options
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

119:   /*
120:      Create eigensolver context
121:   */
122:   QEPCreate(PETSC_COMM_WORLD,&qep);

124:   /*
125:      Set matrices, the problem type and other settings
126:   */
127:   QEPSetOperators(qep,M,C,K);
128:   QEPSetProblemType(qep,QEP_GENERAL);
129:   QEPSetTolerances(qep,PETSC_SMALL,PETSC_DEFAULT);
130:   QEPSetFromOptions(qep);

132:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133:                       Solve the eigensystem
134:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

136:   QEPSolve(qep);

138:   /*
139:      Optional: Get some information from the solver and display it
140:   */
141:   QEPGetType(qep,&type);
142:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
143:   QEPGetDimensions(qep,&nev,NULL,NULL);
144:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
145:   QEPGetTolerances(qep,NULL,&maxit);
146:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: maxit=%D\n",maxit);

148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149:                     Display solution and clean up
150:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

152:   QEPPrintSolution(qep,NULL);

154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:      Free work space
156:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

158:   QEPDestroy(&qep);
159:   MatDestroy(&M);
160:   MatDestroy(&C);
161:   MatDestroy(&K);
162:   SlepcFinalize();
163:   return 0;
164: }