Actual source code: test11.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: */
     10: /*
     11:    Define the function

     13:         f(x) = (exp(x)-1)/x    (the phi_1 function)

     15:    with the following tree:

     17:             f(x)                  f(x)              (combined by division)
     18:            /    \                 p(x) = x          (polynomial)
     19:         a(x)    p(x)              a(x)              (combined by addition)
     20:        /    \                     e(x) = exp(x)     (exponential)
     21:      e(x)   c(x)                  c(x) = -1         (constant)
     22: */

     24: static char help[] = "Another test of a combined function.\n\n";

     26: #include <slepcfn.h>

     28: /*
     29:    Compute matrix function B = A\(exp(A)-I)
     30:  */
     31: PetscErrorCode TestMatCombine(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace)
     32: {
     33:   PetscBool      set,flg;
     34:   PetscInt       n;
     35:   Mat            F,Acopy;
     36:   Vec            v,f0;
     37:   PetscReal      nrm;

     39:   PetscFunctionBeginUser;
     40:   PetscCall(MatGetSize(A,&n,NULL));
     41:   PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&F));
     42:   PetscCall(PetscObjectSetName((PetscObject)F,"F"));
     43:   /* compute matrix function */
     44:   if (inplace) {
     45:     PetscCall(MatCopy(A,F,SAME_NONZERO_PATTERN));
     46:     PetscCall(MatIsHermitianKnown(A,&set,&flg));
     47:     if (set && flg) PetscCall(MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE));
     48:     PetscCall(FNEvaluateFunctionMat(fn,F,NULL));
     49:   } else {
     50:     PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Acopy));
     51:     PetscCall(FNEvaluateFunctionMat(fn,A,F));
     52:     /* check that A has not been modified */
     53:     PetscCall(MatAXPY(Acopy,-1.0,A,SAME_NONZERO_PATTERN));
     54:     PetscCall(MatNorm(Acopy,NORM_1,&nrm));
     55:     if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the input matrix has changed by %g\n",(double)nrm));
     56:     PetscCall(MatDestroy(&Acopy));
     57:   }
     58:   if (verbose) {
     59:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n"));
     60:     PetscCall(MatView(A,viewer));
     61:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) - - - - - - -\n"));
     62:     PetscCall(MatView(F,viewer));
     63:   }
     64:   /* print matrix norm for checking */
     65:   PetscCall(MatNorm(F,NORM_1,&nrm));
     66:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %6.3f\n",(double)nrm));
     67:   /* check FNEvaluateFunctionMatVec() */
     68:   PetscCall(MatCreateVecs(A,&v,&f0));
     69:   PetscCall(MatGetColumnVector(F,f0,0));
     70:   PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
     71:   PetscCall(VecAXPY(v,-1.0,f0));
     72:   PetscCall(VecNorm(v,NORM_2,&nrm));
     73:   if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the norm of f(A)*e_1-v is %g\n",(double)nrm));
     74:   PetscCall(MatDestroy(&F));
     75:   PetscCall(VecDestroy(&v));
     76:   PetscCall(VecDestroy(&f0));
     77:   PetscFunctionReturn(PETSC_SUCCESS);
     78: }

     80: int main(int argc,char **argv)
     81: {
     82:   FN             f,p,a,e,c,f1,f2;
     83:   FNCombineType  ctype;
     84:   Mat            A=NULL;
     85:   PetscInt       i,j,n=10,np;
     86:   PetscScalar    x,y,yp,*As,coeffs[10];
     87:   char           strx[50],str[50];
     88:   PetscViewer    viewer;
     89:   PetscBool      verbose,inplace,matcuda;

     91:   PetscFunctionBeginUser;
     92:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     93:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     94:   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
     95:   PetscCall(PetscOptionsHasName(NULL,NULL,"-inplace",&inplace));
     96:   PetscCall(PetscOptionsHasName(NULL,NULL,"-matcuda",&matcuda));
     97:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Phi1 via a combined function, n=%" PetscInt_FMT ".\n",n));

     99:   /* Create function */

    101:   /* e(x) = exp(x) */
    102:   PetscCall(FNCreate(PETSC_COMM_WORLD,&e));
    103:   PetscCall(PetscObjectSetName((PetscObject)e,"e"));
    104:   PetscCall(FNSetType(e,FNEXP));
    105:   PetscCall(FNSetFromOptions(e));
    106:   /* c(x) = -1 */
    107:   PetscCall(FNCreate(PETSC_COMM_WORLD,&c));
    108:   PetscCall(PetscObjectSetName((PetscObject)c,"c"));
    109:   PetscCall(FNSetType(c,FNRATIONAL));
    110:   PetscCall(FNSetFromOptions(c));
    111:   np = 1;
    112:   coeffs[0] = -1.0;
    113:   PetscCall(FNRationalSetNumerator(c,np,coeffs));
    114:   /* a(x) */
    115:   PetscCall(FNCreate(PETSC_COMM_WORLD,&a));
    116:   PetscCall(PetscObjectSetName((PetscObject)a,"a"));
    117:   PetscCall(FNSetType(a,FNCOMBINE));
    118:   PetscCall(FNSetFromOptions(a));
    119:   PetscCall(FNCombineSetChildren(a,FN_COMBINE_ADD,e,c));
    120:   /* p(x) = x */
    121:   PetscCall(FNCreate(PETSC_COMM_WORLD,&p));
    122:   PetscCall(PetscObjectSetName((PetscObject)p,"p"));
    123:   PetscCall(FNSetType(p,FNRATIONAL));
    124:   PetscCall(FNSetFromOptions(p));
    125:   np = 2;
    126:   coeffs[0] = 1.0; coeffs[1] = 0.0;
    127:   PetscCall(FNRationalSetNumerator(p,np,coeffs));
    128:   /* f(x) */
    129:   PetscCall(FNCreate(PETSC_COMM_WORLD,&f));
    130:   PetscCall(PetscObjectSetName((PetscObject)f,"f"));
    131:   PetscCall(FNSetType(f,FNCOMBINE));
    132:   PetscCall(FNSetFromOptions(f));
    133:   PetscCall(FNCombineSetChildren(f,FN_COMBINE_DIVIDE,a,p));

    135:   /* Set up viewer */
    136:   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
    137:   PetscCall(FNCombineGetChildren(f,&ctype,&f1,&f2));
    138:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Two functions combined with division:\n"));
    139:   PetscCall(FNView(f1,viewer));
    140:   PetscCall(FNView(f2,viewer));
    141:   if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));

    143:   /* Scalar evaluation */
    144:   x = 2.2;
    145:   PetscCall(SlepcSNPrintfScalar(strx,sizeof(strx),x,PETSC_FALSE));
    146:   PetscCall(FNEvaluateFunction(f,x,&y));
    147:   PetscCall(FNEvaluateDerivative(f,x,&yp));
    148:   PetscCall(SlepcSNPrintfScalar(str,sizeof(str),y,PETSC_FALSE));
    149:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  f(%s)=%s\n",strx,str));
    150:   PetscCall(SlepcSNPrintfScalar(str,sizeof(str),yp,PETSC_FALSE));
    151:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  f'(%s)=%s\n",strx,str));

    153:   /* Create matrices */
    154:   if (matcuda) {
    155: #if defined(PETSC_HAVE_CUDA)
    156:     PetscCall(MatCreateSeqDenseCUDA(PETSC_COMM_SELF,n,n,NULL,&A));
    157: #endif
    158:   } else PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
    159:   PetscCall(PetscObjectSetName((PetscObject)A,"A"));

    161:   /* Fill A with 1-D Laplacian matrix */
    162:   PetscCall(MatDenseGetArray(A,&As));
    163:   for (i=0;i<n;i++) As[i+i*n]=2.0;
    164:   j=1;
    165:   for (i=0;i<n-j;i++) { As[i+(i+j)*n]=-1.0; As[(i+j)+i*n]=-1.0; }
    166:   PetscCall(MatDenseRestoreArray(A,&As));
    167:   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
    168:   PetscCall(TestMatCombine(f,A,viewer,verbose,inplace));

    170:   /* Repeat with same matrix as non-symmetric */
    171:   PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
    172:   PetscCall(TestMatCombine(f,A,viewer,verbose,inplace));

    174:   PetscCall(MatDestroy(&A));
    175:   PetscCall(FNDestroy(&f));
    176:   PetscCall(FNDestroy(&p));
    177:   PetscCall(FNDestroy(&a));
    178:   PetscCall(FNDestroy(&e));
    179:   PetscCall(FNDestroy(&c));
    180:   PetscCall(SlepcFinalize());
    181:   return 0;
    182: }

    184: /*TEST

    186:    testset:
    187:       output_file: output/test11_1.out
    188:       test:
    189:          suffix: 1
    190:       test:
    191:          suffix: 1_cuda
    192:          args: -matcuda
    193:          requires: cuda
    194:       test:
    195:          suffix: 2
    196:          args: -inplace
    197:       test:
    198:          suffix: 2_cuda
    199:          args: -inplace -matcuda
    200:          requires: cuda

    202: TEST*/