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 DSGNHEP with upper quasi-triangular matrix pair.\n\n";
      12             : 
      13             : #include <slepcds.h>
      14             : 
      15           1 : int main(int argc,char **argv)
      16             : {
      17           1 :   DS             ds;
      18           1 :   PetscScalar    *A,*B,*X;
      19           1 :   PetscReal      rnorm,aux;
      20           1 :   PetscInt       i,j,n=10,ld;
      21           1 :   PetscViewer    viewer;
      22           1 :   PetscBool      verbose;
      23             : 
      24           1 :   PetscFunctionBeginUser;
      25           1 :   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
      26           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      27           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type GNHEP - dimension %" PetscInt_FMT ".\n",n));
      28           1 :   PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
      29             : 
      30             :   /* Create DS object */
      31           1 :   PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
      32           1 :   PetscCall(DSSetType(ds,DSGNHEP));
      33           1 :   PetscCall(DSSetFromOptions(ds));
      34           1 :   ld = n+2;  /* test leading dimension larger than n */
      35           1 :   PetscCall(DSAllocate(ds,ld));
      36           1 :   PetscCall(DSSetDimensions(ds,n,0,0));
      37             : 
      38             :   /* Set up viewer */
      39           1 :   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
      40           1 :   PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
      41           1 :   PetscCall(DSView(ds,viewer));
      42           1 :   PetscCall(PetscViewerPopFormat(viewer));
      43           1 :   if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
      44             : 
      45             :   /* Fill A,B with upper quasi-triangular matrices */
      46           1 :   PetscCall(DSGetArray(ds,DS_MAT_A,&A));
      47           1 :   PetscCall(DSGetArray(ds,DS_MAT_B,&B));
      48           1 :   PetscCall(PetscArrayzero(A,ld*n));
      49          11 :   for (i=0;i<n;i++) A[i+i*ld]=2.0;
      50           3 :   for (j=1;j<3;j++) {
      51          19 :     for (i=0;i<n-j;i++) A[i+(i+j)*ld]=0.001;
      52             :   }
      53           1 :   PetscCall(PetscArrayzero(B,ld*n));
      54          11 :   for (i=0;i<n;i++) B[i+i*ld]=1.0;
      55           1 :   B[1+0*ld]=B[0+1*ld]=PETSC_MACHINE_EPSILON;
      56           4 :   for (i=1;i<n;i+=3) {
      57           3 :     A[i+(i-1)*ld]=-A[(i-1)+i*ld];
      58             :   }
      59           1 :   PetscCall(DSRestoreArray(ds,DS_MAT_A,&A));
      60           1 :   PetscCall(DSRestoreArray(ds,DS_MAT_B,&B));
      61           1 :   PetscCall(DSSetState(ds,DS_STATE_INTERMEDIATE));
      62             : 
      63           1 :   if (verbose) {
      64           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
      65           0 :     PetscCall(DSView(ds,viewer));
      66             :   }
      67             : 
      68             :   /* Eigenvectors */
      69           1 :   j = 0;
      70           1 :   PetscCall(DSVectors(ds,DS_MAT_X,&j,&rnorm));  /* first eigenvector */
      71           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Value of rnorm for 2nd vector = %.3f\n",(double)rnorm));
      72           1 :   PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));  /* all eigenvectors */
      73           1 :   j = 0;
      74           1 :   rnorm = 0.0;
      75           1 :   PetscCall(DSGetArray(ds,DS_MAT_X,&X));
      76          11 :   for (i=0;i<n;i++) {
      77             : #if defined(PETSC_USE_COMPLEX)
      78             :     aux = PetscAbsScalar(X[i+j*ld]);
      79             : #else
      80          10 :     aux = SlepcAbsEigenvalue(X[i+j*ld],X[i+(j+1)*ld]);
      81             : #endif
      82          10 :     rnorm += aux*aux;
      83             :   }
      84           1 :   PetscCall(DSRestoreArray(ds,DS_MAT_X,&X));
      85           1 :   rnorm = PetscSqrtReal(rnorm);
      86           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of 1st columns = %.3f\n",(double)rnorm));
      87           1 :   if (verbose) {
      88           0 :     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After vectors - - - - - - - - -\n"));
      89           0 :     PetscCall(DSView(ds,viewer));
      90             :   }
      91             : 
      92           1 :   PetscCall(DSDestroy(&ds));
      93           1 :   PetscCall(SlepcFinalize());
      94             :   return 0;
      95             : }
      96             : 
      97             : /*TEST
      98             : 
      99             :    test:
     100             :       suffix: 1
     101             :       filter: sed -e "s/-0\./0./"
     102             : 
     103             : TEST*/