Actual source code: test11.c
slepc-3.11.2 2019-07-30
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: */
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: {
34: PetscBool set,flg;
35: PetscInt n;
36: Mat F;
37: Vec v,f0;
38: PetscReal nrm;
41: MatGetSize(A,&n,NULL);
42: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&F);
43: PetscObjectSetName((PetscObject)F,"F");
44: /* compute matrix function */
45: if (inplace) {
46: MatCopy(A,F,SAME_NONZERO_PATTERN);
47: MatIsHermitianKnown(A,&set,&flg);
48: if (set && flg) { MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE); }
49: FNEvaluateFunctionMat(fn,F,NULL);
50: } else {
51: FNEvaluateFunctionMat(fn,A,F);
52: }
53: if (verbose) {
54: PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n");
55: MatView(A,viewer);
56: PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) - - - - - - -\n");
57: MatView(F,viewer);
58: }
59: /* print matrix norm for checking */
60: MatNorm(F,NORM_1,&nrm);
61: PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %6.3f\n",(double)nrm);
62: /* check FNEvaluateFunctionMatVec() */
63: MatCreateVecs(A,&v,&f0);
64: MatGetColumnVector(F,f0,0);
65: FNEvaluateFunctionMatVec(fn,A,v);
66: VecAXPY(v,-1.0,f0);
67: VecNorm(v,NORM_2,&nrm);
68: if (nrm>100*PETSC_MACHINE_EPSILON) {
69: PetscPrintf(PETSC_COMM_WORLD,"Warning: the norm of f(A)*e_1-v is %g\n",(double)nrm);
70: }
71: MatDestroy(&F);
72: VecDestroy(&v);
73: VecDestroy(&f0);
74: return(0);
75: }
77: int main(int argc,char **argv)
78: {
80: FN f,p,a,e,c,f1,f2;
81: FNCombineType ctype;
82: Mat A;
83: PetscInt i,j,n=10,np;
84: PetscScalar x,y,yp,*As,coeffs[10];
85: char strx[50],str[50];
86: PetscViewer viewer;
87: PetscBool verbose,inplace;
89: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
90: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
91: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
92: PetscOptionsHasName(NULL,NULL,"-inplace",&inplace);
93: PetscPrintf(PETSC_COMM_WORLD,"Phi1 via a combined function, n=%D.\n",n);
95: /* Create function */
97: /* e(x) = exp(x) */
98: FNCreate(PETSC_COMM_WORLD,&e);
99: PetscObjectSetName((PetscObject)e,"e");
100: FNSetType(e,FNEXP);
101: FNSetFromOptions(e);
102: /* c(x) = -1 */
103: FNCreate(PETSC_COMM_WORLD,&c);
104: PetscObjectSetName((PetscObject)c,"c");
105: FNSetType(c,FNRATIONAL);
106: FNSetFromOptions(c);
107: np = 1;
108: coeffs[0] = -1.0;
109: FNRationalSetNumerator(c,np,coeffs);
110: /* a(x) */
111: FNCreate(PETSC_COMM_WORLD,&a);
112: PetscObjectSetName((PetscObject)a,"a");
113: FNSetType(a,FNCOMBINE);
114: FNSetFromOptions(a);
115: FNCombineSetChildren(a,FN_COMBINE_ADD,e,c);
116: /* p(x) = x */
117: FNCreate(PETSC_COMM_WORLD,&p);
118: PetscObjectSetName((PetscObject)p,"p");
119: FNSetType(p,FNRATIONAL);
120: FNSetFromOptions(p);
121: np = 2;
122: coeffs[0] = 1.0; coeffs[1] = 0.0;
123: FNRationalSetNumerator(p,np,coeffs);
124: /* f(x) */
125: FNCreate(PETSC_COMM_WORLD,&f);
126: PetscObjectSetName((PetscObject)f,"f");
127: FNSetType(f,FNCOMBINE);
128: FNSetFromOptions(f);
129: FNCombineSetChildren(f,FN_COMBINE_DIVIDE,a,p);
131: /* Set up viewer */
132: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
133: FNCombineGetChildren(f,&ctype,&f1,&f2);
134: PetscPrintf(PETSC_COMM_WORLD,"Two functions combined with division:\n");
135: FNView(f1,viewer);
136: FNView(f2,viewer);
137: if (verbose) {
138: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
139: }
141: /* Scalar evaluation */
142: x = 2.2;
143: SlepcSNPrintfScalar(strx,50,x,PETSC_FALSE);
144: FNEvaluateFunction(f,x,&y);
145: FNEvaluateDerivative(f,x,&yp);
146: SlepcSNPrintfScalar(str,50,y,PETSC_FALSE);
147: PetscPrintf(PETSC_COMM_WORLD," f(%s)=%s\n",strx,str);
148: SlepcSNPrintfScalar(str,50,yp,PETSC_FALSE);
149: PetscPrintf(PETSC_COMM_WORLD," f'(%s)=%s\n",strx,str);
151: /* Create matrices */
152: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
153: PetscObjectSetName((PetscObject)A,"A");
155: /* Fill A with 1-D Laplacian matrix */
156: MatDenseGetArray(A,&As);
157: for (i=0;i<n;i++) As[i+i*n]=2.0;
158: j=1;
159: for (i=0;i<n-j;i++) { As[i+(i+j)*n]=-1.0; As[(i+j)+i*n]=-1.0; }
160: MatDenseRestoreArray(A,&As);
161: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
162: TestMatCombine(f,A,viewer,verbose,inplace);
164: /* Repeat with same matrix as non-symmetric */
165: MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE);
166: TestMatCombine(f,A,viewer,verbose,inplace);
168: MatDestroy(&A);
169: FNDestroy(&f);
170: FNDestroy(&p);
171: FNDestroy(&a);
172: FNDestroy(&e);
173: FNDestroy(&c);
174: SlepcFinalize();
175: return ierr;
176: }
178: /*TEST
180: test:
181: suffix: 1
182: nsize: 1
184: test:
185: suffix: 2
186: nsize: 1
187: args: -inplace
188: output_file: output/test11_1.out
190: TEST*/