Actual source code: test3.c
slepc-3.15.0 2021-03-31
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: testset:
157: filter: grep -v "computing matrix functions"
158: output_file: output/test3_1.out
159: test:
160: suffix: 1
161: args: -fn_method {{0 1}}
162: test:
163: suffix: 1_subdiagonalpade
164: args: -fn_method {{2 3}}
165: requires: c99_complex !single
166: test:
167: suffix: 1_cuda
168: args: -fn_method 4
169: requires: cuda
170: test:
171: suffix: 1_magma
172: args: -fn_method {{5 6 7 8}}
173: requires: cuda magma
174: test:
175: suffix: 2
176: args: -inplace -fn_method{{0 1}}
177: test:
178: suffix: 2_subdiagonalpade
179: args: -inplace -fn_method{{2 3}}
180: requires: c99_complex !single
181: test:
182: suffix: 2_cuda
183: args: -inplace -fn_method 4
184: requires: cuda
185: test:
186: suffix: 2_magma
187: args: -inplace -fn_method {{5 6 7 8}}
188: requires: cuda magma
190: testset:
191: args: -fn_scale 0.1
192: filter: grep -v "computing matrix functions"
193: output_file: output/test3_3.out
194: test:
195: suffix: 3
196: args: -fn_method {{0 1}}
197: test:
198: suffix: 3_subdiagonalpade
199: args: -fn_method {{2 3}}
200: requires: c99_complex !single
202: testset:
203: args: -n 120 -fn_scale 0.6,1.5
204: filter: grep -v "computing matrix functions"
205: output_file: output/test3_4.out
206: test:
207: suffix: 4
208: args: -fn_method {{0 1}}
209: requires: !single
210: test:
211: suffix: 4_subdiagonalpade
212: args: -fn_method {{2 3}}
213: requires: c99_complex !single
215: test:
216: suffix: 5
217: args: -fn_scale 30 -fn_method {{2 3}}
218: filter: grep -v "computing matrix functions"
219: requires: c99_complex !single
220: output_file: output/test3_5.out
222: test:
223: suffix: 6
224: args: -fn_scale 1e-9 -fn_method {{2 3}}
225: filter: grep -v "computing matrix functions"
226: requires: c99_complex !single
227: output_file: output/test3_6.out
229: TEST*/