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 bi-orthogonalization functions.\n\n";
      12             : 
      13             : #include <slepcbv.h>
      14             : 
      15           8 : int main(int argc,char **argv)
      16             : {
      17           8 :   BV             X,Y;
      18           8 :   Mat            M;
      19           8 :   Vec            v,t;
      20           8 :   PetscInt       i,j,n=20,k=7;
      21           8 :   PetscViewer    view;
      22           8 :   PetscBool      verbose;
      23           8 :   PetscReal      norm,condn=1.0;
      24           8 :   PetscScalar    alpha;
      25             : 
      26           8 :   PetscFunctionBeginUser;
      27           8 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      28           8 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      29           8 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
      30           8 :   PetscCall(PetscOptionsGetReal(NULL,NULL,"-condn",&condn,NULL));
      31           8 :   PetscCheck(condn>=1.0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"The condition number must be > 1");
      32           8 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
      33           8 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV bi-orthogonalization with %" PetscInt_FMT " columns of length %" PetscInt_FMT ".\n",k,n));
      34           8 :   if (condn>1.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," - Using random BVs with condition number = %g\n",(double)condn));
      35             : 
      36             :   /* Create template vector */
      37           8 :   PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
      38           8 :   PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
      39           8 :   PetscCall(VecSetFromOptions(t));
      40             : 
      41             :   /* Create BV object X */
      42           8 :   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
      43           8 :   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
      44           8 :   PetscCall(BVSetSizesFromVec(X,t,k));
      45           8 :   PetscCall(BVSetFromOptions(X));
      46             : 
      47             :   /* Set up viewer */
      48           8 :   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
      49           8 :   if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
      50             : 
      51             :   /* Fill X entries */
      52           8 :   if (condn==1.0) {
      53          64 :     for (j=0;j<k;j++) {
      54          56 :       PetscCall(BVGetColumn(X,j,&v));
      55          56 :       PetscCall(VecSet(v,0.0));
      56         672 :       for (i=0;i<=n/2;i++) {
      57         616 :         if (i+j<n) {
      58         616 :           alpha = (3.0*i+j-2)/(2*(i+j+1));
      59             : #if defined(PETSC_USE_COMPLEX)
      60             :           alpha += (1.2*i+j-2)/(0.1*(i+j+1))*PETSC_i;
      61             : #endif
      62         616 :           PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
      63             :         }
      64             :       }
      65          56 :       PetscCall(VecAssemblyBegin(v));
      66          56 :       PetscCall(VecAssemblyEnd(v));
      67          56 :       PetscCall(BVRestoreColumn(X,j,&v));
      68             :     }
      69           0 :   } else PetscCall(BVSetRandomCond(X,condn));
      70           8 :   if (verbose) PetscCall(BVView(X,view));
      71             : 
      72             :   /* Create Y and fill its entries */
      73           8 :   PetscCall(BVDuplicate(X,&Y));
      74           8 :   PetscCall(PetscObjectSetName((PetscObject)Y,"Y"));
      75           8 :   if (condn==1.0) {
      76          64 :     for (j=0;j<k;j++) {
      77          56 :       PetscCall(BVGetColumn(Y,j,&v));
      78          56 :       PetscCall(VecSet(v,0.0));
      79         384 :       for (i=PetscMin(n,2+(2*j)%6);i<PetscMin(n,6+(3*j)%9);i++) {
      80         272 :         if (i%5 && i!=j) {
      81         200 :           alpha = (1.5*i+j)/(2.2*(i-j));
      82         272 :           PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
      83             :         }
      84             :       }
      85          56 :       PetscCall(VecAssemblyBegin(v));
      86          56 :       PetscCall(VecAssemblyEnd(v));
      87          56 :       PetscCall(BVRestoreColumn(Y,j,&v));
      88             :     }
      89           0 :   } else PetscCall(BVSetRandomCond(Y,condn));
      90           8 :   if (verbose) PetscCall(BVView(Y,view));
      91             : 
      92             :   /* Test BVBiorthonormalizeColumn */
      93          64 :   for (j=0;j<k;j++) PetscCall(BVBiorthonormalizeColumn(X,Y,j,NULL));
      94           8 :   if (verbose) {
      95           0 :     PetscCall(BVView(X,view));
      96           0 :     PetscCall(BVView(Y,view));
      97             :   }
      98             : 
      99             :   /* Check orthogonality */
     100           8 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M));
     101           8 :   PetscCall(PetscObjectSetName((PetscObject)M,"M"));
     102           8 :   PetscCall(BVDot(X,Y,M));
     103           8 :   if (verbose) PetscCall(MatView(M,view));
     104           8 :   PetscCall(MatShift(M,-1.0));
     105           8 :   PetscCall(MatNorm(M,NORM_1,&norm));
     106           8 :   if (norm<200*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of bi-orthogonality < 200*eps\n"));
     107           0 :   else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of bi-orthogonality: %g\n",(double)norm));
     108             : 
     109           8 :   PetscCall(MatDestroy(&M));
     110           8 :   PetscCall(BVDestroy(&X));
     111           8 :   PetscCall(BVDestroy(&Y));
     112           8 :   PetscCall(VecDestroy(&t));
     113           8 :   PetscCall(SlepcFinalize());
     114             :   return 0;
     115             : }
     116             : 
     117             : /*TEST
     118             : 
     119             :    testset:
     120             :       output_file: output/test17_1.out
     121             :       test:
     122             :          suffix: 1
     123             :          args: -bv_type {{vecs contiguous svec mat}} -bv_orthog_type cgs
     124             :          requires: double
     125             :       test:
     126             :          suffix: 1_cuda
     127             :          args: -bv_type {{svec mat}} -vec_type cuda -bv_orthog_type cgs
     128             :          requires: cuda
     129             :       test:
     130             :          suffix: 2
     131             :          args: -bv_type {{vecs contiguous svec mat}} -bv_orthog_type mgs
     132             :       test:
     133             :          suffix: 2_cuda
     134             :          args: -bv_type {{svec mat}} -vec_type cuda -bv_orthog_type mgs
     135             :          requires: cuda
     136             : 
     137             : TEST*/