Actual source code: test3.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2016, 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: */
22: static char help[] = "Test matrix exponential.\n\n";
24: #include <slepcfn.h>
28: /*
29: Compute matrix exponential B = expm(A)
30: */
31: PetscErrorCode TestMatExp(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace)
32: {
34: PetscBool set,flg;
35: PetscInt n;
36: Mat F;
37: Vec v,f0;
38: PetscReal nrm;
41: MatGetSize(A,&n,NULL);
42: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&F);
43: PetscObjectSetName((PetscObject)F,"F");
44: /* compute square root */
45: if (inplace) {
46: MatCopy(A,F,SAME_NONZERO_PATTERN);
47: MatIsHermitianKnown(A,&set,&flg);
48: if (set && flg) { MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE); }
49: FNEvaluateFunctionMat(fn,F,NULL);
50: } else {
51: FNEvaluateFunctionMat(fn,A,F);
52: }
53: if (verbose) {
54: PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n");
55: MatView(A,viewer);
56: PetscPrintf(PETSC_COMM_WORLD,"Computed expm(A) - - - - - - -\n");
57: MatView(F,viewer);
58: }
59: /* print matrix norm for checking */
60: MatNorm(F,NORM_1,&nrm);
61: PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %g\n",(double)nrm);
62: /* check FNEvaluateFunctionMatVec() */
63: MatCreateVecs(A,&v,&f0);
64: MatGetColumnVector(F,f0,0);
65: FNEvaluateFunctionMatVec(fn,A,v);
66: VecAXPY(v,-1.0,f0);
67: VecNorm(v,NORM_2,&nrm);
68: if (nrm>100*PETSC_MACHINE_EPSILON) {
69: PetscPrintf(PETSC_COMM_WORLD,"Warning: the norm of f(A)*e_1-v is %g\n",(double)nrm);
70: }
71: MatDestroy(&F);
72: VecDestroy(&v);
73: VecDestroy(&f0);
74: return(0);
75: }
79: int main(int argc,char **argv)
80: {
82: FN fn;
83: Mat A;
84: PetscInt i,j,n=10;
85: PetscScalar *As,tau=1.0,eta=1.0;
86: PetscViewer viewer;
87: PetscBool verbose,inplace;
89: SlepcInitialize(&argc,&argv,(char*)0,help);
90: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
91: PetscOptionsGetScalar(NULL,NULL,"-tau",&tau,NULL);
92: PetscOptionsGetScalar(NULL,NULL,"-eta",&eta,NULL);
93: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
94: PetscOptionsHasName(NULL,NULL,"-inplace",&inplace);
95: PetscPrintf(PETSC_COMM_WORLD,"Matrix exponential, n=%D.\n",n);
97: /* Create exponential function eta*exp(tau*x) */
98: FNCreate(PETSC_COMM_WORLD,&fn);
99: FNSetType(fn,FNEXP);
100: FNSetScale(fn,tau,eta);
102: /* Set up viewer */
103: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
104: FNView(fn,viewer);
105: if (verbose) {
106: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
107: }
109: /* Create matrices */
110: MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
111: PetscObjectSetName((PetscObject)A,"A");
113: /* Fill A with a symmetric Toeplitz matrix */
114: MatDenseGetArray(A,&As);
115: for (i=0;i<n;i++) As[i+i*n]=2.0;
116: for (j=1;j<3;j++) {
117: for (i=0;i<n-j;i++) { As[i+(i+j)*n]=1.0; As[(i+j)+i*n]=1.0; }
118: }
119: MatDenseRestoreArray(A,&As);
120: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
121: TestMatExp(fn,A,viewer,verbose,inplace);
123: /* Repeat with non-symmetric A */
124: MatDenseGetArray(A,&As);
125: for (j=1;j<3;j++) {
126: for (i=0;i<n-j;i++) { As[(i+j)+i*n]=-1.0; }
127: }
128: MatDenseRestoreArray(A,&As);
129: MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE);
130: TestMatExp(fn,A,viewer,verbose,inplace);
132: MatDestroy(&A);
133: FNDestroy(&fn);
134: SlepcFinalize();
135: return ierr;
136: }