Actual source code: test6.c
slepc-3.6.1 2015-09-03
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
21: /*
22: Define the function
24: f(x) = (1-x^2) exp( -x/(1+x^2) )
26: with the following tree:
28: f(x) f(x) (combined by product)
29: / \ g(x) = 1-x^2 (polynomial)
30: g(x) h(x) h(x) (combined by composition)
31: / \ r(x) = -x/(1+x^2) (rational)
32: r(x) e(x) e(x) = exp(x) (exponential)
33: */
35: static char help[] = "Test combined function.\n\n";
37: #include <slepcfn.h>
41: int main(int argc,char **argv)
42: {
44: FN f,g,h,e,r;
45: Mat A,B;
46: PetscInt i,j,n=10,np,nq;
47: PetscReal nrm;
48: PetscScalar x,y,yp,*As,p[10],q[10];
49: char strx[50],str[50];
50: PetscViewer viewer;
51: PetscBool verbose;
53: SlepcInitialize(&argc,&argv,(char*)0,help);
54: PetscOptionsGetInt(NULL,"-n",&n,NULL);
55: PetscOptionsHasName(NULL,"-verbose",&verbose);
56: PetscPrintf(PETSC_COMM_WORLD,"Combined function, n=%D.\n",n);
58: /* Create function */
60: /* e(x) = exp(x) */
61: FNCreate(PETSC_COMM_WORLD,&e);
62: FNSetType(e,FNEXP);
63: /* r(x) = x/(1+x^2) */
64: FNCreate(PETSC_COMM_WORLD,&r);
65: FNSetType(r,FNRATIONAL);
66: np = 2; nq = 3;
67: p[0] = -1.0; p[1] = 0.0;
68: q[0] = 1.0; q[1] = 0.0; q[2] = 1.0;
69: FNRationalSetNumerator(r,np,p);
70: FNRationalSetDenominator(r,nq,q);
71: /* h(x) */
72: FNCreate(PETSC_COMM_WORLD,&h);
73: FNSetType(h,FNCOMBINE);
74: FNCombineSetChildren(h,FN_COMBINE_COMPOSE,r,e);
75: /* g(x) = 1-x^2 */
76: FNCreate(PETSC_COMM_WORLD,&g);
77: FNSetType(g,FNRATIONAL);
78: np = 3;
79: p[0] = -1.0; p[1] = 0.0; p[2] = 1.0;
80: FNRationalSetNumerator(g,np,p);
81: /* f(x) */
82: FNCreate(PETSC_COMM_WORLD,&f);
83: FNSetType(f,FNCOMBINE);
84: FNCombineSetChildren(f,FN_COMBINE_MULTIPLY,g,h);
86: /* Set up viewer */
87: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
88: FNView(f,viewer);
89: if (verbose) {
90: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
91: }
93: /* Scalar evaluation */
94: x = 2.2;
95: SlepcSNPrintfScalar(strx,50,x,PETSC_FALSE);
96: FNEvaluateFunction(f,x,&y);
97: FNEvaluateDerivative(f,x,&yp);
98: SlepcSNPrintfScalar(str,50,y,PETSC_FALSE);
99: PetscPrintf(PETSC_COMM_WORLD," f(%s)=%s\n",strx,str);
100: SlepcSNPrintfScalar(str,50,yp,PETSC_FALSE);
101: PetscPrintf(PETSC_COMM_WORLD," f'(%s)=%s\n",strx,str);
103: /* Create matrices */
104: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
105: PetscObjectSetName((PetscObject)A,"A");
106: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&B);
107: PetscObjectSetName((PetscObject)B,"B");
109: /* Fill A with a symmetric Toeplitz matrix */
110: MatDenseGetArray(A,&As);
111: for (i=0;i<n;i++) As[i+i*n]=2.0;
112: for (j=1;j<3;j++) {
113: for (i=0;i<n-j;i++) { As[i+(i+j)*n]=1.0; As[(i+j)+i*n]=1.0; }
114: }
115: MatDenseRestoreArray(A,&As);
116: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
117: if (verbose) {
118: PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n");
119: MatView(A,viewer);
120: }
122: /* Evaluate matrix function */
123: FNEvaluateFunctionMat(f,A,B);
124: if (verbose) {
125: PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) - - - - - - -\n");
126: MatView(B,viewer);
127: }
128: MatNorm(B,NORM_1,&nrm);
129: PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %g\n",(double)nrm);
131: /* Repeat with same matrix as non-symmetric */
132: MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE);
134: /* Evaluate matrix function */
135: FNEvaluateFunctionMat(f,A,B);
136: if (verbose) {
137: PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) - - - - - - -\n");
138: MatView(B,viewer);
139: }
140: MatNorm(B,NORM_1,&nrm);
141: PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %g\n",(double)nrm);
143: MatDestroy(&A);
144: MatDestroy(&B);
145: FNDestroy(&f);
146: FNDestroy(&g);
147: FNDestroy(&h);
148: FNDestroy(&e);
149: FNDestroy(&r);
150: SlepcFinalize();
151: return 0;
152: }