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