Actual source code: test15.c
slepc-3.22.1 2024-10-28
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 BVGetSplit().\n\n";
13: #include <slepcbv.h>
15: /*
16: Print the first row of a BV
17: */
18: PetscErrorCode PrintFirstRow(BV X)
19: {
20: PetscMPIInt rank;
21: PetscInt i,ld,k,nc;
22: const PetscScalar *pX;
23: const char *name;
25: PetscFunctionBeginUser;
26: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)X),&rank));
27: if (!rank) {
28: PetscCall(BVGetActiveColumns(X,NULL,&k));
29: PetscCall(BVGetLeadingDimension(X,&ld));
30: PetscCall(BVGetNumConstraints(X,&nc));
31: PetscCall(PetscObjectGetName((PetscObject)X,&name));
32: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"First row of %s =\n",name));
33: PetscCall(BVGetArrayRead(X,&pX));
34: for (i=0;i<nc+k;i++) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"%g ",(double)PetscRealPart(pX[i*ld])));
35: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)X),"\n"));
36: PetscCall(BVRestoreArrayRead(X,&pX));
37: }
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: int main(int argc,char **argv)
42: {
43: Vec t,v,*C;
44: BV X,L,R;
45: PetscInt i,j,n=10,k=5,l=3,nc=0;
46: PetscReal norm;
47: PetscScalar alpha;
48: PetscViewer view;
49: PetscBool verbose;
51: PetscFunctionBeginUser;
52: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
53: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
54: PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
55: PetscCall(PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL));
56: PetscCall(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL));
57: PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
58: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BVGetSplit (length %" PetscInt_FMT ", l=%" PetscInt_FMT ", k=%" PetscInt_FMT ", nc=%" PetscInt_FMT ").\n",n,l,k,nc));
60: /* Create template vector */
61: PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
62: PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
63: PetscCall(VecSetFromOptions(t));
65: /* Create BV object X */
66: PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
67: PetscCall(PetscObjectSetName((PetscObject)X,"X"));
68: PetscCall(BVSetSizesFromVec(X,t,k));
69: PetscCall(BVSetFromOptions(X));
71: /* Generate constraints and attach them to X */
72: if (nc>0) {
73: PetscCall(VecDuplicateVecs(t,nc,&C));
74: for (j=0;j<nc;j++) {
75: for (i=0;i<=j;i++) PetscCall(VecSetValue(C[j],nc-i+1,1.0,INSERT_VALUES));
76: PetscCall(VecAssemblyBegin(C[j]));
77: PetscCall(VecAssemblyEnd(C[j]));
78: }
79: PetscCall(BVInsertConstraints(X,&nc,C));
80: PetscCall(VecDestroyVecs(nc,&C));
81: }
83: /* Set up viewer */
84: PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
85: if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
87: /* Fill X entries */
88: for (j=0;j<k;j++) {
89: PetscCall(BVGetColumn(X,j,&v));
90: PetscCall(VecSet(v,0.0));
91: for (i=0;i<4;i++) {
92: if (i+j<n) PetscCall(VecSetValue(v,i+j,(PetscScalar)(3*i+j-2),INSERT_VALUES));
93: }
94: PetscCall(VecAssemblyBegin(v));
95: PetscCall(VecAssemblyEnd(v));
96: PetscCall(BVRestoreColumn(X,j,&v));
97: }
98: if (verbose) PetscCall(BVView(X,view));
100: /* Get split BVs */
101: PetscCall(BVSetActiveColumns(X,l,k));
102: PetscCall(BVGetSplit(X,&L,&R));
103: PetscCall(PetscObjectSetName((PetscObject)L,"L"));
104: PetscCall(PetscObjectSetName((PetscObject)R,"R"));
106: if (verbose) {
107: PetscCall(BVView(L,view));
108: PetscCall(BVView(R,view));
109: }
111: /* Modify first column of R */
112: PetscCall(BVGetColumn(R,0,&v));
113: PetscCall(VecSet(v,-1.0));
114: PetscCall(BVRestoreColumn(R,0,&v));
116: /* Finished using the split BVs */
117: PetscCall(BVRestoreSplit(X,&L,&R));
118: PetscCall(PrintFirstRow(X));
119: if (verbose) PetscCall(BVView(X,view));
121: /* Get the left split BV only */
122: PetscCall(BVGetSplit(X,&L,NULL));
123: for (j=0;j<l;j++) {
124: PetscCall(BVOrthogonalizeColumn(L,j,NULL,&norm,NULL));
125: alpha = 1.0/norm;
126: PetscCall(BVScaleColumn(L,j,alpha));
127: }
128: PetscCall(BVRestoreSplit(X,&L,NULL));
129: PetscCall(PrintFirstRow(X));
130: if (verbose) PetscCall(BVView(X,view));
132: /* Now get the right split BV after changing the number of leading columns */
133: PetscCall(BVSetActiveColumns(X,l-1,k));
134: PetscCall(BVGetSplit(X,NULL,&R));
135: PetscCall(BVGetColumn(R,0,&v));
136: PetscCall(BVInsertVec(X,0,v));
137: PetscCall(BVRestoreColumn(R,0,&v));
138: PetscCall(BVRestoreSplit(X,NULL,&R));
139: PetscCall(PrintFirstRow(X));
140: if (verbose) PetscCall(BVView(X,view));
142: PetscCall(BVDestroy(&X));
143: PetscCall(VecDestroy(&t));
144: PetscCall(SlepcFinalize());
145: return 0;
146: }
148: /*TEST
150: testset:
151: nsize: 2
152: output_file: output/test15_1.out
153: test:
154: suffix: 1
155: args: -bv_type {{vecs contiguous svec mat}shared output}
156: test:
157: suffix: 1_cuda
158: args: -bv_type {{svec mat}} -vec_type cuda
159: requires: cuda
160: test:
161: suffix: 1_hip
162: args: -bv_type {{svec mat}} -vec_type hip
163: requires: hip
165: testset:
166: nsize: 2
167: output_file: output/test15_2.out
168: test:
169: suffix: 2
170: args: -nc 2 -bv_type {{vecs contiguous svec mat}shared output}
171: test:
172: suffix: 2_cuda
173: args: -nc 2 -bv_type {{svec mat}} -vec_type cuda
174: requires: cuda
175: test:
176: suffix: 2_hip
177: args: -nc 2 -bv_type {{svec mat}} -vec_type hip
178: requires: hip
180: TEST*/