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