Line | Branch | Exec | Source |
---|---|---|---|
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 logarithm.\n\n"; | ||
12 | |||
13 | #include <slepcfn.h> | ||
14 | |||
15 | /* | ||
16 | Compute matrix logarithm B = logm(A) | ||
17 | */ | ||
18 | 25 | PetscErrorCode TestMatLog(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace) | |
19 | { | ||
20 | 25 | PetscBool set,flg; | |
21 | 25 | PetscScalar tau,eta; | |
22 | 25 | PetscInt n; | |
23 | 25 | Mat F,R; | |
24 | 25 | Vec v,f0; | |
25 | 25 | FN fnexp; | |
26 | 25 | PetscReal nrm; | |
27 | |||
28 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
25 | PetscFunctionBeginUser; |
29 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatGetSize(A,&n,NULL)); |
30 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&F)); |
31 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscObjectSetName((PetscObject)F,"F")); |
32 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&R)); |
33 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscObjectSetName((PetscObject)R,"R")); |
34 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNGetScale(fn,&tau,&eta)); |
35 | /* compute matrix logarithm */ | ||
36 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
25 | if (inplace) { |
37 | ✗ | PetscCall(MatCopy(A,F,SAME_NONZERO_PATTERN)); | |
38 | ✗ | PetscCall(MatIsHermitianKnown(A,&set,&flg)); | |
39 | ✗ | if (set && flg) PetscCall(MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE)); | |
40 | ✗ | PetscCall(FNEvaluateFunctionMat(fn,F,NULL)); | |
41 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | } else PetscCall(FNEvaluateFunctionMat(fn,A,F)); |
42 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
25 | if (verbose) { |
43 | ✗ | PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n")); | |
44 | ✗ | PetscCall(MatView(A,viewer)); | |
45 | ✗ | PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed logm(A) - - - - - - -\n")); | |
46 | ✗ | PetscCall(MatView(F,viewer)); | |
47 | } | ||
48 | /* check error ||expm(F)-A||_F */ | ||
49 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNCreate(PETSC_COMM_WORLD,&fnexp)); |
50 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNSetType(fnexp,FNEXP)); |
51 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatCopy(F,R,SAME_NONZERO_PATTERN)); |
52 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
25 | if (eta!=1.0) PetscCall(MatScale(R,1.0/eta)); |
53 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNEvaluateFunctionMat(fnexp,R,NULL)); |
54 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNDestroy(&fnexp)); |
55 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatAXPY(R,-tau,A,SAME_NONZERO_PATTERN)); |
56 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatNorm(R,NORM_FROBENIUS,&nrm)); |
57 |
5/8✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
25 | if (nrm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||expm(F)-A||_F < 100*eps\n")); |
58 | ✗ | else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||expm(F)-A||_F = %g\n",(double)nrm)); | |
59 | /* check FNEvaluateFunctionMatVec() */ | ||
60 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatCreateVecs(A,&v,&f0)); |
61 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatGetColumnVector(F,f0,0)); |
62 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNEvaluateFunctionMatVec(fn,A,v)); |
63 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(VecAXPY(v,-1.0,f0)); |
64 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(VecNorm(v,NORM_2,&nrm)); |
65 |
1/8✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
25 | 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)); |
66 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatDestroy(&F)); |
67 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatDestroy(&R)); |
68 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(VecDestroy(&v)); |
69 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(VecDestroy(&f0)); |
70 |
5/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
5 | PetscFunctionReturn(PETSC_SUCCESS); |
71 | } | ||
72 | |||
73 | 25 | int main(int argc,char **argv) | |
74 | { | ||
75 | 25 | FN fn; | |
76 | 25 | Mat A; | |
77 | 25 | PetscInt i,j,n=10; | |
78 | 25 | PetscScalar *As; | |
79 | 25 | PetscViewer viewer; | |
80 | 25 | PetscBool verbose,inplace,random,triang; | |
81 | |||
82 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
25 | PetscFunctionBeginUser; |
83 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(SlepcInitialize(&argc,&argv,NULL,help)); |
84 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); |
85 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose)); |
86 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscOptionsHasName(NULL,NULL,"-inplace",&inplace)); |
87 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscOptionsHasName(NULL,NULL,"-random",&random)); |
88 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscOptionsHasName(NULL,NULL,"-triang",&triang)); |
89 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix logarithm, n=%" PetscInt_FMT ".\n",n)); |
90 | |||
91 | /* Create logarithm function object */ | ||
92 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNCreate(PETSC_COMM_WORLD,&fn)); |
93 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNSetType(fn,FNLOG)); |
94 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNSetFromOptions(fn)); |
95 | |||
96 | /* Set up viewer */ | ||
97 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer)); |
98 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNView(fn,viewer)); |
99 |
1/8✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
25 | if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB)); |
100 | |||
101 | /* Create matrices */ | ||
102 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A)); |
103 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(PetscObjectSetName((PetscObject)A,"A")); |
104 | |||
105 |
6/8✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
|
25 | if (random) PetscCall(MatSetRandom(A,NULL)); |
106 | else { | ||
107 | /* Fill A with a non-symmetric Toeplitz matrix */ | ||
108 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
20 | PetscCall(MatDenseGetArray(A,&As)); |
109 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1520 | for (i=0;i<n;i++) As[i+i*n]=2.0; |
110 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
60 | for (j=1;j<3;j++) { |
111 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2980 | for (i=0;i<n-j;i++) { |
112 | 2940 | As[i+(i+j)*n]=1.0; | |
113 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
2940 | if (!triang) As[(i+j)+i*n]=-1.0; |
114 | } | ||
115 | } | ||
116 | 20 | As[(n-1)*n] = -5.0; | |
117 | 20 | As[0] = 2.01; | |
118 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
20 | PetscCall(MatDenseRestoreArray(A,&As)); |
119 | } | ||
120 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(TestMatLog(fn,A,viewer,verbose,inplace)); |
121 | |||
122 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(MatDestroy(&A)); |
123 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
25 | PetscCall(FNDestroy(&fn)); |
124 |
3/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
25 | PetscCall(SlepcFinalize()); |
125 | return 0; | ||
126 | } | ||
127 | |||
128 | /*TEST | ||
129 | |||
130 | testset: | ||
131 | filter: grep -v "computing matrix functions" | ||
132 | output_file: output/test13_1.out | ||
133 | test: | ||
134 | suffix: 1 | ||
135 | args: -fn_scale .04,2 -n 75 | ||
136 | requires: c99_complex !__float128 | ||
137 | test: | ||
138 | suffix: 1_triang | ||
139 | args: -fn_scale .04,2 -n 75 -triang | ||
140 | requires: c99_complex !__float128 | ||
141 | test: | ||
142 | suffix: 1_random | ||
143 | args: -fn_scale .02,2 -n 75 -random | ||
144 | requires: complex !__float128 | ||
145 | filter_output: sed -e 's/04/02/' | ||
146 | |||
147 | TEST*/ | ||
148 |