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 BV orthogonalization functions with constraints.\n\n";
      12             : 
      13             : #include <slepcbv.h>
      14             : 
      15          16 : int main(int argc,char **argv)
      16             : {
      17          16 :   BV             X;
      18          16 :   Mat            M;
      19          16 :   Vec            v,t,*C;
      20          16 :   PetscInt       i,j,n=20,k=8,nc=2,nloc;
      21          16 :   PetscViewer    view;
      22          16 :   PetscBool      verbose;
      23          16 :   PetscReal      norm;
      24          16 :   PetscScalar    alpha;
      25             : 
      26          16 :   PetscFunctionBeginUser;
      27          16 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      28          16 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      29          16 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
      30          16 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL));
      31          16 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
      32          16 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV orthogonalization with %" PetscInt_FMT " columns + %" PetscInt_FMT " constraints, of length %" PetscInt_FMT ".\n",k,nc,n));
      33             : 
      34             :   /* Create template vector */
      35          16 :   PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
      36          16 :   PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
      37          16 :   PetscCall(VecSetFromOptions(t));
      38          16 :   PetscCall(VecGetLocalSize(t,&nloc));
      39             : 
      40             :   /* Create BV object X */
      41          16 :   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
      42          16 :   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
      43          16 :   PetscCall(BVSetSizes(X,nloc,n,k));
      44          16 :   PetscCall(BVSetFromOptions(X));
      45             : 
      46             :   /* Generate constraints and attach them to X */
      47          16 :   if (nc>0) {
      48          16 :     PetscCall(VecDuplicateVecs(t,nc,&C));
      49          48 :     for (j=0;j<nc;j++) {
      50          80 :       for (i=0;i<=j;i++) PetscCall(VecSetValue(C[j],i,1.0,INSERT_VALUES));
      51          32 :       PetscCall(VecAssemblyBegin(C[j]));
      52          32 :       PetscCall(VecAssemblyEnd(C[j]));
      53             :     }
      54          16 :     PetscCall(BVInsertConstraints(X,&nc,C));
      55          16 :     PetscCall(VecDestroyVecs(nc,&C));
      56             :   }
      57             : 
      58             :   /* Set up viewer */
      59          16 :   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
      60          16 :   if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
      61             : 
      62             :   /* Fill X entries */
      63         144 :   for (j=0;j<k;j++) {
      64         128 :     PetscCall(BVGetColumn(X,j,&v));
      65         128 :     PetscCall(VecSet(v,0.0));
      66        1536 :     for (i=0;i<=n/2;i++) {
      67        1408 :       if (i+j<n) {
      68        1408 :         alpha = (3.0*i+j-2)/(2*(i+j+1));
      69        1408 :         PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
      70             :       }
      71             :     }
      72         128 :     PetscCall(VecAssemblyBegin(v));
      73         128 :     PetscCall(VecAssemblyEnd(v));
      74         128 :     PetscCall(BVRestoreColumn(X,j,&v));
      75             :   }
      76          16 :   if (verbose) PetscCall(BVView(X,view));
      77             : 
      78             :   /* Test BVOrthogonalizeColumn */
      79         144 :   for (j=0;j<k;j++) {
      80         128 :     PetscCall(BVOrthogonalizeColumn(X,j,NULL,&norm,NULL));
      81         128 :     alpha = 1.0/norm;
      82         128 :     PetscCall(BVScaleColumn(X,j,alpha));
      83             :   }
      84          16 :   if (verbose) PetscCall(BVView(X,view));
      85             : 
      86             :   /* Check orthogonality */
      87          16 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M));
      88          16 :   PetscCall(BVDot(X,X,M));
      89          16 :   PetscCall(MatShift(M,-1.0));
      90          16 :   PetscCall(MatNorm(M,NORM_1,&norm));
      91          16 :   if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n"));
      92           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm));
      93             : 
      94          16 :   PetscCall(MatDestroy(&M));
      95          16 :   PetscCall(BVDestroy(&X));
      96          16 :   PetscCall(VecDestroy(&t));
      97          16 :   PetscCall(SlepcFinalize());
      98             :   return 0;
      99             : }
     100             : 
     101             : /*TEST
     102             : 
     103             :    testset:
     104             :       output_file: output/test6_1.out
     105             :       test:
     106             :          suffix: 1
     107             :          args: -bv_type {{vecs contiguous svec mat}shared output}
     108             :       test:
     109             :          suffix: 1_cuda
     110             :          args: -bv_type {{svec mat}} -vec_type cuda
     111             :          requires: cuda
     112             :       test:
     113             :          suffix: 2
     114             :          nsize: 2
     115             :          args: -bv_type {{vecs contiguous svec mat}shared output}
     116             :       test:
     117             :          suffix: 3
     118             :          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_type mgs
     119             :       test:
     120             :          suffix: 3_cuda
     121             :          args: -bv_type {{svec mat}} -vec_type cuda -bv_orthog_type mgs
     122             :          requires: cuda
     123             : 
     124             : TEST*/