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             : #pragma once
      12             : 
      13             : #include <slepcfn.h>
      14             : #include <slepc/private/slepcimpl.h>
      15             : 
      16             : /* SUBMANSEC = FN */
      17             : 
      18             : SLEPC_EXTERN PetscBool FNRegisterAllCalled;
      19             : SLEPC_EXTERN PetscErrorCode FNRegisterAll(void);
      20             : SLEPC_EXTERN PetscLogEvent FN_Evaluate;
      21             : 
      22             : typedef struct _FNOps *FNOps;
      23             : 
      24             : struct _FNOps {
      25             :   PetscErrorCode (*evaluatefunction)(FN,PetscScalar,PetscScalar*);
      26             :   PetscErrorCode (*evaluatederivative)(FN,PetscScalar,PetscScalar*);
      27             :   PetscErrorCode (*evaluatefunctionmat[FN_MAX_SOLVE])(FN,Mat,Mat);
      28             :   PetscErrorCode (*evaluatefunctionmatcuda[FN_MAX_SOLVE])(FN,Mat,Mat);
      29             :   PetscErrorCode (*evaluatefunctionmatvec[FN_MAX_SOLVE])(FN,Mat,Vec);
      30             :   PetscErrorCode (*evaluatefunctionmatveccuda[FN_MAX_SOLVE])(FN,Mat,Vec);
      31             :   PetscErrorCode (*setfromoptions)(FN,PetscOptionItems*);
      32             :   PetscErrorCode (*view)(FN,PetscViewer);
      33             :   PetscErrorCode (*duplicate)(FN,MPI_Comm,FN*);
      34             :   PetscErrorCode (*destroy)(FN);
      35             : };
      36             : 
      37             : #define FN_MAX_W 6
      38             : 
      39             : struct _p_FN {
      40             :   PETSCHEADER(struct _FNOps);
      41             :   /*------------------------- User parameters --------------------------*/
      42             :   PetscScalar    alpha;          /* inner scaling (argument) */
      43             :   PetscScalar    beta;           /* outer scaling (result) */
      44             :   PetscInt       method;         /* the method to compute matrix functions */
      45             :   FNParallelType pmode;          /* parallel mode (redundant or synchronized) */
      46             : 
      47             :   /*---------------------- Cached data and workspace -------------------*/
      48             :   Mat            W[FN_MAX_W];    /* workspace matrices */
      49             :   PetscInt       nw;             /* number of allocated W matrices */
      50             :   PetscInt       cw;             /* current W matrix */
      51             :   void           *data;
      52             : };
      53             : 
      54             : /*
      55             :   FN_AllocateWorkMat - Allocate a work Mat of the same dimension of A and copy
      56             :   its contents. The work matrix is returned in M and should be freed with
      57             :   FN_FreeWorkMat().
      58             : */
      59         894 : static inline PetscErrorCode FN_AllocateWorkMat(FN fn,Mat A,Mat *M)
      60             : {
      61         894 :   PetscInt       n,na;
      62         894 :   PetscBool      create=PETSC_FALSE;
      63             : 
      64         894 :   PetscFunctionBegin;
      65         894 :   *M = NULL;
      66         894 :   PetscCheck(fn->cw<FN_MAX_W,PETSC_COMM_SELF,PETSC_ERR_SUP,"Too many requested work matrices %" PetscInt_FMT,fn->cw);
      67         894 :   if (fn->nw<=fn->cw) {
      68         129 :     create=PETSC_TRUE;
      69         129 :     fn->nw++;
      70             :   } else {
      71         765 :     PetscCall(MatGetSize(fn->W[fn->cw],&n,NULL));
      72         765 :     PetscCall(MatGetSize(A,&na,NULL));
      73         765 :     if (n!=na) {
      74         407 :       PetscCall(MatDestroy(&fn->W[fn->cw]));
      75             :       create=PETSC_TRUE;
      76             :     }
      77             :   }
      78         536 :   if (create) PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&fn->W[fn->cw]));
      79         358 :   else PetscCall(MatCopy(A,fn->W[fn->cw],SAME_NONZERO_PATTERN));
      80         894 :   *M = fn->W[fn->cw];
      81         894 :   fn->cw++;
      82         894 :   PetscFunctionReturn(PETSC_SUCCESS);
      83             : }
      84             : 
      85             : /*
      86             :   FN_FreeWorkMat - Release a work matrix created with FN_AllocateWorkMat().
      87             : */
      88         894 : static inline PetscErrorCode FN_FreeWorkMat(FN fn,Mat *M)
      89             : {
      90         894 :   PetscFunctionBegin;
      91         894 :   PetscCheck(fn->cw,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"There are no work matrices");
      92         894 :   fn->cw--;
      93         894 :   PetscCheck(fn->W[fn->cw]==*M,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Work matrices must be freed in the reverse order of their creation");
      94         894 :   *M = NULL;
      95         894 :   PetscFunctionReturn(PETSC_SUCCESS);
      96             : }
      97             : 
      98             : SLEPC_INTERN PetscErrorCode FNSqrtmSchur(FN,PetscBLASInt,PetscScalar*,PetscBLASInt,PetscBool);
      99             : SLEPC_INTERN PetscErrorCode FNSqrtmDenmanBeavers(FN,PetscBLASInt,PetscScalar*,PetscBLASInt,PetscBool);
     100             : SLEPC_INTERN PetscErrorCode FNSqrtmNewtonSchulz(FN,PetscBLASInt,PetscScalar*,PetscBLASInt,PetscBool);
     101             : SLEPC_INTERN PetscErrorCode FNSqrtmSadeghi(FN,PetscBLASInt,PetscScalar*,PetscBLASInt);
     102             : SLEPC_INTERN PetscErrorCode SlepcNormAm(PetscBLASInt,PetscScalar*,PetscInt,PetscScalar*,PetscRandom,PetscReal*);
     103             : SLEPC_INTERN PetscErrorCode FNEvaluateFunctionMat_Private(FN,Mat,Mat,PetscBool);
     104             : SLEPC_INTERN PetscErrorCode FNEvaluateFunctionMatVec_Private(FN,Mat,Vec,PetscBool);
     105             : SLEPC_INTERN PetscErrorCode FNEvaluateFunctionMat_Exp_Higham(FN,Mat,Mat); /* used in FNPHI */
     106             : #if defined(PETSC_HAVE_CUDA)
     107             : SLEPC_INTERN PetscErrorCode FNSqrtmDenmanBeavers_CUDAm(FN,PetscBLASInt,PetscScalar*,PetscBLASInt,PetscBool);
     108             : SLEPC_INTERN PetscErrorCode FNSqrtmNewtonSchulz_CUDA(FN,PetscBLASInt,PetscScalar*,PetscBLASInt,PetscBool);
     109             : SLEPC_INTERN PetscErrorCode FNSqrtmSadeghi_CUDAm(FN,PetscBLASInt,PetscScalar*,PetscBLASInt);
     110             : #endif