Actual source code: test3.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 matrix exponential.\n\n";

 13: #include <slepcfn.h>

 15: /*
 16:    Compute matrix exponential B = expm(A)
 17:  */
 18: PetscErrorCode TestMatExp(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace)
 19: {
 21:   PetscBool      set,flg;
 22:   PetscInt       n;
 23:   Mat            F;
 24:   Vec            v,f0;
 25:   PetscReal      nrm;

 28:   MatGetSize(A,&n,NULL);
 29:   MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&F);
 30:   PetscObjectSetName((PetscObject)F,"F");
 31:   /* compute matrix exponential */
 32:   if (inplace) {
 33:     MatCopy(A,F,SAME_NONZERO_PATTERN);
 34:     MatIsHermitianKnown(A,&set,&flg);
 35:     if (set && flg) { MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE); }
 36:     FNEvaluateFunctionMat(fn,F,NULL);
 37:   } else {
 38:     FNEvaluateFunctionMat(fn,A,F);
 39:   }
 40:   if (verbose) {
 41:     PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n");
 42:     MatView(A,viewer);
 43:     PetscPrintf(PETSC_COMM_WORLD,"Computed expm(A) - - - - - - -\n");
 44:     MatView(F,viewer);
 45:   }
 46:   /* print matrix norm for checking */
 47:   MatNorm(F,NORM_1,&nrm);
 48:   PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %g\n",(double)nrm);
 49:   /* check FNEvaluateFunctionMatVec() */
 50:   MatCreateVecs(A,&v,&f0);
 51:   MatGetColumnVector(F,f0,0);
 52:   FNEvaluateFunctionMatVec(fn,A,v);
 53:   VecAXPY(v,-1.0,f0);
 54:   VecNorm(v,NORM_2,&nrm);
 55:   if (nrm>100*PETSC_MACHINE_EPSILON) {
 56:     PetscPrintf(PETSC_COMM_WORLD,"Warning: the norm of f(A)*e_1-v is %g\n",(double)nrm);
 57:   }
 58:   MatDestroy(&F);
 59:   VecDestroy(&v);
 60:   VecDestroy(&f0);
 61:   return(0);
 62: }

 64: int main(int argc,char **argv)
 65: {
 67:   FN             fn;
 68:   Mat            A;
 69:   PetscInt       i,j,n=10;
 70:   PetscScalar    *As;
 71:   PetscViewer    viewer;
 72:   PetscBool      verbose,inplace;

 74:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 75:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 76:   PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
 77:   PetscOptionsHasName(NULL,NULL,"-inplace",&inplace);
 78:   PetscPrintf(PETSC_COMM_WORLD,"Matrix exponential, n=%D.\n",n);

 80:   /* Create exponential function object */
 81:   FNCreate(PETSC_COMM_WORLD,&fn);
 82:   FNSetType(fn,FNEXP);
 83:   FNSetFromOptions(fn);

 85:   /* Set up viewer */
 86:   PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
 87:   FNView(fn,viewer);
 88:   if (verbose) {
 89:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
 90:   }

 92:   /* Create matrices */
 93:   MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
 94:   PetscObjectSetName((PetscObject)A,"A");

 96:   /* Fill A with a symmetric Toeplitz matrix */
 97:   MatDenseGetArray(A,&As);
 98:   for (i=0;i<n;i++) As[i+i*n]=2.0;
 99:   for (j=1;j<3;j++) {
100:     for (i=0;i<n-j;i++) { As[i+(i+j)*n]=1.0; As[(i+j)+i*n]=1.0; }
101:   }
102:   MatDenseRestoreArray(A,&As);
103:   MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
104:   TestMatExp(fn,A,viewer,verbose,inplace);

106:   /* Repeat with non-symmetric A */
107:   MatDenseGetArray(A,&As);
108:   for (j=1;j<3;j++) {
109:     for (i=0;i<n-j;i++) { As[(i+j)+i*n]=-1.0; }
110:   }
111:   MatDenseRestoreArray(A,&As);
112:   MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE);
113:   TestMatExp(fn,A,viewer,verbose,inplace);

115:   MatDestroy(&A);
116:   FNDestroy(&fn);
117:   SlepcFinalize();
118:   return ierr;
119: }

121: /*TEST

123:    test:
124:       suffix: 1
125:       nsize: 1
126:       args: -fn_method {{0 1}shared output}
127:       filter: grep -v "computing matrix functions"
128:       output_file: output/test3_1.out

130:    test:
131:       suffix: 1_subdiagonalpade
132:       nsize: 1
133:       args: -fn_method {{2 3}}
134:       requires: c99_complex !single
135:       filter: grep -v "computing matrix functions"
136:       output_file: output/test3_1.out

138:    test:
139:       suffix: 2
140:       nsize: 1
141:       args: -inplace
142:       filter: grep -v "computing matrix functions"
143:       output_file: output/test3_1.out

145:    test:
146:       suffix: 3
147:       nsize: 1
148:       args: -fn_scale 0.1 -fn_method {{0 1}shared output}
149:       filter: grep -v "computing matrix functions"
150:       output_file: output/test3_3.out

152:    testset:
153:       nsize: 1
154:       filter: grep -v "computing matrix functions"
155:       output_file: output/test3_3.out
156:       test:
157:         suffix: 3_subdiagonalpade_product
158:         args: -fn_scale 0.1 -fn_method 2
159:         requires: c99_complex
160:       test:
161:         suffix: 3_subdiagonalpade_partial
162:         args: -fn_scale 0.1 -fn_method 3
163:         requires: c99_complex !single

165:    test:
166:       suffix: 4
167:       nsize: 1
168:       args: -n 80 -fn_scale 0.6,1.5 -fn_method {{0 1}shared output}
169:       requires: !single
170:       filter: grep -v "computing matrix functions"
171:       output_file: output/test3_4.out

173:    testset:
174:       nsize: 1
175:       filter: grep -v "computing matrix functions"
176:       output_file: output/test3_4.out
177:       test:
178:         suffix: 4_subdiagonalpade_product
179:         args: -n 80 -fn_scale 0.6,1.5 -fn_method 2
180:         requires: c99_complex
181:       test:
182:         suffix: 4_subdiagonalpade_partial
183:         args: -n 80 -fn_scale 0.6,1.5 -fn_method 3
184:         requires: c99_complex !single

186: TEST*/