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 Phi functions.\n\n";
      12             : 
      13             : #include <slepcfn.h>
      14             : 
      15             : /*
      16             :    Evaluates phi_k function on a scalar and on a matrix
      17             :  */
      18           3 : PetscErrorCode TestPhiFunction(FN fn,PetscScalar x,Mat A,PetscBool verbose)
      19             : {
      20           3 :   PetscScalar    y,yp;
      21           3 :   char           strx[50],str[50];
      22           3 :   Vec            v,f;
      23           3 :   PetscReal      nrm;
      24             : 
      25           3 :   PetscFunctionBeginUser;
      26           3 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
      27           3 :   PetscCall(FNView(fn,NULL));
      28           3 :   PetscCall(SlepcSNPrintfScalar(strx,sizeof(strx),x,PETSC_FALSE));
      29           3 :   PetscCall(FNEvaluateFunction(fn,x,&y));
      30           3 :   PetscCall(FNEvaluateDerivative(fn,x,&yp));
      31           3 :   PetscCall(SlepcSNPrintfScalar(str,sizeof(str),y,PETSC_FALSE));
      32           3 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nf(%s)=%s\n",strx,str));
      33           3 :   PetscCall(SlepcSNPrintfScalar(str,sizeof(str),yp,PETSC_FALSE));
      34           3 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"f'(%s)=%s\n",strx,str));
      35             :   /* compute phi_k(A)*e_1 */
      36           3 :   PetscCall(MatCreateVecs(A,&v,&f));
      37           3 :   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
      38           3 :   PetscCall(FNEvaluateFunctionMatVec(fn,A,f));  /* reference result by diagonalization */
      39           3 :   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
      40           3 :   PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
      41           3 :   PetscCall(VecAXPY(v,-1.0,f));
      42           3 :   PetscCall(VecNorm(v,NORM_2,&nrm));
      43           3 :   if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the norm of f(A)*e_1-ref is %g\n",(double)nrm));
      44           3 :   if (verbose) {
      45           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"f(A)*e_1 =\n"));
      46           0 :     PetscCall(VecView(v,NULL));
      47             :   }
      48           3 :   PetscCall(VecDestroy(&v));
      49           3 :   PetscCall(VecDestroy(&f));
      50           3 :   PetscFunctionReturn(PETSC_SUCCESS);
      51             : }
      52             : 
      53           1 : int main(int argc,char **argv)
      54             : {
      55           1 :   FN             phi0,phi1,phik,phicopy;
      56           1 :   Mat            A;
      57           1 :   PetscInt       i,j,n=8,k;
      58           1 :   PetscScalar    tau,eta,*As;
      59           1 :   PetscBool      verbose;
      60             : 
      61           1 :   PetscFunctionBeginUser;
      62           1 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      63           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      64           1 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
      65           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test Phi functions, n=%" PetscInt_FMT ".\n",n));
      66             : 
      67             :   /* Create matrix, fill it with 1-D Laplacian */
      68           1 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
      69           1 :   PetscCall(PetscObjectSetName((PetscObject)A,"A"));
      70           1 :   PetscCall(MatDenseGetArray(A,&As));
      71           9 :   for (i=0;i<n;i++) As[i+i*n]=2.0;
      72           8 :   j=1;
      73           8 :   for (i=0;i<n-j;i++) { As[i+(i+j)*n]=-1.0; As[(i+j)+i*n]=-1.0; }
      74           1 :   PetscCall(MatDenseRestoreArray(A,&As));
      75             : 
      76             :   /* phi_0(x) = exp(x) */
      77           1 :   PetscCall(FNCreate(PETSC_COMM_WORLD,&phi0));
      78           1 :   PetscCall(FNSetType(phi0,FNPHI));
      79           1 :   PetscCall(FNPhiSetIndex(phi0,0));
      80           1 :   PetscCall(TestPhiFunction(phi0,2.2,A,verbose));
      81             : 
      82             :   /* phi_1(x) = (exp(x)-1)/x with scaling factors eta*phi_1(tau*x) */
      83           1 :   PetscCall(FNCreate(PETSC_COMM_WORLD,&phi1));
      84           1 :   PetscCall(FNSetType(phi1,FNPHI));  /* default index should be 1 */
      85           1 :   tau = 0.2;
      86           1 :   eta = 1.3;
      87           1 :   PetscCall(FNSetScale(phi1,tau,eta));
      88           1 :   PetscCall(TestPhiFunction(phi1,2.2,A,verbose));
      89             : 
      90             :   /* phi_k(x) with index set from command-line arguments */
      91           1 :   PetscCall(FNCreate(PETSC_COMM_WORLD,&phik));
      92           1 :   PetscCall(FNSetType(phik,FNPHI));
      93           1 :   PetscCall(FNSetFromOptions(phik));
      94             : 
      95           1 :   PetscCall(FNDuplicate(phik,PETSC_COMM_WORLD,&phicopy));
      96           1 :   PetscCall(FNPhiGetIndex(phicopy,&k));
      97           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Index of phi function is %" PetscInt_FMT "\n",k));
      98           1 :   PetscCall(TestPhiFunction(phicopy,2.2,A,verbose));
      99             : 
     100           1 :   PetscCall(FNDestroy(&phi0));
     101           1 :   PetscCall(FNDestroy(&phi1));
     102           1 :   PetscCall(FNDestroy(&phik));
     103           1 :   PetscCall(FNDestroy(&phicopy));
     104           1 :   PetscCall(MatDestroy(&A));
     105           1 :   PetscCall(SlepcFinalize());
     106             :   return 0;
     107             : }
     108             : 
     109             : /*TEST
     110             : 
     111             :    test:
     112             :       suffix: 1
     113             :       nsize: 1
     114             :       args: -fn_phi_index 3
     115             :       requires: !single
     116             : 
     117             : TEST*/