Actual source code: test6.c
slepc-3.21.1 2024-04-26
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: */
10: /*
11: Define the function
13: f(x) = (1-x^2) exp(-x/(1+x^2))
15: with the following tree:
17: f(x) f(x) (combined by product)
18: / \ g(x) = 1-x^2 (polynomial)
19: g(x) h(x) h(x) (combined by composition)
20: / \ r(x) = -x/(1+x^2) (rational)
21: r(x) e(x) e(x) = exp(x) (exponential)
22: */
24: static char help[] = "Test combined function.\n\n";
26: #include <slepcfn.h>
28: /*
29: Compute matrix function B = (I-A^2) exp(-(I+A^2)\A)
30: */
31: PetscErrorCode TestMatCombine(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace)
32: {
33: PetscBool set,flg;
34: PetscInt n;
35: Mat F,Acopy;
36: Vec v,f0;
37: PetscReal nrm;
39: PetscFunctionBeginUser;
40: PetscCall(MatGetSize(A,&n,NULL));
41: PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&F));
42: PetscCall(PetscObjectSetName((PetscObject)F,"F"));
43: /* compute matrix function */
44: if (inplace) {
45: PetscCall(MatCopy(A,F,SAME_NONZERO_PATTERN));
46: PetscCall(MatIsHermitianKnown(A,&set,&flg));
47: if (set && flg) PetscCall(MatSetOption(F,MAT_HERMITIAN,PETSC_TRUE));
48: PetscCall(FNEvaluateFunctionMat(fn,F,NULL));
49: } else {
50: PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Acopy));
51: PetscCall(FNEvaluateFunctionMat(fn,A,F));
52: /* check that A has not been modified */
53: PetscCall(MatAXPY(Acopy,-1.0,A,SAME_NONZERO_PATTERN));
54: PetscCall(MatNorm(Acopy,NORM_1,&nrm));
55: if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the input matrix has changed by %g\n",(double)nrm));
56: PetscCall(MatDestroy(&Acopy));
57: }
58: if (verbose) {
59: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n"));
60: PetscCall(MatView(A,viewer));
61: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed f(A) - - - - - - -\n"));
62: PetscCall(MatView(F,viewer));
63: }
64: /* print matrix norm for checking */
65: PetscCall(MatNorm(F,NORM_1,&nrm));
66: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"The 1-norm of f(A) is %6.3f\n",(double)nrm));
67: /* check FNEvaluateFunctionMatVec() */
68: PetscCall(MatCreateVecs(A,&v,&f0));
69: PetscCall(MatGetColumnVector(F,f0,0));
70: PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
71: PetscCall(VecAXPY(v,-1.0,f0));
72: PetscCall(VecNorm(v,NORM_2,&nrm));
73: 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));
74: PetscCall(MatDestroy(&F));
75: PetscCall(VecDestroy(&v));
76: PetscCall(VecDestroy(&f0));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: int main(int argc,char **argv)
81: {
82: FN f,g,h,e,r,fcopy;
83: Mat A=NULL;
84: PetscInt i,j,n=10,np,nq;
85: PetscScalar x,y,yp,*As,p[10],q[10];
86: char strx[50],str[50];
87: PetscViewer viewer;
88: PetscBool verbose,inplace,matcuda;
90: PetscFunctionBeginUser;
91: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
92: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
93: PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
94: PetscCall(PetscOptionsHasName(NULL,NULL,"-inplace",&inplace));
95: PetscCall(PetscOptionsHasName(NULL,NULL,"-matcuda",&matcuda));
96: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Combined function, n=%" PetscInt_FMT ".\n",n));
98: /* Create function */
100: /* e(x) = exp(x) */
101: PetscCall(FNCreate(PETSC_COMM_WORLD,&e));
102: PetscCall(FNSetType(e,FNEXP));
103: PetscCall(FNSetFromOptions(e));
104: /* r(x) = x/(1+x^2) */
105: PetscCall(FNCreate(PETSC_COMM_WORLD,&r));
106: PetscCall(FNSetType(r,FNRATIONAL));
107: PetscCall(FNSetFromOptions(r));
108: np = 2; nq = 3;
109: p[0] = -1.0; p[1] = 0.0;
110: q[0] = 1.0; q[1] = 0.0; q[2] = 1.0;
111: PetscCall(FNRationalSetNumerator(r,np,p));
112: PetscCall(FNRationalSetDenominator(r,nq,q));
113: /* h(x) */
114: PetscCall(FNCreate(PETSC_COMM_WORLD,&h));
115: PetscCall(FNSetType(h,FNCOMBINE));
116: PetscCall(FNSetFromOptions(h));
117: PetscCall(FNCombineSetChildren(h,FN_COMBINE_COMPOSE,r,e));
118: /* g(x) = 1-x^2 */
119: PetscCall(FNCreate(PETSC_COMM_WORLD,&g));
120: PetscCall(FNSetType(g,FNRATIONAL));
121: PetscCall(FNSetFromOptions(g));
122: np = 3;
123: p[0] = -1.0; p[1] = 0.0; p[2] = 1.0;
124: PetscCall(FNRationalSetNumerator(g,np,p));
125: /* f(x) */
126: PetscCall(FNCreate(PETSC_COMM_WORLD,&f));
127: PetscCall(FNSetType(f,FNCOMBINE));
128: PetscCall(FNSetFromOptions(f));
129: PetscCall(FNCombineSetChildren(f,FN_COMBINE_MULTIPLY,g,h));
131: /* Set up viewer */
132: PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
133: PetscCall(FNView(f,viewer));
134: if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
136: /* Scalar evaluation */
137: x = 2.2;
138: PetscCall(SlepcSNPrintfScalar(strx,sizeof(strx),x,PETSC_FALSE));
139: PetscCall(FNEvaluateFunction(f,x,&y));
140: PetscCall(FNEvaluateDerivative(f,x,&yp));
141: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),y,PETSC_FALSE));
142: PetscCall(PetscPrintf(PETSC_COMM_WORLD," f(%s)=%s\n",strx,str));
143: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),yp,PETSC_FALSE));
144: PetscCall(PetscPrintf(PETSC_COMM_WORLD," f'(%s)=%s\n",strx,str));
146: /* Test duplication */
147: PetscCall(FNDuplicate(f,PetscObjectComm((PetscObject)f),&fcopy));
149: /* Create matrices */
150: if (matcuda) {
151: #if defined(PETSC_HAVE_CUDA)
152: PetscCall(MatCreateSeqDenseCUDA(PETSC_COMM_SELF,n,n,NULL,&A));
153: #endif
154: } else PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
155: PetscCall(PetscObjectSetName((PetscObject)A,"A"));
157: /* Fill A with a symmetric Toeplitz matrix */
158: PetscCall(MatDenseGetArray(A,&As));
159: for (i=0;i<n;i++) As[i+i*n]=2.0;
160: for (j=1;j<3;j++) {
161: for (i=0;i<n-j;i++) { As[i+(i+j)*n]=1.0; As[(i+j)+i*n]=1.0; }
162: }
163: PetscCall(MatDenseRestoreArray(A,&As));
164: PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
165: PetscCall(TestMatCombine(fcopy,A,viewer,verbose,inplace));
167: /* Repeat with same matrix as non-symmetric */
168: PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
169: PetscCall(TestMatCombine(fcopy,A,viewer,verbose,inplace));
171: PetscCall(MatDestroy(&A));
172: PetscCall(FNDestroy(&f));
173: PetscCall(FNDestroy(&fcopy));
174: PetscCall(FNDestroy(&g));
175: PetscCall(FNDestroy(&h));
176: PetscCall(FNDestroy(&e));
177: PetscCall(FNDestroy(&r));
178: PetscCall(SlepcFinalize());
179: return 0;
180: }
182: /*TEST
184: testset:
185: output_file: output/test6_1.out
186: test:
187: suffix: 1
188: test:
189: suffix: 1_cuda
190: args: -matcuda
191: requires: cuda
192: test:
193: suffix: 2
194: args: -inplace
195: test:
196: suffix: 2_cuda
197: args: -inplace -matcuda
198: requires: cuda
200: TEST*/