Actual source code: test1.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7:
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[] = "Test ST with shell matrices as problem matrices.\n\n";
24: #include <slepceps.h>
26: static PetscErrorCode MatShift_Shell(Mat S,PetscScalar shift);
27: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag);
28: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y);
29: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y);
33: int main(int argc,char **argv)
34: {
35: Mat A, S; /* problem matrix */
36: EPS eps; /* eigenproblem solver context */
37: const EPSType type;
38: PetscScalar value[3];
39: PetscInt n=30,i,Istart,Iend,col[3],nev;
40: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
43: SlepcInitialize(&argc,&argv,(char*)0,help);
45: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
46: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem, n=%D\n\n",n);
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Compute the operator matrix that defines the eigensystem, Ax=kx
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: MatCreate(PETSC_COMM_WORLD,&A);
53: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
54: MatSetFromOptions(A);
55:
56: MatGetOwnershipRange(A,&Istart,&Iend);
57: if (Istart==0) FirstBlock=PETSC_TRUE;
58: if (Iend==n) LastBlock=PETSC_TRUE;
59: value[0]=-1.0; value[1]=2.0; value[2]=-1.0;
60: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
61: col[0]=i-1; col[1]=i; col[2]=i+1;
62: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
63: }
64: if (LastBlock) {
65: i=n-1; col[0]=n-2; col[1]=n-1;
66: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
67: }
68: if (FirstBlock) {
69: i=0; col[0]=0; col[1]=1; value[0]=2.0; value[1]=-1.0;
70: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
71: }
73: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
76: /* Create the shell version of A */
77: MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&A,&S);
78: MatSetFromOptions(S);
79: MatShellSetOperation(S,MATOP_MULT,(void(*)())MatMult_Shell);
80: MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_Shell);
81: MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Shell);
82: MatShellSetOperation(S,MATOP_SHIFT,(void(*)())MatShift_Shell);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create the eigensolver and set various options
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: EPSCreate(PETSC_COMM_WORLD,&eps);
89: EPSSetOperators(eps,S,PETSC_NULL);
90: EPSSetProblemType(eps,EPS_HEP);
91: EPSSetFromOptions(eps);
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Solve the eigensystem
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: EPSSolve(eps);
98: EPSGetType(eps,&type);
99: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
100: EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
101: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Display solution and clean up
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: EPSPrintSolution(eps,PETSC_NULL);
108: EPSDestroy(&eps);
109: MatDestroy(&A);
110: MatDestroy(&S);
111: SlepcFinalize();
112: return 0;
113: }
117: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
118: {
119: PetscErrorCode ierr;
120: Mat *A;
121:
123: MatShellGetContext(S,&A);
124: MatMult(*A,x,y);
125: return(0);
126: }
130: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
131: {
132: PetscErrorCode ierr;
133: Mat *A;
134:
136: MatShellGetContext(S,&A);
137: MatMultTranspose(*A,x,y);
138: return(0);
139: }
143: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
144: {
145: PetscErrorCode ierr;
146: Mat *A;
147:
149: MatShellGetContext(S,&A);
150: MatGetDiagonal(*A,diag);
151: return(0);
152: }
156: static PetscErrorCode MatShift_Shell(Mat S,PetscScalar shift)
157: {
158: PetscErrorCode ierr;
159: Mat *A;
160:
162: MatShellGetContext(S,&A);
163: MatShift(*A,shift);
164: return(0);
165: }
166: