Actual source code: test10.c

slepc-3.18.2 2023-01-26
Report Typos and Errors
  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;

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

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

 67:   /* Create matrix, fill it with 1-D Laplacian */
 68:   MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
 69:   PetscObjectSetName((PetscObject)A,"A");
 70:   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:   MatDenseRestoreArray(A,&As);

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

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

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

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

100:   FNDestroy(&phi0);
101:   FNDestroy(&phi1);
102:   FNDestroy(&phik);
103:   FNDestroy(&phicopy);
104:   MatDestroy(&A);
105:   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*/