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 BVGetSplitRows().\n\n";
12 :
13 : #include <slepcbv.h>
14 :
15 12 : int main(int argc,char **argv)
16 : {
17 12 : Mat M,T,D,block[4];
18 12 : Vec t,v,v1,v2,w,*C;
19 12 : BV X,U,L;
20 12 : IS is[2];
21 12 : PetscInt i,j,n=10,k=5,l=3,nc=0,rstart,rend;
22 12 : PetscViewer view;
23 12 : PetscBool verbose;
24 :
25 12 : PetscFunctionBeginUser;
26 12 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
27 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
28 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
29 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL));
30 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL));
31 12 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
32 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BVGetSplitRows (length %" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT ", nc=%" PetscInt_FMT ").\n",n,l,k,nc));
33 :
34 : /* Create Nest matrix */
35 12 : PetscCall(MatCreate(PETSC_COMM_WORLD,&T));
36 12 : PetscCall(MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n));
37 12 : PetscCall(MatSetFromOptions(T));
38 12 : PetscCall(MatGetOwnershipRange(T,&rstart,&rend));
39 92 : for (i=rstart;i<rend;i++) {
40 80 : if (i>0) PetscCall(MatSetValue(T,i,i-1,-1.0,INSERT_VALUES));
41 80 : PetscCall(MatSetValue(T,i,i,2.0,INSERT_VALUES));
42 80 : if (i<n-1) PetscCall(MatSetValue(T,i,i+1,-1.0,INSERT_VALUES));
43 : }
44 12 : PetscCall(MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY));
45 12 : PetscCall(MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY));
46 12 : PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,&D));
47 :
48 12 : block[0] = T;
49 12 : block[1] = block[2] = NULL;
50 12 : block[3] = D;
51 12 : PetscCall(MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,block,&M));
52 12 : PetscCall(MatDestroy(&T));
53 12 : PetscCall(MatDestroy(&D));
54 12 : PetscCall(MatNestGetISs(M,is,NULL));
55 :
56 12 : PetscCall(MatView(M,NULL));
57 :
58 : /* Create template vector */
59 12 : PetscCall(MatCreateVecs(M,&t,NULL));
60 :
61 : /* Create BV object X */
62 12 : PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
63 12 : PetscCall(PetscObjectSetName((PetscObject)X,"X"));
64 12 : PetscCall(BVSetSizesFromVec(X,t,k));
65 12 : PetscCall(BVSetFromOptions(X));
66 :
67 : /* Generate constraints and attach them to X */
68 12 : if (nc>0) {
69 6 : PetscCall(VecDuplicateVecs(t,nc,&C));
70 18 : for (j=0;j<nc;j++) {
71 30 : for (i=0;i<=j;i++) PetscCall(VecSetValue(C[j],nc-i+1,1.0,INSERT_VALUES));
72 12 : PetscCall(VecAssemblyBegin(C[j]));
73 12 : PetscCall(VecAssemblyEnd(C[j]));
74 : }
75 6 : PetscCall(BVInsertConstraints(X,&nc,C));
76 6 : PetscCall(VecDestroyVecs(nc,&C));
77 : }
78 :
79 : /* Set up viewer */
80 12 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
81 12 : if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
82 :
83 : /* Fill X entries */
84 72 : for (j=0;j<k;j++) {
85 60 : PetscCall(BVGetColumn(X,j,&v));
86 60 : PetscCall(VecSet(v,0.0));
87 60 : PetscCall(VecGetSubVector(v,is[0],&v1));
88 60 : PetscCall(VecGetSubVector(v,is[1],&v2));
89 300 : for (i=0;i<4;i++) {
90 240 : if (i+j>=rstart && i+j<rend) {
91 160 : PetscCall(VecSetValue(v1,i+j,(PetscScalar)(3*i+j-2),INSERT_VALUES));
92 240 : PetscCall(VecSetValue(v2,i+j,(PetscScalar)(-i+2*(j+1)),INSERT_VALUES));
93 : }
94 : }
95 60 : PetscCall(VecAssemblyBegin(v1));
96 60 : PetscCall(VecAssemblyBegin(v2));
97 60 : PetscCall(VecAssemblyEnd(v1));
98 60 : PetscCall(VecAssemblyEnd(v2));
99 60 : PetscCall(VecRestoreSubVector(v,is[0],&v1));
100 60 : PetscCall(VecRestoreSubVector(v,is[1],&v2));
101 60 : PetscCall(BVRestoreColumn(X,j,&v));
102 : }
103 12 : if (verbose) PetscCall(BVView(X,view));
104 :
105 : /* Get split BVs */
106 12 : PetscCall(BVGetSplitRows(X,is[0],is[1],&U,&L));
107 12 : PetscCall(PetscObjectSetName((PetscObject)U,"U"));
108 12 : PetscCall(PetscObjectSetName((PetscObject)L,"L"));
109 :
110 12 : if (verbose) {
111 0 : PetscCall(BVView(U,view));
112 0 : PetscCall(BVView(L,view));
113 : }
114 :
115 : /* Copy l-th column of U to first column of L */
116 12 : PetscCall(BVGetColumn(U,l,&v));
117 12 : PetscCall(BVGetColumn(L,0,&w));
118 12 : PetscCall(VecCopy(v,w));
119 12 : PetscCall(BVRestoreColumn(U,l,&v));
120 12 : PetscCall(BVRestoreColumn(L,0,&w));
121 :
122 : /* Finished using the split BVs */
123 12 : PetscCall(BVRestoreSplitRows(X,is[0],is[1],&U,&L));
124 12 : if (verbose) PetscCall(BVView(X,view));
125 :
126 : /* Check: print bottom part of first column */
127 12 : PetscCall(BVGetColumn(X,0,&v));
128 12 : PetscCall(VecGetSubVector(v,is[1],&v2));
129 12 : PetscCall(VecView(v2,NULL));
130 12 : PetscCall(VecRestoreSubVector(v,is[1],&v2));
131 12 : PetscCall(BVRestoreColumn(X,0,&v));
132 :
133 12 : PetscCall(BVDestroy(&X));
134 12 : PetscCall(VecDestroy(&t));
135 12 : PetscCall(MatDestroy(&M));
136 12 : PetscCall(SlepcFinalize());
137 : return 0;
138 : }
139 :
140 : /*TEST
141 :
142 : testset:
143 : nsize: {{1 2}}
144 : output_file: output/test19_1.out
145 : filter: grep -v Process | grep -v Object | sed -e 's/mpi/seq/' | sed -e 's/seqcuda/seq/' | sed -e 's/seqaijcusparse/seqaij/' | sed -e 's/nc=2/nc=0/'
146 : test:
147 : suffix: 1
148 : args: -nc {{0 2}} -bv_type {{svec mat}}
149 : test:
150 : suffix: 1_cuda
151 : args: -nc {{0 2}} -bv_type {{svec mat}} -mat_type aijcusparse
152 : requires: cuda
153 :
154 : TEST*/
|
