Actual source code: test11.c

slepc-3.8.3 2018-04-03
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2017, 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 block orthogonalization.\n\n";

 13: #include <slepcbv.h>

 15: int main(int argc,char **argv)
 16: {
 17:   PetscErrorCode    ierr;
 18:   BV                X,Y,Z,cached;
 19:   Mat               B,M;
 20:   Vec               v,t;
 21:   PetscInt          i,j,n=20,l=2,k=8,Istart,Iend;
 22:   PetscViewer       view;
 23:   PetscBool         verbose;
 24:   PetscReal         norm;
 25:   PetscScalar       alpha;
 26:   BVOrthogBlockType btype;

 28:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 29:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 30:   PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL);
 31:   PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
 32:   PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
 33:   PetscPrintf(PETSC_COMM_WORLD,"Test BV block orthogonalization (length %D, l=%D, k=%D).\n",n,l,k);

 35:   /* Create template vector */
 36:   VecCreate(PETSC_COMM_WORLD,&t);
 37:   VecSetSizes(t,PETSC_DECIDE,n);
 38:   VecSetFromOptions(t);

 40:   /* Create BV object X */
 41:   BVCreate(PETSC_COMM_WORLD,&X);
 42:   PetscObjectSetName((PetscObject)X,"X");
 43:   BVSetSizesFromVec(X,t,k);
 44:   BVSetFromOptions(X);
 45:   BVGetOrthogonalization(X,NULL,NULL,NULL,&btype);

 47:   /* Set up viewer */
 48:   PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
 49:   if (verbose) {
 50:     PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
 51:   }

 53:   /* Fill X entries */
 54:   for (j=0;j<k;j++) {
 55:     BVGetColumn(X,j,&v);
 56:     VecSet(v,0.0);
 57:     for (i=0;i<=n/2;i++) {
 58:       if (i+j<n) {
 59:         alpha = (3.0*i+j-2)/(2*(i+j+1));
 60:         VecSetValue(v,i+j,alpha,INSERT_VALUES);
 61:       }
 62:     }
 63:     VecAssemblyBegin(v);
 64:     VecAssemblyEnd(v);
 65:     BVRestoreColumn(X,j,&v);
 66:   }
 67:   if (btype==BV_ORTHOG_BLOCK_GS) {  /* GS requires the leading columns to be orthogonal */
 68:     for (j=0;j<l;j++) {
 69:       BVOrthonormalizeColumn(X,j,PETSC_FALSE,NULL,NULL);
 70:     }
 71:   }
 72:   if (verbose) {
 73:     BVView(X,view);
 74:   }

 76:   /* Create copy on Y */
 77:   BVDuplicate(X,&Y);
 78:   PetscObjectSetName((PetscObject)Y,"Y");
 79:   BVCopy(X,Y);
 80:   BVSetActiveColumns(Y,l,k);
 81:   BVSetActiveColumns(X,l,k);

 83:   /* Test BVOrthogonalize */
 84:   BVOrthogonalize(Y,NULL);
 85:   if (verbose) {
 86:     BVView(Y,view);
 87:   }

 89:   /* Check orthogonality */
 90:   MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M);
 91:   MatShift(M,1.0);   /* set leading part to identity */
 92:   BVDot(Y,Y,M);
 93:   MatShift(M,-1.0);
 94:   MatNorm(M,NORM_1,&norm);
 95:   if (norm<100*PETSC_MACHINE_EPSILON) {
 96:     PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
 97:   } else {
 98:     PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm);
 99:   }

101:   /* Create inner product matrix */
102:   MatCreate(PETSC_COMM_WORLD,&B);
103:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
104:   MatSetFromOptions(B);
105:   MatSetUp(B);
106:   PetscObjectSetName((PetscObject)B,"B");

108:   MatGetOwnershipRange(B,&Istart,&Iend);
109:   for (i=Istart;i<Iend;i++) {
110:     if (i>0) { MatSetValue(B,i,i-1,-1.0,INSERT_VALUES); }
111:     if (i<n-1) { MatSetValue(B,i,i+1,-1.0,INSERT_VALUES); }
112:     MatSetValue(B,i,i,2.0,INSERT_VALUES);
113:   }
114:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
115:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

117:   /* Prepare to repeat test, now with a non-standard inner product */
118:   if (btype!=BV_ORTHOG_BLOCK_TSQR) {  /* TSQR does not work with B-matrix */

120:     BVCopy(X,Y);
121:     BVDuplicate(X,&Z);
122:     PetscObjectSetName((PetscObject)Z,"Z");
123:     BVSetActiveColumns(Z,l,k);
124:     BVSetMatrix(X,B,PETSC_FALSE);
125:     BVSetMatrix(Y,B,PETSC_FALSE);
126:     if (btype==BV_ORTHOG_BLOCK_GS) {  /* GS requires the leading columns to be orthogonal */
127:       for (j=0;j<l;j++) {
128:         BVOrthonormalizeColumn(Y,j,PETSC_FALSE,NULL,NULL);
129:       }
130:     }
131:     if (verbose) {
132:       BVView(X,view);
133:     }

135:     /* Test BVOrthogonalize */
136:     BVOrthogonalize(Y,NULL);
137:     if (verbose) {
138:       BVView(Y,view);
139:     }

141:     /* Extract cached BV and check it is equal to B*X */
142:     BVGetCachedBV(Y,&cached);
143:     BVMatMult(X,B,Z);
144:     BVMult(Z,-1.0,1.0,cached,NULL);
145:     BVNorm(Z,NORM_FROBENIUS,&norm);
146:     if (norm<100*PETSC_MACHINE_EPSILON) {
147:       PetscPrintf(PETSC_COMM_WORLD,"Residual ||cached-BX|| < 100*eps\n");
148:     } else {
149:       PetscPrintf(PETSC_COMM_WORLD,"Residual ||cached-BX||: %g\n",(double)norm);
150:     }

152:     /* Check orthogonality */
153:     MatZeroEntries(M);
154:     MatShift(M,1.0);   /* set leading part to identity */
155:     BVDot(Y,Y,M);
156:     MatShift(M,-1.0);
157:     MatNorm(M,NORM_1,&norm);
158:     if (norm<100*PETSC_MACHINE_EPSILON) {
159:       PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
160:     } else {
161:       PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm);
162:     }
163:     BVDestroy(&Z);
164:   }

166:   MatDestroy(&M);
167:   MatDestroy(&B);
168:   BVDestroy(&X);
169:   BVDestroy(&Y);
170:   VecDestroy(&t);
171:   SlepcFinalize();
172:   return ierr;
173: }