Actual source code: test10.c

      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[] = "Test Phi functions.\n\n";

     13: #include <slepcfn.h>

     15: /*
     16:    Evaluates phi_k function on a scalar and on a matrix
     17:  */
     18: PetscErrorCode TestPhiFunction(FN fn,PetscScalar x,Mat A,PetscBool verbose)
     19: {
     20:   PetscScalar    y,yp;
     21:   char           strx[50],str[50];
     22:   Vec            v,f;
     23:   PetscReal      nrm;

     25:   PetscFunctionBeginUser;
     26:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
     27:   PetscCall(FNView(fn,NULL));
     28:   PetscCall(SlepcSNPrintfScalar(strx,sizeof(strx),x,PETSC_FALSE));
     29:   PetscCall(FNEvaluateFunction(fn,x,&y));
     30:   PetscCall(FNEvaluateDerivative(fn,x,&yp));
     31:   PetscCall(SlepcSNPrintfScalar(str,sizeof(str),y,PETSC_FALSE));
     32:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nf(%s)=%s\n",strx,str));
     33:   PetscCall(SlepcSNPrintfScalar(str,sizeof(str),yp,PETSC_FALSE));
     34:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"f'(%s)=%s\n",strx,str));
     35:   /* compute phi_k(A)*e_1 */
     36:   PetscCall(MatCreateVecs(A,&v,&f));
     37:   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
     38:   PetscCall(FNEvaluateFunctionMatVec(fn,A,f));  /* reference result by diagonalization */
     39:   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
     40:   PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
     41:   PetscCall(VecAXPY(v,-1.0,f));
     42:   PetscCall(VecNorm(v,NORM_2,&nrm));
     43:   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:   if (verbose) {
     45:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"f(A)*e_1 =\n"));
     46:     PetscCall(VecView(v,NULL));
     47:   }
     48:   PetscCall(VecDestroy(&v));
     49:   PetscCall(VecDestroy(&f));
     50:   PetscFunctionReturn(PETSC_SUCCESS);
     51: }

     53: int main(int argc,char **argv)
     54: {
     55:   FN             phi0,phi1,phik,phicopy;
     56:   Mat            A;
     57:   PetscInt       i,j,n=8,k;
     58:   PetscScalar    tau,eta,*As;
     59:   PetscBool      verbose;

     61:   PetscFunctionBeginUser;
     62:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     63:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     64:   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
     65:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test Phi functions, n=%" PetscInt_FMT ".\n",n));

     67:   /* Create matrix, fill it with 1-D Laplacian */
     68:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
     69:   PetscCall(PetscObjectSetName((PetscObject)A,"A"));
     70:   PetscCall(MatDenseGetArray(A,&As));
     71:   for (i=0;i<n;i++) As[i+i*n]=2.0;
     72:   j=1;
     73:   for (i=0;i<n-j;i++) { As[i+(i+j)*n]=-1.0; As[(i+j)+i*n]=-1.0; }
     74:   PetscCall(MatDenseRestoreArray(A,&As));

     76:   /* phi_0(x) = exp(x) */
     77:   PetscCall(FNCreate(PETSC_COMM_WORLD,&phi0));
     78:   PetscCall(FNSetType(phi0,FNPHI));
     79:   PetscCall(FNPhiSetIndex(phi0,0));
     80:   PetscCall(TestPhiFunction(phi0,2.2,A,verbose));

     82:   /* phi_1(x) = (exp(x)-1)/x with scaling factors eta*phi_1(tau*x) */
     83:   PetscCall(FNCreate(PETSC_COMM_WORLD,&phi1));
     84:   PetscCall(FNSetType(phi1,FNPHI));  /* default index should be 1 */
     85:   tau = 0.2;
     86:   eta = 1.3;
     87:   PetscCall(FNSetScale(phi1,tau,eta));
     88:   PetscCall(TestPhiFunction(phi1,2.2,A,verbose));

     90:   /* phi_k(x) with index set from command-line arguments */
     91:   PetscCall(FNCreate(PETSC_COMM_WORLD,&phik));
     92:   PetscCall(FNSetType(phik,FNPHI));
     93:   PetscCall(FNSetFromOptions(phik));

     95:   PetscCall(FNDuplicate(phik,PETSC_COMM_WORLD,&phicopy));
     96:   PetscCall(FNPhiGetIndex(phicopy,&k));
     97:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Index of phi function is %" PetscInt_FMT "\n",k));
     98:   PetscCall(TestPhiFunction(phicopy,2.2,A,verbose));

    100:   PetscCall(FNDestroy(&phi0));
    101:   PetscCall(FNDestroy(&phi1));
    102:   PetscCall(FNDestroy(&phik));
    103:   PetscCall(FNDestroy(&phicopy));
    104:   PetscCall(MatDestroy(&A));
    105:   PetscCall(SlepcFinalize());
    106:   return 0;
    107: }

    109: /*TEST

    111:    test:
    112:       suffix: 1
    113:       nsize: 1
    114:       args: -fn_phi_index 3
    115:       requires: !single

    117: TEST*/