Actual source code: slepcimpl.h

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: #pragma once

     13: #include <slepcsys.h>
     14: #include <petsc/private/petscimpl.h>

     16: /* SUBMANSEC = Sys */

     18: SLEPC_INTERN PetscBool SlepcBeganPetsc;

     20: /* SlepcSwap - swap two variables a,b of the same type using a temporary variable t */
     21: #define SlepcSwap(a,b,t) do {t=a;a=b;b=t;} while (0)

     23: /*MC
     24:     SlepcHeaderCreate - Creates a SLEPc object

     26:     Input Parameters:
     27: +   classid    - the classid associated with this object
     28: .   class_name - string name of class; should be static
     29: .   descr      - string containing short description; should be static
     30: .   mansec     - string indicating section in manual pages; should be static
     31: .   comm       - the MPI Communicator
     32: .   destroy    - the destroy routine for this object
     33: -   view       - the view routine for this object

     35:     Output Parameter:
     36: .   h - the newly created object

     38:     Note:
     39:     This is equivalent to `PetscHeaderCreate()` but makes sure that `SlepcInitialize()`
     40:     has been called.

     42:     Level: developer
     43: M*/
     44: #define SlepcHeaderCreate(h,classid,class_name,descr,mansec,comm,destroy,view) \
     45:     ((PetscErrorCode)((!SlepcInitializeCalled && \
     46:                        PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,PETSC_ERR_ORDER,PETSC_ERROR_INITIAL, \
     47:                                   "Must call SlepcInitialize instead of PetscInitialize to use SLEPc classes")) || \
     48:                       PetscHeaderCreate(h,classid,class_name,descr,mansec,comm,destroy,view)))

     50: /* context for structured eigenproblem matrices created via MatCreateXXX */
     51: struct _n_SlepcMatStruct {
     52:   PetscInt    cookie;    /* identify which structured matrix */
     53:   PetscScalar s;         /* in BSE sign of the bottom part of the vector */
     54: };
     55: typedef struct _n_SlepcMatStruct* SlepcMatStruct;

     57: #define SLEPC_MAT_STRUCT_BSE     88101
     58: #define SLEPC_MAT_STRUCT_HAMILT  88102
     59: #define SLEPC_MAT_STRUCT_LREP    88103

     61: /*
     62:   SlepcCheckMatStruct - Check that a given Mat is a structured matrix of the wanted type.

     64:   Returns true/false in flg if it is given, otherwise yields an error if the check fails.
     65:   If cookie==0 it will check for any type.
     66: */
     67: static inline PetscErrorCode SlepcCheckMatStruct(Mat A,PetscInt cookie,PetscBool *flg)
     68: {
     69:   PetscContainer container;
     70:   SlepcMatStruct mctx;

     72:   PetscFunctionBegin;
     73:   if (flg) *flg = PETSC_FALSE;
     74:   PetscCall(PetscObjectQuery((PetscObject)A,"SlepcMatStruct",(PetscObject*)&container));
     75:   if (flg && !container) PetscFunctionReturn(PETSC_SUCCESS);
     76:   PetscCheck(container,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"The Mat is not a structured matrix");
     77:   if (cookie) {
     78:     PetscCall(PetscContainerGetPointer(container,(void**)&mctx));
     79:     if (flg && (!mctx || mctx->cookie!=cookie)) PetscFunctionReturn(PETSC_SUCCESS);
     80:     PetscCheck(mctx && mctx->cookie==cookie,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"The type of structured matrix is different from the expected one");
     81:   }
     82:   if (flg) *flg = PETSC_TRUE;
     83:   PetscFunctionReturn(PETSC_SUCCESS);
     84: }

     86: /*
     87:   SlepcCheckMatLREPReduced - Check whether a given LREP Mat is of reduced form or not.

     89:   True if the matrix has reduced form [0 K; M 0], otherwise it has the form [A B; -B -A].
     90: */
     91: static inline PetscErrorCode SlepcCheckMatLREPReduced(Mat A,PetscBool *flg)
     92: {
     93:   Mat A00;

     95:   PetscFunctionBegin;
     96:   PetscAssertPointer(flg,2);
     97:   *flg = PETSC_FALSE;
     98:   PetscCall(SlepcCheckMatStruct(A,SLEPC_MAT_STRUCT_LREP,NULL));  /* error out if not LREP */
     99:   PetscCall(MatNestGetSubMat(A,0,0,&A00));
    100:   if (A00==NULL) *flg = PETSC_TRUE;
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: /*
    105:   SlepcPrintEigenvalueASCII - Print an eigenvalue on an ASCII viewer.
    106: */
    107: static inline PetscErrorCode SlepcPrintEigenvalueASCII(PetscViewer viewer,PetscScalar eigr,PetscScalar eigi)
    108: {
    109:   PetscReal      re,im;

    111:   PetscFunctionBegin;
    112: #if defined(PETSC_USE_COMPLEX)
    113:   re = PetscRealPart(eigr);
    114:   im = PetscImaginaryPart(eigr);
    115:   (void)eigi;
    116: #else
    117:   re = eigr;
    118:   im = eigi;
    119: #endif
    120:   /* print zero instead of tiny value */
    121:   if (PetscAbs(im) && PetscAbs(re)/PetscAbs(im)<PETSC_SMALL) re = 0.0;
    122:   if (PetscAbs(re) && PetscAbs(im)/PetscAbs(re)<PETSC_SMALL) im = 0.0;
    123:   /* print as real if imaginary part is zero */
    124:   if (im!=(PetscReal)0.0) PetscCall(PetscViewerASCIIPrintf(viewer,"%.5f%+.5fi",(double)re,(double)im));
    125:   else PetscCall(PetscViewerASCIIPrintf(viewer,"%.5f",(double)re));
    126:   PetscFunctionReturn(PETSC_SUCCESS);
    127: }

    129: /*
    130:   SlepcViewEigenvector - Outputs an eigenvector xr,xi to a viewer.
    131:   In complex scalars only xr is written.
    132:   The name of xr,xi is set before writing, based on the label, the index, and the name of obj.
    133: */
    134: static inline PetscErrorCode SlepcViewEigenvector(PetscViewer viewer,Vec xr,Vec xi,const char *label,PetscInt index,PetscObject obj)
    135: {
    136:   size_t         count;
    137:   char           vname[30];
    138:   const char     *pname;

    140:   PetscFunctionBegin;
    141:   PetscCall(PetscObjectGetName(obj,&pname));
    142:   PetscCall(PetscSNPrintfCount(vname,sizeof(vname),"%s%s",&count,label,PetscDefined(USE_COMPLEX)?"":"r"));
    143:   count--;
    144:   PetscCall(PetscSNPrintf(vname+count,sizeof(vname)-count,"%" PetscInt_FMT "_%s",index,pname));
    145:   PetscCall(PetscObjectSetName((PetscObject)xr,vname));
    146:   PetscCall(VecView(xr,viewer));
    147:   if (!PetscDefined(USE_COMPLEX)) {
    148:     vname[count-1] = 'i';
    149:     PetscCall(PetscObjectSetName((PetscObject)xi,vname));
    150:     PetscCall(VecView(xi,viewer));
    151:   }
    152:   PetscFunctionReturn(PETSC_SUCCESS);
    153: }

    155: /* Macros for strings with different value in real and complex */
    156: #if defined(PETSC_USE_COMPLEX)
    157: #define SLEPC_STRING_HERMITIAN "hermitian"
    158: #else
    159: #define SLEPC_STRING_HERMITIAN "symmetric"
    160: #endif

    162: /* Private functions that are shared by several classes */
    163: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcBasisReference_Private(PetscInt,Vec*,PetscInt*,Vec**);
    164: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcBasisDestroy_Private(PetscInt*,Vec**);
    165: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode SlepcMonitorMakeKey_Internal(const char[],PetscViewerType,PetscViewerFormat,char[]);
    166: SLEPC_SINGLE_LIBRARY_INTERN PetscErrorCode PetscViewerAndFormatCreate_Internal(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**);

    168: SLEPC_INTERN PetscErrorCode SlepcCitationsInitialize(void);
    169: SLEPC_INTERN PetscErrorCode SlepcInitialize_DynamicLibraries(void);
    170: SLEPC_INTERN PetscErrorCode SlepcInitialize_Packages(void);

    172: /* Macro to check a sequential Mat (including GPU) */
    173: #if !defined(PETSC_USE_DEBUG)
    174: #define SlepcMatCheckSeq(h) do {(void)(h);} while (0)
    175: #else
    176: #if defined(PETSC_HAVE_CUDA)
    177: #define SlepcMatCheckSeq(h) do { PetscCheckTypeNames((h),MATSEQDENSE,MATSEQDENSECUDA); } while (0)
    178: #elif defined(PETSC_HAVE_HIP)
    179: #define SlepcMatCheckSeq(h) do { PetscCheckTypeNames((h),MATSEQDENSE,MATSEQDENSEHIP); } while (0)
    180: #else
    181: #define SlepcMatCheckSeq(h) do { PetscCheckTypeName((h),MATSEQDENSE); } while (0)
    182: #endif
    183: #endif

    185: /* Definitions needed to work with GPU kernels */
    186: #if defined(PETSC_HAVE_CUPM)
    187: #include <petscdevice_cupm.h>

    189: #define X_AXIS 0
    190: #define Y_AXIS 1

    192: #define SLEPC_TILE_SIZE_X  32
    193: #define SLEPC_BLOCK_SIZE_X 128
    194: #define SLEPC_TILE_SIZE_Y  32
    195: #define SLEPC_BLOCK_SIZE_Y 128

    197: static inline PetscErrorCode SlepcKernelSetGrid1D(PetscInt rows,dim3 *dimGrid,dim3 *dimBlock,PetscInt *dimGrid_xcount)
    198: {
    199:   int card;
    200: #if defined(PETSC_HAVE_CUDA)
    201:   struct cudaDeviceProp devprop;
    202: #elif defined(PETSC_HAVE_HIP)
    203:   hipDeviceProp_t devprop;
    204: #endif

    206:   PetscFunctionBegin;
    207: #if defined(PETSC_HAVE_CUDA)
    208:   PetscCallCUDA(cudaGetDevice(&card));
    209:   PetscCallCUDA(cudaGetDeviceProperties(&devprop,card));
    210: #elif defined(PETSC_HAVE_HIP)
    211:   PetscCallHIP(hipGetDevice(&card));
    212:   PetscCallHIP(hipGetDeviceProperties(&devprop,card));
    213: #endif
    214:   *dimGrid_xcount = 1;

    216:   /* X axis */
    217:   dimGrid->x  = 1;
    218:   dimBlock->x = SLEPC_BLOCK_SIZE_X;
    219:   if (rows>SLEPC_BLOCK_SIZE_X) dimGrid->x = (rows+SLEPC_BLOCK_SIZE_X-1)/SLEPC_BLOCK_SIZE_X;
    220:   else dimBlock->x = rows;
    221:   if (dimGrid->x>(unsigned)devprop.maxGridSize[X_AXIS]) {
    222:     *dimGrid_xcount = (dimGrid->x+(devprop.maxGridSize[X_AXIS]-1))/devprop.maxGridSize[X_AXIS];
    223:     dimGrid->x = devprop.maxGridSize[X_AXIS];
    224:   }
    225:   PetscFunctionReturn(PETSC_SUCCESS);
    226: }

    228: static inline PetscErrorCode SlepcKernelSetGrid2DTiles(PetscInt rows,PetscInt cols,dim3 *dimGrid,dim3 *dimBlock,PetscInt *dimGrid_xcount,PetscInt *dimGrid_ycount)
    229: {
    230:   int card;
    231: #if defined(PETSC_HAVE_CUDA)
    232:   struct cudaDeviceProp devprop;
    233: #elif defined(PETSC_HAVE_HIP)
    234:   hipDeviceProp_t devprop;
    235: #endif

    237:   PetscFunctionBegin;
    238: #if defined(PETSC_HAVE_CUDA)
    239:   PetscCallCUDA(cudaGetDevice(&card));
    240:   PetscCallCUDA(cudaGetDeviceProperties(&devprop,card));
    241: #elif defined(PETSC_HAVE_HIP)
    242:   PetscCallHIP(hipGetDevice(&card));
    243:   PetscCallHIP(hipGetDeviceProperties(&devprop,card));
    244: #endif
    245:   *dimGrid_xcount = *dimGrid_ycount = 1;

    247:   /* X axis */
    248:   dimGrid->x  = 1;
    249:   dimBlock->x = SLEPC_BLOCK_SIZE_X;
    250:   if (rows>SLEPC_BLOCK_SIZE_X*SLEPC_TILE_SIZE_X) dimGrid->x = (rows+SLEPC_BLOCK_SIZE_X*SLEPC_TILE_SIZE_X-1)/(SLEPC_BLOCK_SIZE_X*SLEPC_TILE_SIZE_X);
    251:   else dimBlock->x = (rows+SLEPC_TILE_SIZE_X-1)/SLEPC_TILE_SIZE_X;
    252:   if (dimGrid->x>(unsigned)devprop.maxGridSize[X_AXIS]) {
    253:     *dimGrid_xcount = (dimGrid->x+(devprop.maxGridSize[X_AXIS]-1))/devprop.maxGridSize[X_AXIS];
    254:     dimGrid->x = devprop.maxGridSize[X_AXIS];
    255:   }

    257:   /* Y axis */
    258:   dimGrid->y  = 1;
    259:   dimBlock->y = SLEPC_BLOCK_SIZE_Y;
    260:   if (cols>SLEPC_BLOCK_SIZE_Y*SLEPC_TILE_SIZE_Y) dimGrid->y = (cols+SLEPC_BLOCK_SIZE_Y*SLEPC_TILE_SIZE_Y-1)/(SLEPC_BLOCK_SIZE_Y*SLEPC_TILE_SIZE_Y);
    261:   else dimBlock->y = (cols+SLEPC_TILE_SIZE_Y-1)/SLEPC_TILE_SIZE_Y;
    262:   if (dimGrid->y>(unsigned)devprop.maxGridSize[Y_AXIS]) {
    263:     *dimGrid_ycount = (dimGrid->y+(devprop.maxGridSize[Y_AXIS]-1))/devprop.maxGridSize[Y_AXIS];
    264:     dimGrid->y = devprop.maxGridSize[Y_AXIS];
    265:   }
    266:   PetscFunctionReturn(PETSC_SUCCESS);
    267: }
    268: #undef X_AXIS
    269: #undef Y_AXIS
    270: #endif