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[] = "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";
      15             : 
      16             : #include <slepcmfn.h>
      17             : 
      18           1 : int main(int argc,char **argv)
      19             : {
      20           1 :   Mat                A;           /* problem matrix */
      21           1 :   MFN                mfn;
      22           1 :   FN                 f;
      23           1 :   BV                 V;
      24           1 :   PetscInt           k=8;
      25           1 :   PetscReal          norm;
      26           1 :   PetscScalar        t=2.0;
      27           1 :   Vec                v,y;
      28           1 :   PetscViewer        viewer;
      29           1 :   PetscBool          flg;
      30           1 :   char               filename[PETSC_MAX_PATH_LEN];
      31           1 :   MFNConvergedReason reason;
      32             : 
      33           1 :   PetscFunctionBeginUser;
      34           1 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      35             : 
      36           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
      37           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
      38           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMatrix exponential y=exp(t*A)*e, loaded from file\n\n"));
      39             : 
      40             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      41             :                 Load matrix A from binary file
      42             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      43             : 
      44           1 :   PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
      45           1 :   PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name with the -file option");
      46             : 
      47             : #if defined(PETSC_USE_COMPLEX)
      48           1 :   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           1 :   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
      53           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      54           1 :   PetscCall(MatSetFromOptions(A));
      55           1 :   PetscCall(MatLoad(A,viewer));
      56           1 :   PetscCall(PetscViewerDestroy(&viewer));
      57             : 
      58             :   /* set v = ones(n,1) */
      59           1 :   PetscCall(MatCreateVecs(A,&v,&y));
      60           1 :   PetscCall(VecSet(v,1.0));
      61             : 
      62             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      63             :                         Create the BV object
      64             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      65             : 
      66           1 :   PetscCall(BVCreate(PETSC_COMM_WORLD,&V));
      67           1 :   PetscCall(PetscObjectSetName((PetscObject)V,"V"));
      68           1 :   PetscCall(BVSetSizesFromVec(V,v,k));
      69           1 :   PetscCall(BVSetFromOptions(V));
      70             : 
      71             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      72             :                 Create the solver and set various options
      73             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      74             : 
      75           1 :   PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
      76           1 :   PetscCall(MFNSetOperator(mfn,A));
      77           1 :   PetscCall(MFNSetBV(mfn,V));
      78           1 :   PetscCall(MFNGetFN(mfn,&f));
      79           1 :   PetscCall(FNSetType(f,FNEXP));
      80           1 :   PetscCall(FNSetScale(f,t,1.0));
      81           1 :   PetscCall(MFNSetDimensions(mfn,k));
      82           1 :   PetscCall(MFNSetFromOptions(mfn));
      83             : 
      84             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      85             :                       Solve the problem, y=exp(t*A)*v
      86             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      87             : 
      88           1 :   PetscCall(MFNSolve(mfn,v,y));
      89           1 :   PetscCall(MFNGetConvergedReason(mfn,&reason));
      90           1 :   PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
      91           1 :   PetscCall(VecNorm(y,NORM_2,&norm));
      92           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n\n",(double)PetscRealPart(t),(double)norm));
      93             : 
      94             :   /*
      95             :      Free work space
      96             :   */
      97           1 :   PetscCall(BVDestroy(&V));
      98           1 :   PetscCall(MFNDestroy(&mfn));
      99           1 :   PetscCall(MatDestroy(&A));
     100           1 :   PetscCall(VecDestroy(&v));
     101           1 :   PetscCall(VecDestroy(&y));
     102           1 :   PetscCall(SlepcFinalize());
     103             :   return 0;
     104             : }
     105             : 
     106             : /*TEST
     107             : 
     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)
     112             : 
     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)
     117             : 
     118             : TEST*/