Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      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";
      15             : 
      16             : #include <slepcmfn.h>
      17             : 
      18           2 : int main(int argc,char **argv)
      19             : {
      20           2 :   Mat                A;           /* problem matrix */
      21           2 :   MFN                mfn;
      22           2 :   FN                 f;
      23           2 :   PetscReal          norm;
      24           2 :   PetscScalar        t=2.0;
      25           2 :   Vec                v,y;
      26           2 :   PetscViewer        viewer;
      27           2 :   PetscBool          flg;
      28           2 :   char               filename[PETSC_MAX_PATH_LEN];
      29           2 :   MFNConvergedReason reason;
      30             : 
      31           2 :   PetscFunctionBeginUser;
      32           2 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      33             : 
      34           2 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
      35           2 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMatrix exponential y=exp(t*A)*e, loaded from file\n\n"));
      36             : 
      37             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      38             :                 Load matrix A from binary file
      39             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      40             : 
      41           2 :   PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
      42           2 :   PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name with the -file option");
      43             : 
      44             : #if defined(PETSC_USE_COMPLEX)
      45           2 :   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           2 :   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
      50           2 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      51           2 :   PetscCall(MatSetFromOptions(A));
      52           2 :   PetscCall(MatLoad(A,viewer));
      53           2 :   PetscCall(PetscViewerDestroy(&viewer));
      54             : 
      55             :   /* set v = ones(n,1) */
      56           2 :   PetscCall(MatCreateVecs(A,&v,&y));
      57           2 :   PetscCall(VecSet(v,1.0));
      58             : 
      59             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      60             :                 Create the solver and set various options
      61             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      62             : 
      63           2 :   PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
      64           2 :   PetscCall(MFNSetOperator(mfn,A));
      65           2 :   PetscCall(MFNGetFN(mfn,&f));
      66           2 :   PetscCall(FNSetType(f,FNEXP));
      67           2 :   PetscCall(FNSetScale(f,t,1.0));
      68           2 :   PetscCall(MFNSetFromOptions(mfn));
      69             : 
      70             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      71             :                       Solve the problem, y=exp(t*A)*v
      72             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      73             : 
      74           2 :   PetscCall(MFNSolve(mfn,v,y));
      75           2 :   PetscCall(MFNGetConvergedReason(mfn,&reason));
      76           2 :   PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
      77           2 :   PetscCall(VecNorm(y,NORM_2,&norm));
      78           2 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n",(double)PetscRealPart(t),(double)norm));
      79             : 
      80             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      81             :                       Solve the problem, y=exp(t*A^T)*v
      82             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      83             : 
      84           2 :   PetscCall(MFNSolveTranspose(mfn,v,y));
      85           2 :   PetscCall(MFNGetConvergedReason(mfn,&reason));
      86           2 :   PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
      87           2 :   PetscCall(VecNorm(y,NORM_2,&norm));
      88           2 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," With transpose: computed vector has norm %g\n\n",(double)norm));
      89             : 
      90             :   /*
      91             :      Free work space
      92             :   */
      93           2 :   PetscCall(MFNDestroy(&mfn));
      94           2 :   PetscCall(MatDestroy(&A));
      95           2 :   PetscCall(VecDestroy(&v));
      96           2 :   PetscCall(VecDestroy(&y));
      97           2 :   PetscCall(SlepcFinalize());
      98             :   return 0;
      99             : }
     100             : 
     101             : /*TEST
     102             : 
     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
     113             : 
     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
     124             : 
     125             : TEST*/