Actual source code: test4.c
slepc-3.21.2 2024-09-25
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[] = "Illustrates use of MFNSetBV().\n\n"
12: "The command line options are:\n"
13: " -t <sval>, where <sval> = scalar value that multiplies the argument.\n"
14: " -file <filename>, where <filename> = matrix file in PETSc binary form.\n\n";
16: #include <slepcmfn.h>
18: int main(int argc,char **argv)
19: {
20: Mat A; /* problem matrix */
21: MFN mfn;
22: FN f;
23: BV V;
24: PetscInt k=8;
25: PetscReal norm;
26: PetscScalar t=2.0;
27: Vec v,y;
28: PetscViewer viewer;
29: PetscBool flg;
30: char filename[PETSC_MAX_PATH_LEN];
31: MFNConvergedReason reason;
33: PetscFunctionBeginUser;
34: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
36: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
37: PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
38: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMatrix exponential y=exp(t*A)*e, loaded from file\n\n"));
40: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: Load matrix A from binary file
42: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44: PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
45: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name with the -file option");
47: #if defined(PETSC_USE_COMPLEX)
48: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
49: #else
50: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
51: #endif
52: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
53: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
54: PetscCall(MatSetFromOptions(A));
55: PetscCall(MatLoad(A,viewer));
56: PetscCall(PetscViewerDestroy(&viewer));
58: /* set v = ones(n,1) */
59: PetscCall(MatCreateVecs(A,&v,&y));
60: PetscCall(VecSet(v,1.0));
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create the BV object
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: PetscCall(BVCreate(PETSC_COMM_WORLD,&V));
67: PetscCall(PetscObjectSetName((PetscObject)V,"V"));
68: PetscCall(BVSetSizesFromVec(V,v,k));
69: PetscCall(BVSetFromOptions(V));
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Create the solver and set various options
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
76: PetscCall(MFNSetOperator(mfn,A));
77: PetscCall(MFNSetBV(mfn,V));
78: PetscCall(MFNGetFN(mfn,&f));
79: PetscCall(FNSetType(f,FNEXP));
80: PetscCall(FNSetScale(f,t,1.0));
81: PetscCall(MFNSetDimensions(mfn,k));
82: PetscCall(MFNSetFromOptions(mfn));
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Solve the problem, y=exp(t*A)*v
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: PetscCall(MFNSolve(mfn,v,y));
89: PetscCall(MFNGetConvergedReason(mfn,&reason));
90: PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
91: PetscCall(VecNorm(y,NORM_2,&norm));
92: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n\n",(double)PetscRealPart(t),(double)norm));
94: /*
95: Free work space
96: */
97: PetscCall(BVDestroy(&V));
98: PetscCall(MFNDestroy(&mfn));
99: PetscCall(MatDestroy(&A));
100: PetscCall(VecDestroy(&v));
101: PetscCall(VecDestroy(&y));
102: PetscCall(SlepcFinalize());
103: return 0;
104: }
106: /*TEST
108: test:
109: suffix: 1
110: args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -k 12
111: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
113: test:
114: suffix: 2
115: args: -file ${DATAFILESPATH}/matrices/complex/qc324.petsc -k 12
116: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
118: TEST*/