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