Actual source code: test6.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) = (1-x^2) exp( -x/(1+x^2) )
15: with the following tree:
17: f(x) f(x) (combined by product)
18: / \ g(x) = 1-x^2 (polynomial)
19: g(x) h(x) h(x) (combined by composition)
20: / \ r(x) = -x/(1+x^2) (rational)
21: r(x) e(x) e(x) = exp(x) (exponential)
22: */
24: static char help[] = "Test combined function.\n\n";
26: #include <slepcfn.h>
28: /*
29: Compute matrix function B = (I-A^2) exp( -(I+A^2)\A )
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,g,h,e,r,fcopy;
81: Mat A;
82: PetscInt i,j,n=10,np,nq;
83: PetscScalar x,y,yp,*As,p[10],q[10];
84: char strx[50],str[50];
85: PetscViewer viewer;
86: PetscBool verbose,inplace;
88: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
89: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
90: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
91: PetscOptionsHasName(NULL,NULL,"-inplace",&inplace);
92: PetscPrintf(PETSC_COMM_WORLD,"Combined function, n=%D.\n",n);
94: /* Create function */
96: /* e(x) = exp(x) */
97: FNCreate(PETSC_COMM_WORLD,&e);
98: FNSetType(e,FNEXP);
99: FNSetFromOptions(e);
100: /* r(x) = x/(1+x^2) */
101: FNCreate(PETSC_COMM_WORLD,&r);
102: FNSetType(r,FNRATIONAL);
103: FNSetFromOptions(r);
104: np = 2; nq = 3;
105: p[0] = -1.0; p[1] = 0.0;
106: q[0] = 1.0; q[1] = 0.0; q[2] = 1.0;
107: FNRationalSetNumerator(r,np,p);
108: FNRationalSetDenominator(r,nq,q);
109: /* h(x) */
110: FNCreate(PETSC_COMM_WORLD,&h);
111: FNSetType(h,FNCOMBINE);
112: FNSetFromOptions(h);
113: FNCombineSetChildren(h,FN_COMBINE_COMPOSE,r,e);
114: /* g(x) = 1-x^2 */
115: FNCreate(PETSC_COMM_WORLD,&g);
116: FNSetType(g,FNRATIONAL);
117: FNSetFromOptions(g);
118: np = 3;
119: p[0] = -1.0; p[1] = 0.0; p[2] = 1.0;
120: FNRationalSetNumerator(g,np,p);
121: /* f(x) */
122: FNCreate(PETSC_COMM_WORLD,&f);
123: FNSetType(f,FNCOMBINE);
124: FNSetFromOptions(f);
125: FNCombineSetChildren(f,FN_COMBINE_MULTIPLY,g,h);
127: /* Set up viewer */
128: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
129: FNView(f,viewer);
130: if (verbose) {
131: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
132: }
134: /* Scalar evaluation */
135: x = 2.2;
136: SlepcSNPrintfScalar(strx,50,x,PETSC_FALSE);
137: FNEvaluateFunction(f,x,&y);
138: FNEvaluateDerivative(f,x,&yp);
139: SlepcSNPrintfScalar(str,50,y,PETSC_FALSE);
140: PetscPrintf(PETSC_COMM_WORLD," f(%s)=%s\n",strx,str);
141: SlepcSNPrintfScalar(str,50,yp,PETSC_FALSE);
142: PetscPrintf(PETSC_COMM_WORLD," f'(%s)=%s\n",strx,str);
144: /* Test duplication */
145: FNDuplicate(f,PetscObjectComm((PetscObject)f),&fcopy);
147: /* Create matrices */
148: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
149: PetscObjectSetName((PetscObject)A,"A");
151: /* Fill A with a symmetric Toeplitz matrix */
152: MatDenseGetArray(A,&As);
153: for (i=0;i<n;i++) As[i+i*n]=2.0;
154: for (j=1;j<3;j++) {
155: for (i=0;i<n-j;i++) { As[i+(i+j)*n]=1.0; As[(i+j)+i*n]=1.0; }
156: }
157: MatDenseRestoreArray(A,&As);
158: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
159: TestMatCombine(fcopy,A,viewer,verbose,inplace);
161: /* Repeat with same matrix as non-symmetric */
162: MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE);
163: TestMatCombine(fcopy,A,viewer,verbose,inplace);
165: MatDestroy(&A);
166: FNDestroy(&f);
167: FNDestroy(&fcopy);
168: FNDestroy(&g);
169: FNDestroy(&h);
170: FNDestroy(&e);
171: FNDestroy(&r);
172: SlepcFinalize();
173: return ierr;
174: }
176: /*TEST
178: test:
179: suffix: 1
180: nsize: 1
182: test:
183: suffix: 2
184: nsize: 1
185: args: -inplace
186: output_file: output/test6_1.out
188: TEST*/