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,(char*)0,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*/