Actual source code: test10.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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: {
 21:   PetscScalar    y,yp;
 22:   char           strx[50],str[50];
 23:   Vec            v,f;
 24:   PetscReal      nrm;

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

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

 65:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 66:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 67:   PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
 68:   PetscPrintf(PETSC_COMM_WORLD,"Test Phi functions, n=%D.\n",n);

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

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

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

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

 98:   FNDuplicate(phik,PETSC_COMM_WORLD,&phicopy);
 99:   FNPhiGetIndex(phicopy,&k);
100:   PetscPrintf(PETSC_COMM_WORLD,"Index of phi function is %D\n",k);
101:   TestPhiFunction(phicopy,2.2,A,verbose);

103:   FNDestroy(&phi0);
104:   FNDestroy(&phi1);
105:   FNDestroy(&phik);
106:   FNDestroy(&phicopy);
107:   MatDestroy(&A);
108:   SlepcFinalize();
109:   return ierr;
110: }

112: /*TEST

114:    test:
115:       suffix: 1
116:       nsize: 1
117:       args: -fn_phi_index 3

119: TEST*/