Actual source code: test1.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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Test ST with shell matrices as problem matrices.\n\n";
24: #include <slepceps.h>
26: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag);
27: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y);
28: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y);
29: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M);
33: static PetscErrorCode MyShellMatCreate(Mat *A,Mat *M)
34: {
36: MPI_Comm comm;
37: PetscInt n;
40: MatGetSize(*A,&n,NULL);
41: PetscObjectGetComm((PetscObject)*A,&comm);
42: MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,n,n,A,M);
43: MatSetFromOptions(*M);
44: MatShellSetOperation(*M,MATOP_MULT,(void(*)())MatMult_Shell);
45: MatShellSetOperation(*M,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_Shell);
46: MatShellSetOperation(*M,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Shell);
47: MatShellSetOperation(*M,MATOP_DUPLICATE,(void(*)())MatDuplicate_Shell);
48: return(0);
49: }
53: int main(int argc,char **argv)
54: {
55: Mat A,S; /* problem matrix */
56: EPS eps; /* eigenproblem solver context */
57: EPSType type;
58: PetscScalar value[3];
59: PetscInt n=30,i,Istart,Iend,col[3],nev;
60: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
63: SlepcInitialize(&argc,&argv,(char*)0,help);
65: PetscOptionsGetInt(NULL,"-n",&n,NULL);
66: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem, n=%D\n\n",n);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Compute the operator matrix that defines the eigensystem, Ax=kx
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: MatCreate(PETSC_COMM_WORLD,&A);
73: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
74: MatSetFromOptions(A);
75: MatSetUp(A);
77: MatGetOwnershipRange(A,&Istart,&Iend);
78: if (Istart==0) FirstBlock=PETSC_TRUE;
79: if (Iend==n) LastBlock=PETSC_TRUE;
80: value[0]=-1.0; value[1]=2.0; value[2]=-1.0;
81: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
82: col[0]=i-1; col[1]=i; col[2]=i+1;
83: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
84: }
85: if (LastBlock) {
86: i=n-1; col[0]=n-2; col[1]=n-1;
87: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
88: }
89: if (FirstBlock) {
90: i=0; col[0]=0; col[1]=1; value[0]=2.0; value[1]=-1.0;
91: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
92: }
94: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
97: /* Create the shell version of A */
98: MyShellMatCreate(&A,&S);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create the eigensolver and set various options
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: EPSCreate(PETSC_COMM_WORLD,&eps);
105: EPSSetOperators(eps,S,NULL);
106: EPSSetProblemType(eps,EPS_HEP);
107: EPSSetFromOptions(eps);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Solve the eigensystem
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: EPSSolve(eps);
114: EPSGetType(eps,&type);
115: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
116: EPSGetDimensions(eps,&nev,NULL,NULL);
117: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Display solution and clean up
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: EPSPrintSolution(eps,NULL);
124: EPSDestroy(&eps);
125: MatDestroy(&A);
126: MatDestroy(&S);
127: SlepcFinalize();
128: return 0;
129: }
133: static PetscErrorCode MatMult_Shell(Mat S,Vec x,Vec y)
134: {
135: PetscErrorCode ierr;
136: Mat *A;
139: MatShellGetContext(S,&A);
140: MatMult(*A,x,y);
141: return(0);
142: }
146: static PetscErrorCode MatMultTranspose_Shell(Mat S,Vec x,Vec y)
147: {
148: PetscErrorCode ierr;
149: Mat *A;
152: MatShellGetContext(S,&A);
153: MatMultTranspose(*A,x,y);
154: return(0);
155: }
159: static PetscErrorCode MatGetDiagonal_Shell(Mat S,Vec diag)
160: {
161: PetscErrorCode ierr;
162: Mat *A;
165: MatShellGetContext(S,&A);
166: MatGetDiagonal(*A,diag);
167: return(0);
168: }
172: static PetscErrorCode MatDuplicate_Shell(Mat S,MatDuplicateOption op,Mat *M)
173: {
175: Mat *A;
178: MatShellGetContext(S,&A);
179: MyShellMatCreate(A,M);
180: return(0);
181: }