Actual source code: ex57.c

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

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: static char help[] = "Another eigenvalue problem with Hamiltonian structure.\n\n"
 12:   "Position and velocity control of a string of high-speed vehicles.\n"
 13:   "The command line options are:\n"
 14:   "  -m <m>, where <n> = number of vehicles.\n\n";

 16: #include <slepceps.h>

 18: /*
 19:    This example computes eigenvalues of a matrix

 21:         H = [ A    N
 22:               K  -A^* ],

 24:    where A, N and K are n-by-n matrices, with n=2*m-1 and

 26:         N = diag( 1, 0, 1, 0,..., 0, 1)
 27:         K = diag( 0,10, 0,10,...,10, 0)
 28:         A = tridiag(b,d,c)
 29:         d = [-1, 0,-1, 0,...,-1, 0,-1]
 30:         b = [ 1, 0, 1, 0,..., 1, 0]
 31:         c = [ 0,-1, 0,-1,..., 0,-1]

 33:    References:

 35:    [1] W.R. Ferng, Wen-Wei Lin, Chern-Shuh Wang, The shift-inverted J-Lanczos algorithm
 36:        for the numerical solutions of large sparse algebraic Riccati equations,
 37:        Computers & Mathematics with Applications, Volume 33, Issue 10, 1997,
 38:        https://doi.org/10.1016/S0898-1221(97)00074-6.

 40:    */

 42: static PetscErrorCode eigenCompare (PetscScalar ar,PetscScalar ai,PetscScalar br,PetscScalar bi,PetscInt *res,void *ctx)
 43: {
 44:   PetscReal abs1, abs2, tol=1e-12, r_ar, r_ai, r_br, r_bi;

 46:   PetscFunctionBegin;
 47:   #if defined(PETSC_USE_COMPLEX)
 48:     r_ar = PetscRealPart(ar);
 49:     r_ai = PetscImaginaryPart(ar);
 50:     r_br = PetscRealPart(br);
 51:     r_bi = PetscImaginaryPart(br);
 52:   #else
 53:     r_ar = ar;
 54:     r_ai = ai;
 55:     r_br = br;
 56:     r_bi = bi;
 57:   #endif
 58:   if (PetscAbs(PetscAbs(r_ar)-PetscAbs(r_br)) < tol && PetscAbs(PetscAbs(r_ai)-PetscAbs(r_bi)) < tol) {
 59:     /* Negative real part first*/
 60:     if (r_ar<0 && r_br>0)
 61:       *res = -1;
 62:     else if (r_ar>0 && r_br<0)
 63:       *res = 1;
 64:     /* Positive imaginary part first*/
 65:     else if (r_ai>0 && r_bi<0)
 66:       *res = -1;
 67:     else if (r_ai<0 && r_bi>0)
 68:       *res = 1;
 69:   }
 70:   else {
 71:     abs1 = SlepcAbs(r_ar,r_ai);
 72:     abs2 = SlepcAbs(r_br,r_bi);
 73:     if (abs1 > abs2)
 74:       *res = -1;
 75:     else if (abs1 < abs2)
 76:       *res = 1;
 77:     else
 78:       *res = 0;
 79:   }
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

 83: int main(int argc,char **argv)
 84: {
 85:   Mat            H,A,N,K;    /* problem matrices */
 86:   EPS            eps;        /* eigenproblem solver context */
 87:   PetscInt       m=12,n,Istart,Iend,i;
 88:   PetscBool      terse, sort_hamilt;

 90:   PetscFunctionBeginUser;
 91:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));

 93:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
 94:   n = 2*m-1;
 95:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHamiltonian eigenproblem, n=%" PetscInt_FMT
 96:     " (m=%" PetscInt_FMT " vehicles)\n\n",n,m));

 98:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:                Compute the problem matrices A, N and K
100:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

102:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
103:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
104:   PetscCall(MatSetFromOptions(A));

106:   PetscCall(MatCreate(PETSC_COMM_WORLD,&N));
107:   PetscCall(MatSetSizes(N,PETSC_DECIDE,PETSC_DECIDE,n,n));
108:   PetscCall(MatSetFromOptions(N));

110:   PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
111:   PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
112:   PetscCall(MatSetFromOptions(K));

114:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
115:   for (i=Istart;i<Iend;i++) {
116:     if (i%2==0)
117:       PetscCall(MatSetValue(A,i,i,-1.0,INSERT_VALUES));
118:     else {
119:       if (i>0) PetscCall(MatSetValue(A,i,i-1,1.0,INSERT_VALUES));
120:       if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
121:     }
122:   }

124:   PetscCall(MatGetOwnershipRange(N,&Istart,&Iend));
125:   for (i=Istart+Istart%2;i<Iend;i+=2) {
126:     PetscCall(MatSetValue(N,i,i,1.0,INSERT_VALUES));
127:   }

129:   PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
130:   for (i=Istart+(Istart+1)%2;i<Iend;i+=2) {
131:     PetscCall(MatSetValue(K,i,i,10.0,INSERT_VALUES));
132:   }

134:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
135:   PetscCall(MatAssemblyBegin(N,MAT_FINAL_ASSEMBLY));
136:   PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
137:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
138:   PetscCall(MatAssemblyEnd(N,MAT_FINAL_ASSEMBLY));
139:   PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));

141:   PetscCall(MatCreateHamiltonian(A,N,K,&H));

143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:                 Create the eigensolver and set various options
145:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

147:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
148:   PetscCall(EPSSetOperators(eps,H,NULL));
149:   PetscCall(EPSSetProblemType(eps,EPS_HAMILT));

151:   PetscCall(PetscOptionsHasName(NULL,NULL,"-sort_hamilt",&sort_hamilt));
152:   if (sort_hamilt) {
153:     /* Adjust ordering of non-hermitian solver so that it is the same as in as in EPS_HAMILT solver */
154:     PetscCall(EPSSetEigenvalueComparison(eps,eigenCompare,NULL));
155:     PetscCall(EPSSetWhichEigenpairs(eps,EPS_WHICH_USER));
156:   }
157:   PetscCall(EPSSetFromOptions(eps));

159:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160:                  Solve the eigensystem and display solution
161:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

163:   PetscCall(EPSSolve(eps));

165:   /* show detailed info unless -terse option is given by user */
166:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
167:   if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
168:   else {
169:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
170:     PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
171:     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
172:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
173:   }

175:   PetscCall(EPSDestroy(&eps));
176:   PetscCall(MatDestroy(&A));
177:   PetscCall(MatDestroy(&N));
178:   PetscCall(MatDestroy(&K));
179:   PetscCall(MatDestroy(&H));
180:   PetscCall(SlepcFinalize());
181:   return 0;
182: }

184: /*TEST

186:    testset:
187:       args: -eps_nev 8 -eps_ncv 28 -terse
188:       nsize: {{1 2}}
189:       output_file: output/ex57_1.out
190:       test:
191:          requires: double !complex
192:          suffix: 1
193:       test:
194:          requires: double !complex
195:          args: -eps_non_hermitian -sort_hamilt
196:          suffix: 1_nhep
197:       test:
198:          requires: double complex
199:          TODO: no support for complex scalars yet
200:          suffix: 1_complex
201:       test:
202:          requires: double complex
203:          args: -eps_non_hermitian -sort_hamilt
204:          suffix: 1_complex_nhep

206:    testset:
207:       args: -eps_nev 4 -eps_smallest_magnitude -eps_ncv 32 -terse
208:       output_file: output/ex57_2.out
209:       test:
210:          requires: double !complex
211:          suffix: 2

213: TEST*/