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[] = "Test matrix function evaluation via diagonalization.\n\n";
      12             : 
      13             : #include <slepcfn.h>
      14             : 
      15           3 : int main(int argc,char **argv)
      16             : {
      17           3 :   FN             fn;
      18           3 :   Mat            A,F,G;
      19           3 :   PetscInt       i,j,n=10;
      20           3 :   PetscReal      nrm;
      21           3 :   PetscScalar    *As,alpha,beta;
      22           3 :   PetscViewer    viewer;
      23           3 :   PetscBool      verbose;
      24             : 
      25           3 :   PetscFunctionBeginUser;
      26           3 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      27           3 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      28           3 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
      29           3 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix function of symmetric/Hermitian matrix, n=%" PetscInt_FMT ".\n",n));
      30             : 
      31             :   /* Create function object */
      32           3 :   PetscCall(FNCreate(PETSC_COMM_WORLD,&fn));
      33           3 :   PetscCall(FNSetType(fn,FNEXP));   /* default to exponential */
      34             : #if defined(PETSC_USE_COMPLEX)
      35           3 :   alpha = PetscCMPLX(0.3,0.8);
      36           3 :   beta  = PetscCMPLX(1.1,-0.1);
      37             : #else
      38             :   alpha = 0.3;
      39             :   beta  = 1.1;
      40             : #endif
      41           3 :   PetscCall(FNSetScale(fn,alpha,beta));
      42           3 :   PetscCall(FNSetFromOptions(fn));
      43             : 
      44             :   /* Set up viewer */
      45           3 :   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
      46           3 :   if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
      47             : 
      48             :   /* Create a symmetric/Hermitian Toeplitz matrix */
      49           3 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
      50           3 :   PetscCall(PetscObjectSetName((PetscObject)A,"A"));
      51           3 :   PetscCall(MatDenseGetArray(A,&As));
      52          33 :   for (i=0;i<n;i++) As[i+i*n]=2.0;
      53           9 :   for (j=1;j<3;j++) {
      54          57 :     for (i=0;i<n-j;i++) {
      55             : #if defined(PETSC_USE_COMPLEX)
      56          51 :       As[i+(i+j)*n]=PetscCMPLX(1.0,0.1); As[(i+j)+i*n]=PetscCMPLX(1.0,-0.1);
      57             : #else
      58             :       As[i+(i+j)*n]=0.5; As[(i+j)+i*n]=0.5;
      59             : #endif
      60             :     }
      61             :   }
      62           3 :   PetscCall(MatDenseRestoreArray(A,&As));
      63           3 :   if (verbose) {
      64           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n"));
      65           0 :     PetscCall(MatView(A,viewer));
      66             :   }
      67             : 
      68             :   /* compute matrix function */
      69           3 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&F));
      70           3 :   PetscCall(PetscObjectSetName((PetscObject)F,"F"));
      71           3 :   PetscCall(FNEvaluateFunctionMat(fn,A,F));
      72           3 :   if (verbose) {
      73           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) - - - - - - -\n"));
      74           0 :     PetscCall(MatView(F,viewer));
      75             :   }
      76             : 
      77             :   /* Repeat with MAT_HERMITIAN flag set */
      78           3 :   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
      79           3 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&G));
      80           3 :   PetscCall(PetscObjectSetName((PetscObject)G,"G"));
      81           3 :   PetscCall(FNEvaluateFunctionMat(fn,A,G));
      82           3 :   if (verbose) {
      83           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) symm - - - - - - -\n"));
      84           0 :     PetscCall(MatView(G,viewer));
      85             :   }
      86             : 
      87             :   /* compare the two results */
      88           3 :   PetscCall(MatAXPY(F,-1.0,G,SAME_NONZERO_PATTERN));
      89           3 :   PetscCall(MatNorm(F,NORM_FROBENIUS,&nrm));
      90           3 :   if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the norm of F-G is %g\n",(double)nrm));
      91           3 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed results match.\n"));
      92             : 
      93           3 :   PetscCall(MatDestroy(&A));
      94           3 :   PetscCall(MatDestroy(&F));
      95           3 :   PetscCall(MatDestroy(&G));
      96           3 :   PetscCall(FNDestroy(&fn));
      97           3 :   PetscCall(SlepcFinalize());
      98             :   return 0;
      99             : }
     100             : 
     101             : /*TEST
     102             : 
     103             :    test:
     104             :       suffix: 1
     105             :       nsize: 1
     106             :       args: -fn_type {{exp sqrt}shared output}
     107             :       output_file: output/test12_1.out
     108             : 
     109             :    test:
     110             :       suffix: 1_rational
     111             :       nsize: 1
     112             :       args: -fn_type rational -fn_rational_numerator 2,-1.5 -fn_rational_denominator 1,0.8
     113             :       output_file: output/test12_1.out
     114             : 
     115             : TEST*/