Actual source code: ex24.c

slepc-3.7.4 2017-05-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2016, 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:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 20: */

 22: static char help[] = "Spectrum folding for a standard symmetric eigenproblem.\n\n"
 23:   "The problem matrix is the 2-D Laplacian.\n\n"
 24:   "The command line options are:\n"
 25:   "  -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
 26:   "  -m <m>, where <m> = number of grid subdivisions in y dimension.\n";

 28: #include <slepceps.h>

 30: /*
 31:    User context for spectrum folding
 32: */
 33: typedef struct {
 34:   Mat       A;
 35:   Vec       w;
 36:   PetscReal target;
 37: } CTX_FOLD;

 39: /*
 40:    Auxiliary routines
 41: */
 42: PetscErrorCode MatMult_Fold(Mat,Vec,Vec);
 43: PetscErrorCode RayleighQuotient(Mat,Vec,PetscScalar*);
 44: PetscErrorCode ComputeResidualNorm(Mat,PetscScalar,Vec,PetscReal*);

 48: int main(int argc,char **argv)
 49: {
 50:   Mat            A,M,P;       /* problem matrix, shell matrix and preconditioner */
 51:   Vec            x;           /* eigenvector */
 52:   EPS            eps;         /* eigenproblem solver context */
 53:   ST             st;          /* spectral transformation */
 54:   KSP            ksp;
 55:   PC             pc;
 56:   EPSType        type;
 57:   CTX_FOLD       *ctx;
 58:   PetscInt       nconv,N,n=10,m,Istart,Iend,II,i,j;
 59:   PetscReal      *error,*evals,target=2.1,tol;
 60:   PetscScalar    lambda;
 61:   PetscBool      flag,terse,errok;

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

 66:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 67:   PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag);
 68:   if (!flag) m=n;
 69:   PetscOptionsGetReal(NULL,NULL,"-target",&target,NULL);
 70:   N = n*m;
 71:   PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum Folding, N=%D (%Dx%D grid) target=%f\n\n",N,n,m,(double)target);

 73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:      Compute the 5-point stencil Laplacian
 75:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 77:   MatCreate(PETSC_COMM_WORLD,&A);
 78:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
 79:   MatSetFromOptions(A);
 80:   MatSetUp(A);

 82:   MatGetOwnershipRange(A,&Istart,&Iend);
 83:   for (II=Istart;II<Iend;II++) {
 84:     i = II/n; j = II-i*n;
 85:     if (i>0) { MatSetValue(A,II,II-n,-1.0,INSERT_VALUES); }
 86:     if (i<m-1) { MatSetValue(A,II,II+n,-1.0,INSERT_VALUES); }
 87:     if (j>0) { MatSetValue(A,II,II-1,-1.0,INSERT_VALUES); }
 88:     if (j<n-1) { MatSetValue(A,II,II+1,-1.0,INSERT_VALUES); }
 89:     MatSetValue(A,II,II,4.0,INSERT_VALUES);
 90:   }

 92:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 93:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 94:   MatCreateVecs(A,&x,NULL);

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:                 Create shell matrix to perform spectrum folding
 98:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 99:   PetscNew(&ctx);
100:   ctx->A = A;
101:   ctx->target = target;
102:   VecDuplicate(x,&ctx->w);

104:   MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,ctx,&M);
105:   MatSetFromOptions(M);
106:   MatShellSetOperation(M,MATOP_MULT,(void(*)())MatMult_Fold);

108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109:                 Create the eigensolver and set various options
110:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

112:   EPSCreate(PETSC_COMM_WORLD,&eps);
113:   EPSSetOperators(eps,M,NULL);
114:   EPSSetProblemType(eps,EPS_HEP);
115:   EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
116:   EPSSetFromOptions(eps);

118:   PetscObjectTypeCompareAny((PetscObject)eps,&flag,EPSBLOPEX,EPSRQCG,"");
119:   if (flag) {
120:     /*
121:        Build preconditioner specific for this application (diagonal of A^2)
122:     */
123:     MatGetRowSum(A,x);
124:     VecScale(x,-1.0);
125:     VecShift(x,20.0);
126:     MatCreate(PETSC_COMM_WORLD,&P);
127:     MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,N);
128:     MatSetFromOptions(P);
129:     MatSetUp(P);
130:     MatDiagonalSet(P,x,INSERT_VALUES);
131:     /*
132:        Set diagonal preconditioner
133:     */
134:     EPSGetST(eps,&st);
135:     STSetType(st,STPRECOND);
136:     STPrecondSetMatForPC(st,P);
137:     MatDestroy(&P);
138:     STGetKSP(st,&ksp);
139:     KSPGetPC(ksp,&pc);
140:     PCSetType(pc,PCJACOBI);
141:   }

143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:                       Solve the eigensystem
145:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

147:   EPSSolve(eps);
148:   EPSGetType(eps,&type);
149:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
150:   EPSGetTolerances(eps,&tol,NULL);

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:                     Display solution and clean up
154:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

156:   EPSGetConverged(eps,&nconv);
157:   PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);
158:   if (nconv>0) {
159:     PetscMalloc2(nconv,&evals,nconv,&error);
160:     for (i=0;i<nconv;i++) {
161:       /*  Get i-th eigenvector, compute eigenvalue approximation from
162:           Rayleigh quotient and compute residual norm */
163:       EPSGetEigenpair(eps,i,NULL,NULL,x,NULL);
164:       RayleighQuotient(A,x,&lambda);
165:       ComputeResidualNorm(A,lambda,x,&error[i]);
166: #if defined(PETSC_USE_COMPLEX)
167:       evals[i] = PetscRealPart(lambda);
168: #else
169:       evals[i] = lambda;
170: #endif
171:     }
172:     PetscOptionsHasName(NULL,NULL,"-terse",&terse);
173:     if (!terse) {
174:       PetscPrintf(PETSC_COMM_WORLD,
175:            "           k              ||Ax-kx||\n"
176:            "   ----------------- ------------------\n");
177:       for (i=0;i<nconv;i++) {
178:         PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12.2g\n",(double)evals[i],(double)error[i]);
179:       }
180:     } else {
181:       errok = PETSC_TRUE;
182:       for (i=0;i<nconv;i++) errok = (errok && error[i]<5.0*tol)? PETSC_TRUE: PETSC_FALSE;
183:       if (!errok) {
184:         PetscPrintf(PETSC_COMM_WORLD," Problem: some of the first %D relative errors are higher than the tolerance\n\n",nconv);
185:       } else {
186:         PetscPrintf(PETSC_COMM_WORLD," nconv=%D eigenvalues computed up to the required tolerance:",nconv);
187:         for (i=0;i<nconv;i++) {
188:           PetscPrintf(PETSC_COMM_WORLD," %.5f",(double)evals[i]);
189:         }
190:       }
191:     }
192:     PetscPrintf(PETSC_COMM_WORLD,"\n");
193:     PetscFree2(evals,error);
194:   }

196:   EPSDestroy(&eps);
197:   MatDestroy(&A);
198:   MatDestroy(&M);
199:   VecDestroy(&ctx->w);
200:   VecDestroy(&x);
201:   PetscFree(ctx);
202:   SlepcFinalize();
203:   return ierr;
204: }

208: /*
209:     Matrix-vector product subroutine for the spectrum folding.
210:        y <-- (A-t*I)^2*x
211:  */
212: PetscErrorCode MatMult_Fold(Mat M,Vec x,Vec y)
213: {
214:   CTX_FOLD       *ctx;
215:   PetscScalar    sigma;

219:   MatShellGetContext(M,&ctx);
220:   sigma = -ctx->target;
221:   MatMult(ctx->A,x,ctx->w);
222:   VecAXPY(ctx->w,sigma,x);
223:   MatMult(ctx->A,ctx->w,y);
224:   VecAXPY(y,sigma,ctx->w);
225:   return(0);
226: }

230: /*
231:     Computes the Rayleigh quotient of a vector x
232:        r <-- x^T*A*x       (assumes x has unit norm)
233:  */
234: PetscErrorCode RayleighQuotient(Mat A,Vec x,PetscScalar *r)
235: {
236:   Vec            Ax;

240:   VecDuplicate(x,&Ax);
241:   MatMult(A,x,Ax);
242:   VecDot(Ax,x,r);
243:   VecDestroy(&Ax);
244:   return(0);
245: }

249: /*
250:     Computes the residual norm of an approximate eigenvector x, |A*x-lambda*x|
251:  */
252: PetscErrorCode ComputeResidualNorm(Mat A,PetscScalar lambda,Vec x,PetscReal *r)
253: {
254:   Vec            Ax;

258:   VecDuplicate(x,&Ax);
259:   MatMult(A,x,Ax);
260:   VecAXPY(Ax,-lambda,x);
261:   VecNorm(Ax,NORM_2,r);
262:   VecDestroy(&Ax);
263:   return(0);
264: }