Actual source code: test1.c
slepc-3.7.2 2016-07-19
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[] = "Computes exp(A)*v for a matrix loaded from file.\n\n"
23: "The command line options are:\n"
24: " -file <filename>, where <filename> = matrix file in PETSc binary form.\n\n";
26: #include <slepcmfn.h>
30: int main(int argc,char **argv)
31: {
32: Mat A; /* problem matrix */
33: MFN mfn;
34: FN f;
35: PetscReal norm;
36: PetscScalar t=2.0;
37: Vec v,y;
38: PetscErrorCode ierr;
39: PetscViewer viewer;
40: PetscBool flg;
41: char filename[PETSC_MAX_PATH_LEN];
42: MFNConvergedReason reason;
44: SlepcInitialize(&argc,&argv,(char*)0,help);
46: PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL);
47: PetscPrintf(PETSC_COMM_WORLD,"\nMatrix exponential y=exp(t*A)*e, loaded from file\n\n");
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Load matrix A from binary file
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: PetscOptionsGetString(NULL,NULL,"-file",filename,PETSC_MAX_PATH_LEN,&flg);
54: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name with the -file option");
56: #if defined(PETSC_USE_COMPLEX)
57: PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");
58: #else
59: PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");
60: #endif
61: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
62: MatCreate(PETSC_COMM_WORLD,&A);
63: MatSetFromOptions(A);
64: MatLoad(A,viewer);
65: PetscViewerDestroy(&viewer);
67: /* set v = ones(n,1) */
68: MatCreateVecs(A,NULL,&y);
69: MatCreateVecs(A,NULL,&v);
70: VecSet(v,1.0);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create the solver and set various options
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: MFNCreate(PETSC_COMM_WORLD,&mfn);
77: MFNSetOperator(mfn,A);
78: MFNGetFN(mfn,&f);
79: FNSetType(f,FNEXP);
80: FNSetScale(f,t,1.0);
81: MFNSetFromOptions(mfn);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Solve the problem, y=exp(t*A)*v
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: MFNSolve(mfn,v,y);
88: MFNGetConvergedReason(mfn,&reason);
89: if (reason<0) SETERRQ(PETSC_COMM_WORLD,1,"Solver did not converge");
90: VecNorm(y,NORM_2,&norm);
91: PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n\n",(double)PetscRealPart(t),(double)norm);
92:
93: /*
94: Free work space
95: */
96: MFNDestroy(&mfn);
97: MatDestroy(&A);
98: VecDestroy(&v);
99: VecDestroy(&y);
100: SlepcFinalize();
101: return ierr;
102: }