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 matrix projection.\n\n";
12 :
13 : #include <slepcbv.h>
14 :
15 15 : int main(int argc,char **argv)
16 : {
17 15 : Vec t,v;
18 15 : Mat B,G,H0,H1;
19 15 : BV X,Y,Z;
20 15 : PetscInt i,j,n=20,kx=6,lx=3,ky=5,ly=2,Istart,Iend,col[5];
21 15 : PetscScalar alpha,value[] = { -1, 1, 1, 1, 1 };
22 15 : PetscViewer view;
23 15 : PetscReal norm;
24 15 : PetscBool verbose;
25 :
26 15 : PetscFunctionBeginUser;
27 15 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
28 15 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
29 15 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-kx",&kx,NULL));
30 15 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-lx",&lx,NULL));
31 15 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-ky",&ky,NULL));
32 15 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-ly",&ly,NULL));
33 15 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
34 15 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV projection (n=%" PetscInt_FMT ").\n",n));
35 15 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"X has %" PetscInt_FMT " active columns (%" PetscInt_FMT " leading columns).\n",kx,lx));
36 15 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Y has %" PetscInt_FMT " active columns (%" PetscInt_FMT " leading columns).\n",ky,ly));
37 :
38 : /* Set up viewer */
39 15 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
40 15 : if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
41 :
42 : /* Create non-symmetric matrix G (Toeplitz) */
43 15 : PetscCall(MatCreate(PETSC_COMM_WORLD,&G));
44 15 : PetscCall(MatSetSizes(G,PETSC_DECIDE,PETSC_DECIDE,n,n));
45 15 : PetscCall(MatSetFromOptions(G));
46 15 : PetscCall(PetscObjectSetName((PetscObject)G,"G"));
47 :
48 15 : PetscCall(MatGetOwnershipRange(G,&Istart,&Iend));
49 215 : for (i=Istart;i<Iend;i++) {
50 200 : col[0]=i-1; col[1]=i; col[2]=i+1; col[3]=i+2; col[4]=i+3;
51 200 : if (i==0) PetscCall(MatSetValues(G,1,&i,PetscMin(4,n-i),col+1,value+1,INSERT_VALUES));
52 200 : else PetscCall(MatSetValues(G,1,&i,PetscMin(5,n-i+1),col,value,INSERT_VALUES));
53 : }
54 15 : PetscCall(MatAssemblyBegin(G,MAT_FINAL_ASSEMBLY));
55 15 : PetscCall(MatAssemblyEnd(G,MAT_FINAL_ASSEMBLY));
56 15 : if (verbose) PetscCall(MatView(G,view));
57 :
58 : /* Create symmetric matrix B (1-D Laplacian) */
59 15 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
60 15 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
61 15 : PetscCall(MatSetFromOptions(B));
62 15 : PetscCall(PetscObjectSetName((PetscObject)B,"B"));
63 :
64 15 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
65 215 : for (i=Istart;i<Iend;i++) {
66 200 : if (i>0) PetscCall(MatSetValue(B,i,i-1,-1.0,INSERT_VALUES));
67 200 : if (i<n-1) PetscCall(MatSetValue(B,i,i+1,-1.0,INSERT_VALUES));
68 200 : PetscCall(MatSetValue(B,i,i,2.0,INSERT_VALUES));
69 : }
70 15 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
71 15 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
72 15 : PetscCall(MatCreateVecs(B,&t,NULL));
73 15 : if (verbose) PetscCall(MatView(B,view));
74 :
75 : /* Create BV object X */
76 15 : PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
77 15 : PetscCall(PetscObjectSetName((PetscObject)X,"X"));
78 15 : PetscCall(BVSetSizesFromVec(X,t,kx+2)); /* two extra columns to test active columns */
79 15 : PetscCall(BVSetFromOptions(X));
80 :
81 : /* Fill X entries */
82 135 : for (j=0;j<kx+2;j++) {
83 120 : PetscCall(BVGetColumn(X,j,&v));
84 120 : PetscCall(VecSet(v,0.0));
85 600 : for (i=0;i<4;i++) {
86 480 : if (i+j<n) {
87 : #if defined(PETSC_USE_COMPLEX)
88 480 : alpha = PetscCMPLX((PetscReal)(3*i+j-2),(PetscReal)(2*i));
89 : #else
90 : alpha = (PetscReal)(3*i+j-2);
91 : #endif
92 480 : PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
93 : }
94 : }
95 120 : PetscCall(VecAssemblyBegin(v));
96 120 : PetscCall(VecAssemblyEnd(v));
97 120 : PetscCall(BVRestoreColumn(X,j,&v));
98 : }
99 15 : if (verbose) PetscCall(BVView(X,view));
100 :
101 : /* Duplicate BV object and store Z=G*X */
102 15 : PetscCall(BVDuplicate(X,&Z));
103 15 : PetscCall(PetscObjectSetName((PetscObject)Z,"Z"));
104 15 : PetscCall(BVSetActiveColumns(X,0,kx));
105 15 : PetscCall(BVSetActiveColumns(Z,0,kx));
106 15 : PetscCall(BVMatMult(X,G,Z));
107 15 : PetscCall(BVSetActiveColumns(X,lx,kx));
108 15 : PetscCall(BVSetActiveColumns(Z,lx,kx));
109 :
110 : /* Create BV object Y */
111 15 : PetscCall(BVCreate(PETSC_COMM_WORLD,&Y));
112 15 : PetscCall(PetscObjectSetName((PetscObject)Y,"Y"));
113 15 : PetscCall(BVSetSizesFromVec(Y,t,ky+1));
114 15 : PetscCall(BVSetFromOptions(Y));
115 15 : PetscCall(BVSetActiveColumns(Y,ly,ky));
116 :
117 : /* Fill Y entries */
118 105 : for (j=0;j<ky+1;j++) {
119 90 : PetscCall(BVGetColumn(Y,j,&v));
120 : #if defined(PETSC_USE_COMPLEX)
121 90 : alpha = PetscCMPLX((PetscReal)(j+1)/4.0,-(PetscReal)j);
122 : #else
123 : alpha = (PetscReal)(j+1)/4.0;
124 : #endif
125 90 : PetscCall(VecSet(v,(PetscScalar)(j+1)/4.0));
126 90 : PetscCall(BVRestoreColumn(Y,j,&v));
127 : }
128 15 : if (verbose) PetscCall(BVView(Y,view));
129 :
130 : /* Test BVMatProject for non-symmetric matrix G */
131 15 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H0));
132 15 : PetscCall(PetscObjectSetName((PetscObject)H0,"H0"));
133 15 : PetscCall(BVMatProject(X,G,Y,H0));
134 15 : if (verbose) PetscCall(MatView(H0,view));
135 :
136 : /* Test BVMatProject with previously stored G*X */
137 15 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H1));
138 15 : PetscCall(PetscObjectSetName((PetscObject)H1,"H1"));
139 15 : PetscCall(BVMatProject(Z,NULL,Y,H1));
140 15 : if (verbose) PetscCall(MatView(H1,view));
141 :
142 : /* Check that H0 and H1 are equal */
143 15 : PetscCall(MatAXPY(H0,-1.0,H1,SAME_NONZERO_PATTERN));
144 15 : PetscCall(MatNorm(H0,NORM_1,&norm));
145 15 : if (norm<10*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||H0-H1|| < 10*eps\n"));
146 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||H0-H1||=%g\n",(double)norm));
147 15 : PetscCall(MatDestroy(&H0));
148 15 : PetscCall(MatDestroy(&H1));
149 :
150 : /* Test BVMatProject for symmetric matrix B with orthogonal projection */
151 15 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,kx,kx,NULL,&H0));
152 15 : PetscCall(PetscObjectSetName((PetscObject)H0,"H0"));
153 15 : PetscCall(BVMatProject(X,B,X,H0));
154 15 : if (verbose) PetscCall(MatView(H0,view));
155 :
156 : /* Repeat previous test with symmetry flag set */
157 15 : PetscCall(MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE));
158 15 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,kx,kx,NULL,&H1));
159 15 : PetscCall(PetscObjectSetName((PetscObject)H1,"H1"));
160 15 : PetscCall(BVMatProject(X,B,X,H1));
161 15 : if (verbose) PetscCall(MatView(H1,view));
162 :
163 : /* Check that H0 and H1 are equal */
164 15 : PetscCall(MatAXPY(H0,-1.0,H1,SAME_NONZERO_PATTERN));
165 15 : PetscCall(MatNorm(H0,NORM_1,&norm));
166 15 : if (norm<10*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||H0-H1|| < 10*eps\n"));
167 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"||H0-H1||=%g\n",(double)norm));
168 15 : PetscCall(MatDestroy(&H0));
169 15 : PetscCall(MatDestroy(&H1));
170 :
171 15 : PetscCall(BVDestroy(&X));
172 15 : PetscCall(BVDestroy(&Y));
173 15 : PetscCall(BVDestroy(&Z));
174 15 : PetscCall(MatDestroy(&B));
175 15 : PetscCall(MatDestroy(&G));
176 15 : PetscCall(VecDestroy(&t));
177 15 : PetscCall(SlepcFinalize());
178 : return 0;
179 : }
180 :
181 : /*TEST
182 :
183 : testset:
184 : output_file: output/test9_1.out
185 : test:
186 : suffix: 1
187 : args: -bv_type {{vecs contiguous svec mat}shared output}
188 : test:
189 : suffix: 1_svec_vecs
190 : args: -bv_type svec -bv_matmult vecs
191 : test:
192 : suffix: 1_cuda
193 : args: -bv_type {{svec mat}} -mat_type aijcusparse
194 : requires: cuda
195 : test:
196 : suffix: 2
197 : nsize: 2
198 : args: -bv_type {{vecs contiguous svec mat}shared output}
199 : test:
200 : suffix: 2_svec_vecs
201 : nsize: 2
202 : args: -bv_type svec -bv_matmult vecs
203 :
204 : TEST*/
|
