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 : static char help[] = "Test matrix exponential.\n\n";
12 :
13 : #include <slepcfn.h>
14 :
15 : /*
16 : Compute matrix exponential B = expm(A)
17 : */
18 40 : PetscErrorCode TestMatExp(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace,PetscBool checkerror)
19 : {
20 40 : PetscScalar tau,eta;
21 40 : PetscBool set,flg;
22 40 : PetscInt n;
23 40 : Mat F,R,Finv,Acopy;
24 40 : Vec v,f0;
25 40 : FN finv;
26 40 : PetscReal nrm,nrmf;
27 :
28 40 : PetscFunctionBeginUser;
29 40 : PetscCall(MatGetSize(A,&n,NULL));
30 40 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&F));
31 40 : PetscCall(PetscObjectSetName((PetscObject)F,"F"));
32 : /* compute matrix exponential */
33 40 : if (inplace) {
34 8 : PetscCall(MatCopy(A,F,SAME_NONZERO_PATTERN));
35 8 : PetscCall(MatIsHermitianKnown(A,&set,&flg));
36 8 : if (set && flg) PetscCall(MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE));
37 8 : PetscCall(FNEvaluateFunctionMat(fn,F,NULL));
38 : } else {
39 32 : PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Acopy));
40 32 : PetscCall(FNEvaluateFunctionMat(fn,A,F));
41 : /* check that A has not been modified */
42 32 : PetscCall(MatAXPY(Acopy,-1.0,A,SAME_NONZERO_PATTERN));
43 32 : PetscCall(MatNorm(Acopy,NORM_1,&nrm));
44 32 : if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the input matrix has changed by %g\n",(double)nrm));
45 32 : PetscCall(MatDestroy(&Acopy));
46 : }
47 40 : if (verbose) {
48 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n"));
49 0 : PetscCall(MatView(A,viewer));
50 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed expm(A) - - - - - - -\n"));
51 0 : PetscCall(MatView(F,viewer));
52 : }
53 : /* print matrix norm for checking */
54 40 : PetscCall(MatNorm(F,NORM_1,&nrmf));
55 40 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %g\n",(double)nrmf));
56 40 : if (checkerror) {
57 0 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Finv));
58 0 : PetscCall(PetscObjectSetName((PetscObject)Finv,"Finv"));
59 0 : PetscCall(FNGetScale(fn,&tau,&eta));
60 : /* compute inverse exp(-tau*A)/eta */
61 0 : PetscCall(FNCreate(PETSC_COMM_WORLD,&finv));
62 0 : PetscCall(FNSetType(finv,FNEXP));
63 0 : PetscCall(FNSetFromOptions(finv));
64 0 : PetscCall(FNSetScale(finv,-tau,1.0/eta));
65 0 : if (inplace) {
66 0 : PetscCall(MatCopy(A,Finv,SAME_NONZERO_PATTERN));
67 0 : PetscCall(MatIsHermitianKnown(A,&set,&flg));
68 0 : if (set && flg) PetscCall(MatSetOption(Finv,MAT_HERMITIAN,PETSC_TRUE));
69 0 : PetscCall(FNEvaluateFunctionMat(finv,Finv,NULL));
70 0 : } else PetscCall(FNEvaluateFunctionMat(finv,A,Finv));
71 0 : PetscCall(FNDestroy(&finv));
72 : /* check error ||F*Finv-I||_F */
73 0 : PetscCall(MatMatMult(F,Finv,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&R));
74 0 : PetscCall(MatShift(R,-1.0));
75 0 : PetscCall(MatNorm(R,NORM_FROBENIUS,&nrm));
76 0 : if (nrm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||exp(A)*exp(-A)-I||_F < 100*eps\n"));
77 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||exp(A)*exp(-A)-I||_F = %g\n",(double)nrm));
78 0 : PetscCall(MatDestroy(&R));
79 0 : PetscCall(MatDestroy(&Finv));
80 : }
81 : /* check FNEvaluateFunctionMatVec() */
82 40 : PetscCall(MatCreateVecs(A,&v,&f0));
83 40 : PetscCall(MatGetColumnVector(F,f0,0));
84 40 : PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
85 40 : PetscCall(VecAXPY(v,-1.0,f0));
86 40 : PetscCall(VecNorm(v,NORM_2,&nrm));
87 40 : if (nrm/nrmf>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the norm of f(A)*e_1-v is %g\n",(double)nrm));
88 40 : PetscCall(MatDestroy(&F));
89 40 : PetscCall(VecDestroy(&v));
90 40 : PetscCall(VecDestroy(&f0));
91 40 : PetscFunctionReturn(PETSC_SUCCESS);
92 : }
93 :
94 20 : int main(int argc,char **argv)
95 : {
96 20 : FN fn;
97 20 : Mat A=NULL;
98 20 : PetscInt i,j,n=10;
99 20 : PetscScalar *As;
100 20 : PetscViewer viewer;
101 20 : PetscBool verbose,inplace,checkerror,matcuda;
102 :
103 20 : PetscFunctionBeginUser;
104 20 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
105 20 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
106 20 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
107 20 : PetscCall(PetscOptionsHasName(NULL,NULL,"-inplace",&inplace));
108 20 : PetscCall(PetscOptionsHasName(NULL,NULL,"-checkerror",&checkerror));
109 20 : PetscCall(PetscOptionsHasName(NULL,NULL,"-matcuda",&matcuda));
110 20 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix exponential, n=%" PetscInt_FMT ".\n",n));
111 :
112 : /* Create exponential function object */
113 20 : PetscCall(FNCreate(PETSC_COMM_WORLD,&fn));
114 20 : PetscCall(FNSetType(fn,FNEXP));
115 20 : PetscCall(FNSetFromOptions(fn));
116 :
117 : /* Set up viewer */
118 20 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
119 20 : PetscCall(FNView(fn,viewer));
120 20 : if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
121 :
122 : /* Create matrices */
123 20 : if (matcuda) {
124 : #if defined(PETSC_HAVE_CUDA)
125 : PetscCall(MatCreateSeqDenseCUDA(PETSC_COMM_SELF,n,n,NULL,&A));
126 : #endif
127 20 : } else PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
128 20 : PetscCall(PetscObjectSetName((PetscObject)A,"A"));
129 :
130 : /* Fill A with a symmetric Toeplitz matrix */
131 20 : PetscCall(MatDenseGetArray(A,&As));
132 660 : for (i=0;i<n;i++) As[i+i*n]=2.0;
133 60 : for (j=1;j<3;j++) {
134 1260 : for (i=0;i<n-j;i++) { As[i+(i+j)*n]=1.0; As[(i+j)+i*n]=1.0; }
135 : }
136 20 : PetscCall(MatDenseRestoreArray(A,&As));
137 20 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
138 20 : PetscCall(TestMatExp(fn,A,viewer,verbose,inplace,checkerror));
139 :
140 : /* Repeat with non-symmetric A */
141 20 : PetscCall(MatDenseGetArray(A,&As));
142 60 : for (j=1;j<3;j++) {
143 1260 : for (i=0;i<n-j;i++) { As[(i+j)+i*n]=-1.0; }
144 : }
145 20 : PetscCall(MatDenseRestoreArray(A,&As));
146 20 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
147 20 : PetscCall(TestMatExp(fn,A,viewer,verbose,inplace,checkerror));
148 :
149 20 : PetscCall(MatDestroy(&A));
150 20 : PetscCall(FNDestroy(&fn));
151 20 : PetscCall(SlepcFinalize());
152 : return 0;
153 : }
154 :
155 : /*TEST
156 :
157 : testset:
158 : filter: grep -v "computing matrix functions"
159 : output_file: output/test3_1.out
160 : test:
161 : suffix: 1
162 : args: -fn_method {{0 1}}
163 : test:
164 : suffix: 1_subdiagonalpade
165 : args: -fn_method {{2 3}}
166 : requires: c99_complex !single
167 : test:
168 : suffix: 1_cuda
169 : args: -fn_method 1 -matcuda
170 : requires: cuda !magma
171 : test:
172 : suffix: 1_magma
173 : args: -fn_method {{0 1 2 3}} -matcuda
174 : requires: cuda magma
175 : test:
176 : suffix: 2
177 : args: -inplace -fn_method{{0 1}}
178 : test:
179 : suffix: 2_subdiagonalpade
180 : args: -inplace -fn_method{{2 3}}
181 : requires: c99_complex !single
182 : test:
183 : suffix: 2_cuda
184 : args: -inplace -fn_method 1 -matcuda
185 : requires: cuda !magma
186 : test:
187 : suffix: 2_magma
188 : args: -inplace -fn_method {{0 1 2 3}} -matcuda
189 : requires: cuda magma
190 :
191 : testset:
192 : args: -fn_scale 0.1
193 : filter: grep -v "computing matrix functions"
194 : output_file: output/test3_3.out
195 : test:
196 : suffix: 3
197 : args: -fn_method {{0 1}}
198 : test:
199 : suffix: 3_subdiagonalpade
200 : args: -fn_method {{2 3}}
201 : requires: c99_complex !single
202 :
203 : testset:
204 : args: -n 120 -fn_scale 0.6,1.5
205 : filter: grep -v "computing matrix functions"
206 : output_file: output/test3_4.out
207 : test:
208 : suffix: 4
209 : args: -fn_method {{0 1}}
210 : requires: !single
211 : test:
212 : suffix: 4_subdiagonalpade
213 : args: -fn_method {{2 3}}
214 : requires: c99_complex !single
215 :
216 : test:
217 : suffix: 5
218 : args: -fn_scale 30 -fn_method {{2 3}}
219 : filter: grep -v "computing matrix functions"
220 : requires: c99_complex !single
221 : output_file: output/test3_5.out
222 :
223 : test:
224 : suffix: 6
225 : args: -fn_scale 1e-9 -fn_method {{2 3}}
226 : filter: grep -v "computing matrix functions"
227 : requires: c99_complex !single
228 : output_file: output/test3_6.out
229 :
230 : TEST*/
|
