Actual source code: test1.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[] = "Computes exp(t*A)*v for a matrix loaded from file.\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: PetscReal norm;
24: PetscScalar t=2.0;
25: Vec v,y;
26: PetscViewer viewer;
27: PetscBool flg;
28: char filename[PETSC_MAX_PATH_LEN];
29: MFNConvergedReason reason;
31: PetscFunctionBeginUser;
32: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
34: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
35: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMatrix exponential y=exp(t*A)*e, loaded from file\n\n"));
37: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: Load matrix A from binary file
39: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41: PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
42: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name with the -file option");
44: #if defined(PETSC_USE_COMPLEX)
45: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
46: #else
47: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
48: #endif
49: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
50: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
51: PetscCall(MatSetFromOptions(A));
52: PetscCall(MatLoad(A,viewer));
53: PetscCall(PetscViewerDestroy(&viewer));
55: /* set v = ones(n,1) */
56: PetscCall(MatCreateVecs(A,&v,&y));
57: PetscCall(VecSet(v,1.0));
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Create the solver and set various options
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
64: PetscCall(MFNSetOperator(mfn,A));
65: PetscCall(MFNGetFN(mfn,&f));
66: PetscCall(FNSetType(f,FNEXP));
67: PetscCall(FNSetScale(f,t,1.0));
68: PetscCall(MFNSetFromOptions(mfn));
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Solve the problem, y=exp(t*A)*v
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: PetscCall(MFNSolve(mfn,v,y));
75: PetscCall(MFNGetConvergedReason(mfn,&reason));
76: PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
77: PetscCall(VecNorm(y,NORM_2,&norm));
78: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n",(double)PetscRealPart(t),(double)norm));
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Solve the problem, y=exp(t*A^T)*v
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: PetscCall(MFNSolveTranspose(mfn,v,y));
85: PetscCall(MFNGetConvergedReason(mfn,&reason));
86: PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
87: PetscCall(VecNorm(y,NORM_2,&norm));
88: PetscCall(PetscPrintf(PETSC_COMM_WORLD," With transpose: computed vector has norm %g\n\n",(double)norm));
90: /*
91: Free work space
92: */
93: PetscCall(MFNDestroy(&mfn));
94: PetscCall(MatDestroy(&A));
95: PetscCall(VecDestroy(&v));
96: PetscCall(VecDestroy(&y));
97: PetscCall(SlepcFinalize());
98: return 0;
99: }
101: /*TEST
103: testset:
104: args: -file ${DATAFILESPATH}/matrices/real/bfw782a.petsc -mfn_type {{krylov expokit}} -t 0.05
105: requires: double !complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
106: output_file: output/test1_1.out
107: test:
108: suffix: 1
109: test:
110: suffix: 1_cuda
111: args: -mat_type aijcusparse
112: requires: cuda
114: testset:
115: args: -file ${DATAFILESPATH}/matrices/complex/qc324.petsc -mfn_type {{krylov expokit}}
116: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
117: output_file: output/test1_2.out
118: test:
119: suffix: 2
120: test:
121: suffix: 2_cuda
122: args: -mat_type aijcusparse
123: requires: cuda
125: TEST*/