Line data Source code
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 inverse square root.\n\n";
12 :
13 : #include <slepcfn.h>
14 :
15 : /*
16 : Compute matrix inverse square root B = inv(sqrtm(A))
17 : Check result as norm(B*B*A-I)
18 : */
19 24 : PetscErrorCode TestMatInvSqrt(FN fn,Mat A,PetscViewer viewer,PetscBool verbose,PetscBool inplace)
20 : {
21 24 : PetscScalar tau,eta;
22 24 : PetscReal nrm;
23 24 : PetscBool set,flg;
24 24 : PetscInt n;
25 24 : Mat S,R,Acopy;
26 24 : Vec v,f0;
27 :
28 24 : PetscFunctionBeginUser;
29 24 : PetscCall(MatGetSize(A,&n,NULL));
30 24 : PetscCall(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&S));
31 24 : PetscCall(PetscObjectSetName((PetscObject)S,"S"));
32 24 : PetscCall(FNGetScale(fn,&tau,&eta));
33 : /* compute inverse square root */
34 24 : if (inplace) {
35 12 : PetscCall(MatCopy(A,S,SAME_NONZERO_PATTERN));
36 12 : PetscCall(MatIsHermitianKnown(A,&set,&flg));
37 12 : if (set && flg) PetscCall(MatSetOption(S,MAT_HERMITIAN,PETSC_TRUE));
38 12 : PetscCall(FNEvaluateFunctionMat(fn,S,NULL));
39 : } else {
40 12 : PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&Acopy));
41 12 : PetscCall(FNEvaluateFunctionMat(fn,A,S));
42 : /* check that A has not been modified */
43 12 : PetscCall(MatAXPY(Acopy,-1.0,A,SAME_NONZERO_PATTERN));
44 12 : PetscCall(MatNorm(Acopy,NORM_1,&nrm));
45 12 : if (nrm>100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: the input matrix has changed by %g\n",(double)nrm));
46 12 : PetscCall(MatDestroy(&Acopy));
47 : }
48 24 : if (verbose) {
49 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix A - - - - - - - -\n"));
50 0 : PetscCall(MatView(A,viewer));
51 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed inv(sqrtm(A)) - - - - - - -\n"));
52 0 : PetscCall(MatView(S,viewer));
53 : }
54 : /* check error ||S*S*A-I||_F */
55 24 : PetscCall(MatMatMult(S,S,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&R));
56 24 : if (eta!=1.0) PetscCall(MatScale(R,1.0/(eta*eta)));
57 24 : PetscCall(MatCreateVecs(A,&v,&f0));
58 24 : PetscCall(MatGetColumnVector(S,f0,0));
59 24 : PetscCall(MatCopy(R,S,SAME_NONZERO_PATTERN));
60 24 : PetscCall(MatDestroy(&R));
61 24 : if (tau!=1.0) PetscCall(MatScale(S,tau));
62 24 : PetscCall(MatMatMult(S,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&R));
63 24 : PetscCall(MatShift(R,-1.0));
64 24 : PetscCall(MatNorm(R,NORM_FROBENIUS,&nrm));
65 24 : PetscCall(MatDestroy(&R));
66 24 : if (nrm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||S*S*A-I||_F < 100*eps\n"));
67 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||S*S*A-I||_F = %g\n",(double)nrm));
68 : /* check FNEvaluateFunctionMatVec() */
69 24 : PetscCall(FNEvaluateFunctionMatVec(fn,A,v));
70 24 : PetscCall(VecAXPY(v,-1.0,f0));
71 24 : PetscCall(VecNorm(v,NORM_2,&nrm));
72 24 : 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));
73 24 : PetscCall(MatDestroy(&S));
74 24 : PetscCall(VecDestroy(&v));
75 24 : PetscCall(VecDestroy(&f0));
76 24 : PetscFunctionReturn(PETSC_SUCCESS);
77 : }
78 :
79 8 : int main(int argc,char **argv)
80 : {
81 8 : FN fn;
82 8 : Mat A=NULL;
83 8 : PetscInt i,j,n=10;
84 8 : PetscScalar x,y,yp,*As;
85 8 : PetscViewer viewer;
86 8 : PetscBool verbose,inplace,matcuda;
87 8 : PetscRandom myrand;
88 8 : PetscReal v;
89 8 : char strx[50],str[50];
90 :
91 8 : PetscFunctionBeginUser;
92 8 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
93 8 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
94 8 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
95 8 : PetscCall(PetscOptionsHasName(NULL,NULL,"-inplace",&inplace));
96 8 : PetscCall(PetscOptionsHasName(NULL,NULL,"-matcuda",&matcuda));
97 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix inverse square root, n=%" PetscInt_FMT ".\n",n));
98 :
99 : /* Create function object */
100 8 : PetscCall(FNCreate(PETSC_COMM_WORLD,&fn));
101 8 : PetscCall(FNSetType(fn,FNINVSQRT));
102 8 : PetscCall(FNSetFromOptions(fn));
103 :
104 : /* Set up viewer */
105 8 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
106 8 : PetscCall(FNView(fn,viewer));
107 8 : if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
108 :
109 : /* Scalar evaluation */
110 8 : x = 2.2;
111 8 : PetscCall(SlepcSNPrintfScalar(strx,sizeof(strx),x,PETSC_FALSE));
112 8 : PetscCall(FNEvaluateFunction(fn,x,&y));
113 8 : PetscCall(FNEvaluateDerivative(fn,x,&yp));
114 8 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),y,PETSC_FALSE));
115 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," f(%s)=%s\n",strx,str));
116 8 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),yp,PETSC_FALSE));
117 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," f'(%s)=%s\n",strx,str));
118 :
119 : /* Create matrix */
120 8 : if (matcuda) {
121 : #if defined(PETSC_HAVE_CUDA)
122 : PetscCall(MatCreateSeqDenseCUDA(PETSC_COMM_SELF,n,n,NULL,&A));
123 : #endif
124 8 : } else PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A));
125 8 : PetscCall(PetscObjectSetName((PetscObject)A,"A"));
126 :
127 : /* Compute square root of a symmetric matrix A */
128 8 : PetscCall(MatDenseGetArray(A,&As));
129 88 : for (i=0;i<n;i++) As[i+i*n]=2.5;
130 24 : for (j=1;j<3;j++) {
131 152 : for (i=0;i<n-j;i++) { As[i+(i+j)*n]=1.0; As[(i+j)+i*n]=1.0; }
132 : }
133 8 : PetscCall(MatDenseRestoreArray(A,&As));
134 8 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
135 8 : PetscCall(TestMatInvSqrt(fn,A,viewer,verbose,inplace));
136 :
137 : /* Repeat with upper triangular A */
138 8 : PetscCall(MatDenseGetArray(A,&As));
139 24 : for (j=1;j<3;j++) {
140 152 : for (i=0;i<n-j;i++) As[(i+j)+i*n]=0.0;
141 : }
142 8 : PetscCall(MatDenseRestoreArray(A,&As));
143 8 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
144 8 : PetscCall(TestMatInvSqrt(fn,A,viewer,verbose,inplace));
145 :
146 : /* Repeat with non-symmetic A */
147 8 : PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&myrand));
148 8 : PetscCall(PetscRandomSetFromOptions(myrand));
149 8 : PetscCall(PetscRandomSetInterval(myrand,0.0,1.0));
150 8 : PetscCall(MatDenseGetArray(A,&As));
151 24 : for (j=1;j<3;j++) {
152 152 : for (i=0;i<n-j;i++) {
153 136 : PetscCall(PetscRandomGetValueReal(myrand,&v));
154 136 : As[(i+j)+i*n]=v;
155 : }
156 : }
157 8 : PetscCall(MatDenseRestoreArray(A,&As));
158 8 : PetscCall(PetscRandomDestroy(&myrand));
159 8 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_FALSE));
160 8 : PetscCall(TestMatInvSqrt(fn,A,viewer,verbose,inplace));
161 :
162 8 : PetscCall(MatDestroy(&A));
163 8 : PetscCall(FNDestroy(&fn));
164 8 : PetscCall(SlepcFinalize());
165 : return 0;
166 : }
167 :
168 : /*TEST
169 :
170 : testset:
171 : args: -fn_scale 0.9,0.5 -n 10
172 : filter: grep -v "computing matrix functions"
173 : requires: !__float128
174 : output_file: output/test8_1.out
175 : test:
176 : suffix: 1
177 : args: -fn_method {{0 1 2 3}}
178 : test:
179 : suffix: 1_cuda
180 : args: -fn_method 2 -matcuda
181 : requires: cuda
182 : test:
183 : suffix: 1_magma
184 : args: -fn_method {{1 3}} -matcuda
185 : requires: cuda magma
186 : test:
187 : suffix: 2
188 : args: -inplace -fn_method {{0 1 2 3}}
189 : test:
190 : suffix: 2_cuda
191 : args: -inplace -fn_method 2 -matcuda
192 : requires: cuda
193 : test:
194 : suffix: 2_magma
195 : args: -inplace -fn_method {{1 3}} -matcuda
196 : requires: cuda magma
197 :
198 : TEST*/
|
