Actual source code: test8.c
slepc-3.17.0 2022-03-31
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 BV orthogonalization with selected columns.\n\n";
13: #include <slepcbv.h>
15: int main(int argc,char **argv)
16: {
17: BV X;
18: Vec v,t,z;
19: PetscInt i,j,n=20,k=8;
20: PetscViewer view;
21: PetscBool verbose,*which;
22: PetscReal norm;
23: PetscScalar alpha,*pz;
25: SlepcInitialize(&argc,&argv,(char*)0,help);
26: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
27: PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
28: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
29: PetscPrintf(PETSC_COMM_WORLD,"Test BV orthogonalization with selected columns of length %" PetscInt_FMT ".\n",n);
31: /* Create template vector */
32: VecCreate(PETSC_COMM_WORLD,&t);
33: VecSetSizes(t,PETSC_DECIDE,n);
34: VecSetFromOptions(t);
36: /* Create BV object X */
37: BVCreate(PETSC_COMM_WORLD,&X);
38: PetscObjectSetName((PetscObject)X,"X");
39: BVSetSizesFromVec(X,t,k);
40: BVSetOrthogonalization(X,BV_ORTHOG_MGS,BV_ORTHOG_REFINE_IFNEEDED,PETSC_DEFAULT,BV_ORTHOG_BLOCK_GS);
41: BVSetFromOptions(X);
43: /* Set up viewer */
44: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
45: if (verbose) PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
47: /* Fill X entries */
48: for (j=0;j<k;j++) {
49: BVGetColumn(X,j,&v);
50: VecSet(v,0.0);
51: for (i=0;i<=n/2;i++) {
52: if (i+j<n) {
53: alpha = (3.0*i+j-2)/(2*(i+j+1));
54: VecSetValue(v,i+j,alpha,INSERT_VALUES);
55: }
56: }
57: VecAssemblyBegin(v);
58: VecAssemblyEnd(v);
59: BVRestoreColumn(X,j,&v);
60: }
61: if (verbose) BVView(X,view);
63: /* Orthonormalize first k-1 columns */
64: for (j=0;j<k-1;j++) {
65: BVOrthogonalizeColumn(X,j,NULL,&norm,NULL);
66: alpha = 1.0/norm;
67: BVScaleColumn(X,j,alpha);
68: }
69: if (verbose) BVView(X,view);
71: /* Select odd columns and orthogonalize last column against those only */
72: PetscMalloc1(k,&which);
73: for (i=0;i<k;i++) which[i] = (i%2)? PETSC_TRUE: PETSC_FALSE;
74: BVOrthogonalizeSomeColumn(X,k-1,which,NULL,NULL,NULL);
75: PetscFree(which);
76: if (verbose) BVView(X,view);
78: /* Check orthogonality */
79: PetscPrintf(PETSC_COMM_WORLD,"Orthogonalization coefficients:\n");
80: VecCreateSeq(PETSC_COMM_SELF,k-1,&z);
81: PetscObjectSetName((PetscObject)z,"z");
82: VecGetArray(z,&pz);
83: BVDotColumn(X,k-1,pz);
84: for (i=0;i<k-1;i++) {
85: if (PetscAbsScalar(pz[i])<5.0*PETSC_MACHINE_EPSILON) pz[i]=0.0;
86: }
87: VecRestoreArray(z,&pz);
88: VecView(z,view);
89: VecDestroy(&z);
91: BVDestroy(&X);
92: VecDestroy(&t);
93: SlepcFinalize();
94: return 0;
95: }
97: /*TEST
99: testset:
100: output_file: output/test8_1.out
101: test:
102: suffix: 1
103: args: -bv_type {{vecs contiguous svec mat}shared output}
104: test:
105: suffix: 1_cuda
106: args: -bv_type svec -vec_type cuda
107: requires: cuda
108: test:
109: suffix: 2
110: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_refine never
111: requires: !single
112: test:
113: suffix: 3
114: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_refine always
115: test:
116: suffix: 4
117: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_type mgs
118: test:
119: suffix: 4_cuda
120: args: -bv_type svec -vec_type cuda -bv_orthog_type mgs
121: requires: cuda
123: TEST*/