Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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
12 :
13 : f(x) = (exp(x)-1)/x (the phi_1 function)
14 :
15 : with the following tree:
16 :
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 : */
23 :
24 : static char help[] = "Another test of a combined function.\n\n";
25 :
26 : #include <slepcfn.h>
27 :
28 : /*
29 : Compute matrix function B = A\(exp(A)-I)
30 : */
31 4 : PetscErrorCode TestMatCombine(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace)
32 : {
33 4 : PetscBool set,flg;
34 4 : PetscInt n;
35 4 : Mat F,Acopy;
36 4 : Vec v,f0;
37 4 : PetscReal nrm;
38 :
39 4 : PetscFunctionBeginUser;
40 4 : PetscCall(MatGetSize(A,&n,NULL));
41 4 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&F));
42 4 : PetscCall(PetscObjectSetName((PetscObject)F,"F"));
43 : /* compute matrix function */
44 4 : if (inplace) {
45 2 : PetscCall(MatCopy(A,F,SAME_NONZERO_PATTERN));
46 2 : PetscCall(MatIsHermitianKnown(A,&set,&flg));
47 2 : if (set && flg) PetscCall(MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE));
48 2 : PetscCall(FNEvaluateFunctionMat(fn,F,NULL));
49 : } else {
50 2 : PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Acopy));
51 2 : PetscCall(FNEvaluateFunctionMat(fn,A,F));
52 : /* check that A has not been modified */
53 2 : PetscCall(MatAXPY(Acopy,-1.0,A,SAME_NONZERO_PATTERN));
54 2 : PetscCall(MatNorm(Acopy,NORM_1,&nrm));
55 2 : if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the input matrix has changed by %g\n",(double)nrm));
56 2 : PetscCall(MatDestroy(&Acopy));
57 : }
58 4 : if (verbose) {
59 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n"));
60 0 : PetscCall(MatView(A,viewer));
61 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) - - - - - - -\n"));
62 0 : PetscCall(MatView(F,viewer));
63 : }
64 : /* print matrix norm for checking */
65 4 : PetscCall(MatNorm(F,NORM_1,&nrm));
66 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %6.3f\n",(double)nrm));
67 : /* check FNEvaluateFunctionMatVec() */
68 4 : PetscCall(MatCreateVecs(A,&v,&f0));
69 4 : PetscCall(MatGetColumnVector(F,f0,0));
70 4 : PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
71 4 : PetscCall(VecAXPY(v,-1.0,f0));
72 4 : PetscCall(VecNorm(v,NORM_2,&nrm));
73 4 : 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 4 : PetscCall(MatDestroy(&F));
75 4 : PetscCall(VecDestroy(&v));
76 4 : PetscCall(VecDestroy(&f0));
77 4 : PetscFunctionReturn(PETSC_SUCCESS);
78 : }
79 :
80 2 : int main(int argc,char **argv)
81 : {
82 2 : FN f,p,a,e,c,f1,f2;
83 2 : FNCombineType ctype;
84 2 : Mat A=NULL;
85 2 : PetscInt i,j,n=10,np;
86 2 : PetscScalar x,y,yp,*As,coeffs[10];
87 2 : char strx[50],str[50];
88 2 : PetscViewer viewer;
89 2 : PetscBool verbose,inplace,matcuda;
90 :
91 2 : PetscFunctionBeginUser;
92 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
93 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
94 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
95 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-inplace",&inplace));
96 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-matcuda",&matcuda));
97 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Phi1 via a combined function, n=%" PetscInt_FMT ".\n",n));
98 :
99 : /* Create function */
100 :
101 : /* e(x) = exp(x) */
102 2 : PetscCall(FNCreate(PETSC_COMM_WORLD,&e));
103 2 : PetscCall(PetscObjectSetName((PetscObject)e,"e"));
104 2 : PetscCall(FNSetType(e,FNEXP));
105 2 : PetscCall(FNSetFromOptions(e));
106 : /* c(x) = -1 */
107 2 : PetscCall(FNCreate(PETSC_COMM_WORLD,&c));
108 2 : PetscCall(PetscObjectSetName((PetscObject)c,"c"));
109 2 : PetscCall(FNSetType(c,FNRATIONAL));
110 2 : PetscCall(FNSetFromOptions(c));
111 2 : np = 1;
112 2 : coeffs[0] = -1.0;
113 2 : PetscCall(FNRationalSetNumerator(c,np,coeffs));
114 : /* a(x) */
115 2 : PetscCall(FNCreate(PETSC_COMM_WORLD,&a));
116 2 : PetscCall(PetscObjectSetName((PetscObject)a,"a"));
117 2 : PetscCall(FNSetType(a,FNCOMBINE));
118 2 : PetscCall(FNSetFromOptions(a));
119 2 : PetscCall(FNCombineSetChildren(a,FN_COMBINE_ADD,e,c));
120 : /* p(x) = x */
121 2 : PetscCall(FNCreate(PETSC_COMM_WORLD,&p));
122 2 : PetscCall(PetscObjectSetName((PetscObject)p,"p"));
123 2 : PetscCall(FNSetType(p,FNRATIONAL));
124 2 : PetscCall(FNSetFromOptions(p));
125 2 : np = 2;
126 2 : coeffs[0] = 1.0; coeffs[1] = 0.0;
127 2 : PetscCall(FNRationalSetNumerator(p,np,coeffs));
128 : /* f(x) */
129 2 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f));
130 2 : PetscCall(PetscObjectSetName((PetscObject)f,"f"));
131 2 : PetscCall(FNSetType(f,FNCOMBINE));
132 2 : PetscCall(FNSetFromOptions(f));
133 2 : PetscCall(FNCombineSetChildren(f,FN_COMBINE_DIVIDE,a,p));
134 :
135 : /* Set up viewer */
136 2 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
137 2 : PetscCall(FNCombineGetChildren(f,&ctype,&f1,&f2));
138 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Two functions combined with division:\n"));
139 2 : PetscCall(FNView(f1,viewer));
140 2 : PetscCall(FNView(f2,viewer));
141 2 : if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
142 :
143 : /* Scalar evaluation */
144 2 : x = 2.2;
145 2 : PetscCall(SlepcSNPrintfScalar(strx,sizeof(strx),x,PETSC_FALSE));
146 2 : PetscCall(FNEvaluateFunction(f,x,&y));
147 2 : PetscCall(FNEvaluateDerivative(f,x,&yp));
148 2 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),y,PETSC_FALSE));
149 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," f(%s)=%s\n",strx,str));
150 2 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),yp,PETSC_FALSE));
151 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," f'(%s)=%s\n",strx,str));
152 :
153 : /* Create matrices */
154 2 : if (matcuda) {
155 : #if defined(PETSC_HAVE_CUDA)
156 : PetscCall(MatCreateSeqDenseCUDA(PETSC_COMM_SELF,n,n,NULL,&A));
157 : #endif
158 2 : } else PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
159 2 : PetscCall(PetscObjectSetName((PetscObject)A,"A"));
160 :
161 : /* Fill A with 1-D Laplacian matrix */
162 2 : PetscCall(MatDenseGetArray(A,&As));
163 22 : for (i=0;i<n;i++) As[i+i*n]=2.0;
164 20 : j=1;
165 20 : for (i=0;i<n-j;i++) { As[i+(i+j)*n]=-1.0; As[(i+j)+i*n]=-1.0; }
166 2 : PetscCall(MatDenseRestoreArray(A,&As));
167 2 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
168 2 : PetscCall(TestMatCombine(f,A,viewer,verbose,inplace));
169 :
170 : /* Repeat with same matrix as non-symmetric */
171 2 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
172 2 : PetscCall(TestMatCombine(f,A,viewer,verbose,inplace));
173 :
174 2 : PetscCall(MatDestroy(&A));
175 2 : PetscCall(FNDestroy(&f));
176 2 : PetscCall(FNDestroy(&p));
177 2 : PetscCall(FNDestroy(&a));
178 2 : PetscCall(FNDestroy(&e));
179 2 : PetscCall(FNDestroy(&c));
180 2 : PetscCall(SlepcFinalize());
181 : return 0;
182 : }
183 :
184 : /*TEST
185 :
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
201 :
202 : TEST*/
|
