Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : static char help[] = "Test block orthogonalization on a rank-deficient BV.\n\n";
      12             : 
      13             : #include <slepcbv.h>
      14             : 
      15           4 : int main(int argc,char **argv)
      16             : {
      17           4 :   BV             X,Z;
      18           4 :   Mat            M,R;
      19           4 :   Vec            v,w,t;
      20           4 :   PetscInt       i,j,n=20,k=8;
      21           4 :   PetscViewer    view;
      22           4 :   PetscBool      verbose;
      23           4 :   PetscReal      norm;
      24           4 :   PetscScalar    alpha;
      25             : 
      26           4 :   PetscFunctionBeginUser;
      27           4 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      28           4 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      29           4 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
      30           4 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
      31           4 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV block orthogonalization (length %" PetscInt_FMT ", k=%" PetscInt_FMT ").\n",n,k));
      32           4 :   PetscCheck(k>5,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"k must be at least 6");
      33             : 
      34             :   /* Create template vector */
      35           4 :   PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
      36           4 :   PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
      37           4 :   PetscCall(VecSetFromOptions(t));
      38             : 
      39             :   /* Create BV object X */
      40           4 :   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
      41           4 :   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
      42           4 :   PetscCall(BVSetSizesFromVec(X,t,k));
      43           4 :   PetscCall(BVSetFromOptions(X));
      44             : 
      45             :   /* Set up viewer */
      46           4 :   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
      47           4 :   if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
      48             : 
      49             :   /* Fill X entries (first half) */
      50          20 :   for (j=0;j<k/2;j++) {
      51          16 :     PetscCall(BVGetColumn(X,j,&v));
      52          16 :     PetscCall(VecSet(v,0.0));
      53         192 :     for (i=0;i<=n/2;i++) {
      54         176 :       if (i+j<n) {
      55         176 :         alpha = (3.0*i+j-2)/(2*(i+j+1));
      56         176 :         PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
      57             :       }
      58             :     }
      59          16 :     PetscCall(VecAssemblyBegin(v));
      60          16 :     PetscCall(VecAssemblyEnd(v));
      61          16 :     PetscCall(BVRestoreColumn(X,j,&v));
      62             :   }
      63             : 
      64             :   /* make middle column linearly dependent wrt columns 0 and 1 */
      65           4 :   PetscCall(BVCopyColumn(X,0,j));
      66           4 :   PetscCall(BVGetColumn(X,j,&v));
      67           4 :   PetscCall(BVGetColumn(X,1,&w));
      68           4 :   PetscCall(VecAXPY(v,0.5,w));
      69           4 :   PetscCall(BVRestoreColumn(X,1,&w));
      70           4 :   PetscCall(BVRestoreColumn(X,j,&v));
      71           4 :   j++;
      72             : 
      73             :   /* Fill X entries (second half) */
      74          12 :   for (;j<k-1;j++) {
      75           8 :     PetscCall(BVGetColumn(X,j,&v));
      76           8 :     PetscCall(VecSet(v,0.0));
      77          96 :     for (i=0;i<=n/2;i++) {
      78          88 :       if (i+j<n) {
      79          88 :         alpha = (3.0*i+j-2)/(2*(i+j+1));
      80          88 :         PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
      81             :       }
      82             :     }
      83           8 :     PetscCall(VecAssemblyBegin(v));
      84           8 :     PetscCall(VecAssemblyEnd(v));
      85           8 :     PetscCall(BVRestoreColumn(X,j,&v));
      86             :   }
      87             : 
      88             :   /* make middle column linearly dependent wrt columns 1 and k/2+1 */
      89           4 :   PetscCall(BVCopyColumn(X,1,j));
      90           4 :   PetscCall(BVGetColumn(X,j,&v));
      91           4 :   PetscCall(BVGetColumn(X,k/2+1,&w));
      92           4 :   PetscCall(VecAXPY(v,-1.2,w));
      93           4 :   PetscCall(BVRestoreColumn(X,k/2+1,&w));
      94           4 :   PetscCall(BVRestoreColumn(X,j,&v));
      95             : 
      96           4 :   if (verbose) PetscCall(BVView(X,view));
      97             : 
      98             :   /* Create a copy on Z */
      99           4 :   PetscCall(BVDuplicate(X,&Z));
     100           4 :   PetscCall(PetscObjectSetName((PetscObject)Z,"Z"));
     101           4 :   PetscCall(BVCopy(X,Z));
     102             : 
     103             :   /* Test BVOrthogonalize */
     104           4 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&R));
     105           4 :   PetscCall(PetscObjectSetName((PetscObject)R,"R"));
     106           4 :   PetscCall(BVOrthogonalize(X,R));
     107           4 :   if (verbose) {
     108           0 :     PetscCall(BVView(X,view));
     109           0 :     PetscCall(MatView(R,view));
     110             :   }
     111             : 
     112             :   /* Check orthogonality */
     113           4 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M));
     114           4 :   PetscCall(MatShift(M,1.0));   /* set leading part to identity */
     115           4 :   PetscCall(BVDot(X,X,M));
     116           4 :   PetscCall(MatShift(M,-1.0));
     117           4 :   PetscCall(MatNorm(M,NORM_1,&norm));
     118           4 :   if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n"));
     119           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm));
     120             : 
     121             :   /* Check residual */
     122           4 :   PetscCall(BVMult(Z,-1.0,1.0,X,R));
     123           4 :   PetscCall(BVNorm(Z,NORM_FROBENIUS,&norm));
     124           4 :   if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR|| < 100*eps\n"));
     125           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR||: %g\n",(double)norm));
     126             : 
     127           4 :   PetscCall(MatDestroy(&R));
     128           4 :   PetscCall(MatDestroy(&M));
     129           4 :   PetscCall(BVDestroy(&X));
     130           4 :   PetscCall(BVDestroy(&Z));
     131           4 :   PetscCall(VecDestroy(&t));
     132           4 :   PetscCall(SlepcFinalize());
     133             :   return 0;
     134             : }
     135             : 
     136             : /*TEST
     137             : 
     138             :    test:
     139             :       suffix: 1
     140             :       nsize: 1
     141             :       args: -bv_orthog_block gs -bv_type {{vecs contiguous svec mat}shared output}
     142             :       output_file: output/test12_1.out
     143             : 
     144             : TEST*/