Actual source code: test2.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2016, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Test BV orthogonalization functions.\n\n";
24: #include <slepcbv.h>
28: int main(int argc,char **argv)
29: {
31: BV X,Y,Z;
32: Mat M,R;
33: Vec v,t,e;
34: PetscInt i,j,n=20,k=8;
35: PetscViewer view;
36: PetscBool verbose;
37: PetscReal norm;
38: PetscScalar alpha;
40: SlepcInitialize(&argc,&argv,(char*)0,help);
41: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
42: PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
43: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
44: PetscPrintf(PETSC_COMM_WORLD,"Test BV orthogonalization with %D columns of length %D.\n",k,n);
46: /* Create template vector */
47: VecCreate(PETSC_COMM_WORLD,&t);
48: VecSetSizes(t,PETSC_DECIDE,n);
49: VecSetFromOptions(t);
51: /* Create BV object X */
52: BVCreate(PETSC_COMM_WORLD,&X);
53: PetscObjectSetName((PetscObject)X,"X");
54: BVSetSizesFromVec(X,t,k);
55: BVSetFromOptions(X);
57: /* Set up viewer */
58: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
59: if (verbose) {
60: PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
61: }
63: /* Fill X entries */
64: for (j=0;j<k;j++) {
65: BVGetColumn(X,j,&v);
66: VecSet(v,0.0);
67: for (i=0;i<=n/2;i++) {
68: if (i+j<n) {
69: alpha = (3.0*i+j-2)/(2*(i+j+1));
70: VecSetValue(v,i+j,alpha,INSERT_VALUES);
71: }
72: }
73: VecAssemblyBegin(v);
74: VecAssemblyEnd(v);
75: BVRestoreColumn(X,j,&v);
76: }
77: if (verbose) {
78: BVView(X,view);
79: }
81: /* Create copies on Y and Z */
82: BVDuplicate(X,&Y);
83: PetscObjectSetName((PetscObject)Y,"Y");
84: BVCopy(X,Y);
85: BVDuplicate(X,&Z);
86: PetscObjectSetName((PetscObject)Z,"Z");
87: BVCopy(X,Z);
89: /* Test BVOrthogonalizeColumn */
90: for (j=0;j<k;j++) {
91: BVOrthogonalizeColumn(X,j,NULL,&norm,NULL);
92: alpha = 1.0/norm;
93: BVScaleColumn(X,j,alpha);
94: }
95: if (verbose) {
96: BVView(X,view);
97: }
99: /* Check orthogonality */
100: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M);
101: BVDot(X,X,M);
102: MatShift(M,-1.0);
103: MatNorm(M,NORM_1,&norm);
104: if (norm<100*PETSC_MACHINE_EPSILON) {
105: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
106: } else {
107: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm);
108: }
110: /* Test BVOrthogonalize */
111: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&R);
112: PetscObjectSetName((PetscObject)R,"R");
113: BVOrthogonalize(Y,R);
114: if (verbose) {
115: BVView(Y,view);
116: MatView(R,view);
117: }
119: /* Check orthogonality */
120: BVDot(Y,Y,M);
121: MatShift(M,-1.0);
122: MatNorm(M,NORM_1,&norm);
123: if (norm<100*PETSC_MACHINE_EPSILON) {
124: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
125: } else {
126: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm);
127: }
129: /* Check residual */
130: BVMult(Z,-1.0,1.0,Y,R);
131: BVNorm(Z,NORM_FROBENIUS,&norm);
132: if (norm<100*PETSC_MACHINE_EPSILON) {
133: PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR|| < 100*eps\n");
134: } else {
135: PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR||: %g\n",(double)norm);
136: }
138: /* Test BVOrthogonalizeVec */
139: VecDuplicate(t,&e);
140: VecSet(e,1.0);
141: BVOrthogonalizeVec(X,e,NULL,&norm,NULL);
142: PetscPrintf(PETSC_COMM_WORLD,"Norm of ones(n,1) after orthogonalizing against X: %g\n",(double)norm);
144: MatDestroy(&M);
145: MatDestroy(&R);
146: BVDestroy(&X);
147: BVDestroy(&Y);
148: BVDestroy(&Z);
149: VecDestroy(&e);
150: VecDestroy(&t);
151: SlepcFinalize();
152: return ierr;
153: }