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 BV bi-orthogonalization functions.\n\n";
12 :
13 : #include <slepcbv.h>
14 :
15 8 : int main(int argc,char **argv)
16 : {
17 8 : BV X,Y;
18 8 : Mat M;
19 8 : Vec v,t;
20 8 : PetscInt i,j,n=20,k=7;
21 8 : PetscViewer view;
22 8 : PetscBool verbose;
23 8 : PetscReal norm,condn=1.0;
24 8 : PetscScalar alpha;
25 :
26 8 : PetscFunctionBeginUser;
27 8 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
28 8 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
29 8 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
30 8 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-condn",&condn,NULL));
31 8 : PetscCheck(condn>=1.0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The condition number must be > 1");
32 8 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
33 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV bi-orthogonalization with %" PetscInt_FMT " columns of length %" PetscInt_FMT ".\n",k,n));
34 8 : if (condn>1.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," - Using random BVs with condition number = %g\n",(double)condn));
35 :
36 : /* Create template vector */
37 8 : PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
38 8 : PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
39 8 : PetscCall(VecSetFromOptions(t));
40 :
41 : /* Create BV object X */
42 8 : PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
43 8 : PetscCall(PetscObjectSetName((PetscObject)X,"X"));
44 8 : PetscCall(BVSetSizesFromVec(X,t,k));
45 8 : PetscCall(BVSetFromOptions(X));
46 :
47 : /* Set up viewer */
48 8 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
49 8 : if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
50 :
51 : /* Fill X entries */
52 8 : if (condn==1.0) {
53 64 : for (j=0;j<k;j++) {
54 56 : PetscCall(BVGetColumn(X,j,&v));
55 56 : PetscCall(VecSet(v,0.0));
56 672 : for (i=0;i<=n/2;i++) {
57 616 : if (i+j<n) {
58 616 : alpha = (3.0*i+j-2)/(2*(i+j+1));
59 : #if defined(PETSC_USE_COMPLEX)
60 : alpha += (1.2*i+j-2)/(0.1*(i+j+1))*PETSC_i;
61 : #endif
62 616 : PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
63 : }
64 : }
65 56 : PetscCall(VecAssemblyBegin(v));
66 56 : PetscCall(VecAssemblyEnd(v));
67 56 : PetscCall(BVRestoreColumn(X,j,&v));
68 : }
69 0 : } else PetscCall(BVSetRandomCond(X,condn));
70 8 : if (verbose) PetscCall(BVView(X,view));
71 :
72 : /* Create Y and fill its entries */
73 8 : PetscCall(BVDuplicate(X,&Y));
74 8 : PetscCall(PetscObjectSetName((PetscObject)Y,"Y"));
75 8 : if (condn==1.0) {
76 64 : for (j=0;j<k;j++) {
77 56 : PetscCall(BVGetColumn(Y,j,&v));
78 56 : PetscCall(VecSet(v,0.0));
79 384 : for (i=PetscMin(n,2+(2*j)%6);i<PetscMin(n,6+(3*j)%9);i++) {
80 272 : if (i%5 && i!=j) {
81 200 : alpha = (1.5*i+j)/(2.2*(i-j));
82 272 : PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
83 : }
84 : }
85 56 : PetscCall(VecAssemblyBegin(v));
86 56 : PetscCall(VecAssemblyEnd(v));
87 56 : PetscCall(BVRestoreColumn(Y,j,&v));
88 : }
89 0 : } else PetscCall(BVSetRandomCond(Y,condn));
90 8 : if (verbose) PetscCall(BVView(Y,view));
91 :
92 : /* Test BVBiorthonormalizeColumn */
93 64 : for (j=0;j<k;j++) PetscCall(BVBiorthonormalizeColumn(X,Y,j,NULL));
94 8 : if (verbose) {
95 0 : PetscCall(BVView(X,view));
96 0 : PetscCall(BVView(Y,view));
97 : }
98 :
99 : /* Check orthogonality */
100 8 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M));
101 8 : PetscCall(PetscObjectSetName((PetscObject)M,"M"));
102 8 : PetscCall(BVDot(X,Y,M));
103 8 : if (verbose) PetscCall(MatView(M,view));
104 8 : PetscCall(MatShift(M,-1.0));
105 8 : PetscCall(MatNorm(M,NORM_1,&norm));
106 8 : if (norm<200*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of bi-orthogonality < 200*eps\n"));
107 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of bi-orthogonality: %g\n",(double)norm));
108 :
109 8 : PetscCall(MatDestroy(&M));
110 8 : PetscCall(BVDestroy(&X));
111 8 : PetscCall(BVDestroy(&Y));
112 8 : PetscCall(VecDestroy(&t));
113 8 : PetscCall(SlepcFinalize());
114 : return 0;
115 : }
116 :
117 : /*TEST
118 :
119 : testset:
120 : output_file: output/test17_1.out
121 : test:
122 : suffix: 1
123 : args: -bv_type {{vecs contiguous svec mat}} -bv_orthog_type cgs
124 : requires: double
125 : test:
126 : suffix: 1_cuda
127 : args: -bv_type {{svec mat}} -vec_type cuda -bv_orthog_type cgs
128 : requires: cuda
129 : test:
130 : suffix: 1_hip
131 : args: -bv_type {{svec mat}} -vec_type hip -bv_orthog_type cgs
132 : requires: hip
133 : test:
134 : suffix: 2
135 : args: -bv_type {{vecs contiguous svec mat}} -bv_orthog_type mgs
136 : test:
137 : suffix: 2_cuda
138 : args: -bv_type {{svec mat}} -vec_type cuda -bv_orthog_type mgs
139 : requires: cuda
140 : test:
141 : suffix: 2_hip
142 : args: -bv_type {{svec mat}} -vec_type hip -bv_orthog_type mgs
143 : requires: hip
144 :
145 : TEST*/
|
