Actual source code: test13.c
slepc-3.21.2 2024-09-25
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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 logarithm.\n\n";
13: #include <slepcfn.h>
15: /*
16: Compute matrix logarithm B = logm(A)
17: */
18: PetscErrorCode TestMatLog(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace)
19: {
20: PetscBool set,flg;
21: PetscScalar tau,eta;
22: PetscInt n;
23: Mat F,R;
24: Vec v,f0;
25: FN fnexp;
26: PetscReal nrm;
28: PetscFunctionBeginUser;
29: PetscCall(MatGetSize(A,&n,NULL));
30: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&F));
31: PetscCall(PetscObjectSetName((PetscObject)F,"F"));
32: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&R));
33: PetscCall(PetscObjectSetName((PetscObject)R,"R"));
34: PetscCall(FNGetScale(fn,&tau,&eta));
35: /* compute matrix logarithm */
36: 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: } else PetscCall(FNEvaluateFunctionMat(fn,A,F));
42: 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: PetscCall(FNCreate(PETSC_COMM_WORLD,&fnexp));
50: PetscCall(FNSetType(fnexp,FNEXP));
51: PetscCall(MatCopy(F,R,SAME_NONZERO_PATTERN));
52: if (eta!=1.0) PetscCall(MatScale(R,1.0/eta));
53: PetscCall(FNEvaluateFunctionMat(fnexp,R,NULL));
54: PetscCall(FNDestroy(&fnexp));
55: PetscCall(MatAXPY(R,-tau,A,SAME_NONZERO_PATTERN));
56: PetscCall(MatNorm(R,NORM_FROBENIUS,&nrm));
57: 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: PetscCall(MatCreateVecs(A,&v,&f0));
61: PetscCall(MatGetColumnVector(F,f0,0));
62: PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
63: PetscCall(VecAXPY(v,-1.0,f0));
64: PetscCall(VecNorm(v,NORM_2,&nrm));
65: 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: PetscCall(MatDestroy(&F));
67: PetscCall(MatDestroy(&R));
68: PetscCall(VecDestroy(&v));
69: PetscCall(VecDestroy(&f0));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: int main(int argc,char **argv)
74: {
75: FN fn;
76: Mat A;
77: PetscInt i,j,n=10;
78: PetscScalar *As;
79: PetscViewer viewer;
80: PetscBool verbose,inplace,random,triang;
82: PetscFunctionBeginUser;
83: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
84: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
85: PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
86: PetscCall(PetscOptionsHasName(NULL,NULL,"-inplace",&inplace));
87: PetscCall(PetscOptionsHasName(NULL,NULL,"-random",&random));
88: PetscCall(PetscOptionsHasName(NULL,NULL,"-triang",&triang));
89: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix logarithm, n=%" PetscInt_FMT ".\n",n));
91: /* Create logarithm function object */
92: PetscCall(FNCreate(PETSC_COMM_WORLD,&fn));
93: PetscCall(FNSetType(fn,FNLOG));
94: PetscCall(FNSetFromOptions(fn));
96: /* Set up viewer */
97: PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
98: PetscCall(FNView(fn,viewer));
99: if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
101: /* Create matrices */
102: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
103: PetscCall(PetscObjectSetName((PetscObject)A,"A"));
105: if (random) PetscCall(MatSetRandom(A,NULL));
106: else {
107: /* Fill A with a non-symmetric Toeplitz matrix */
108: PetscCall(MatDenseGetArray(A,&As));
109: for (i=0;i<n;i++) As[i+i*n]=2.0;
110: for (j=1;j<3;j++) {
111: for (i=0;i<n-j;i++) {
112: As[i+(i+j)*n]=1.0;
113: if (!triang) As[(i+j)+i*n]=-1.0;
114: }
115: }
116: As[(n-1)*n] = -5.0;
117: As[0] = 2.01;
118: PetscCall(MatDenseRestoreArray(A,&As));
119: }
120: PetscCall(TestMatLog(fn,A,viewer,verbose,inplace));
122: PetscCall(MatDestroy(&A));
123: PetscCall(FNDestroy(&fn));
124: PetscCall(SlepcFinalize());
125: return 0;
126: }
128: /*TEST
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/'
147: TEST*/