scalapack.templates.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dftfe_scalapack_templates_h
#define dftfe_scalapack_templates_h
 
 
#include "headers.h"
 
 
#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
#  include <cfenv>
#endif
 
// useful examples:
// https://stackoverflow.com/questions/14147705/cholesky-decomposition-scalapack-error/14203864
// http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=139   // second post by
// Julien Langou
// https://andyspiros.wordpress.com/2011/07/08/an-example-of-blacs-with-c/
// http://qboxcode.org/trac/browser/qb/tags/rel1_63_4/src/Matrix.C
// https://gitlab.phys.ethz.ch/lwossnig/lecture/blob/a534f562dfb2ad5c564abe5c2356d5d956fb7218/examples/mpi/scalapack.cpp
// https://github.com/elemental/Elemental/blob/master/src/core/imports/scalapack.cpp
// https://scicomp.stackexchange.com/questions/7766/performance-optimization-or-tuning-possible-for-scalapack-gemm
//
// info:
// http://www.netlib.org/scalapack/slug/index.html       // User guide
// http://www.netlib.org/scalapack/slug/node135.html // How to Measure Errors
// NOTE: Non-templated functions are chosen over templated functions if their
// names match and template function is not explicitly called
namespace dftfe
{
  extern "C"
  {
    /* Basic Linear Algebra Communication Subprograms (BLACS) declarations */
    // https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dinitb.htm#dinitb
 
    /**
     * Determine how many processes are available and the current process rank.
     *
     * https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbpnf.htm
     */
    void
    Cblacs_pinfo(int *rank, int *nprocs);
 
    /**
     * Return internal BLACS value in @p val based on the input @p what and @p icontxt.
     * The most common use is in retrieving a default system context (@p what = 0, @p icontxt is ignored)
     * to be used in BLACS_GRIDINIT or BLACS_GRIDMAP.
     *
     * https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbget.htm
     */
    void
    Cblacs_get(int icontxt, int what, int *val);
 
    /**
     * Map the processes sequentially in row-major or column-major order
     * into the process grid. Input arguments must be the same on every process.
     *
     * On return, @p context is the integer handle to the BLACS context,
     * whereas on entry it is a system context to be used in creating the
     * BLACS context.
     *
     * https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbint.htm
     */
    void
    Cblacs_gridinit(int        *context,
                    const char *order,
                    int         grid_height,
                    int         grid_width);
 
    /**
     * Return the process row and column index.
     *
     * https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dbinfo.htm
     */
    void
    Cblacs_gridinfo(int  context,
                    int *grid_height,
                    int *grid_width,
                    int *grid_row,
                    int *grid_col);
 
    /**
     * Given the system process number, return the row and column coordinates in
     * the BLACS' process grid.
     */
    void
    Cblacs_pcoord(int ictxt, int pnum, int *prow, int *pcol);
 
    /**
     * Release a BLACS context.
     */
    void
    Cblacs_gridexit(int context);
 
    /**
     * This routines holds up execution of all processes within the indicated
     * scope until they have all called the routine.
     */
    void
    Cblacs_barrier(int, const char *);
 
    /**
     * Free all BLACS contexts and releases all allocated memory.
     */
    void
    Cblacs_exit(int error_code);
 
    /**
     * Receives a message from a process @prsrc, @p csrc into a general rectangular matrix.
     *
     * https://software.intel.com/en-us/mkl-developer-reference-c-gerv2d
     */
    void
    Cdgerv2d(int context, int M, int N, double *A, int lda, int rsrc, int csrc);
    void
    Csgerv2d(int context, int M, int N, float *A, int lda, int rsrc, int csrc);
 
    /**
     * Sends the general rectangular matrix A to the destination
     * process @p rdest @p cdest in the process grid.
     *
     * https://software.intel.com/en-us/mkl-developer-reference-c-2018-beta-gesd2d
     */
    void
    Cdgesd2d(int     context,
             int     M,
             int     N,
             double *A,
             int     lda,
             int     rdest,
             int     cdest);
    void
    Csgesd2d(int    context,
             int    M,
             int    N,
             float *A,
             int    lda,
             int    rdest,
             int    cdest);
 
    /**
     * Get BLACS context from MPI @p comm.
     */
    int
    Csys2blacs_handle(MPI_Comm comm);
 
    /**
     * Compute how many rows and columns each process owns (NUMber of Rows Or
     * Columns).
     *
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_dnumy.htm
     */
    int
    numroc_(const int *n,
            const int *nb,
            const int *iproc,
            const int *isproc,
            const int *nprocs);
    /**
     * Compute complex conjugate
     */
    void
    pzlacgv_(const int            *N,
             std::complex<double> *A,
             const int            *IA,
             const int            *JA,
             const int            *DESCA,
             const int            *INCX);
 
    /**
     * Compute the Cholesky factorization of an N-by-N real
     * symmetric positive definite distributed matrix sub( A ) denoting
     * A(IA:IA+N-1, JA:JA+N-1).
     *
     * http://www.netlib.org/scalapack/explore-html/d5/d9e/pdpotrf_8f_source.html
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lpotrf.htm
     */
    void
    pdpotrf_(const char *UPLO,
             const int  *N,
             double     *A,
             const int  *IA,
             const int  *JA,
             const int  *DESCA,
             int        *INFO);
    void
    pspotrf_(const char *UPLO,
             const int  *N,
             float      *A,
             const int  *IA,
             const int  *JA,
             const int  *DESCA,
             int        *INFO);
    void
    pzpotrf_(const char           *UPLO,
             const int            *N,
             std::complex<double> *A,
             const int            *IA,
             const int            *JA,
             const int            *DESCA,
             int                  *INFO);
 
    /**
     * Computes an LU factorization of a general distributed matrix sub( A )
     * using partial pivoting with row interchanges.
     *
     * http://www.netlib.org/scalapack/explore-html/df/dfe/pdgetrf_8f_source.html
     * https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgetrf.htm
     */
    void
    pdgetrf_(const int *m,
             const int *n,
             double    *A,
             const int *IA,
             const int *JA,
             const int *DESCA,
             int       *ipiv,
             int       *INFO);
    void
    psgetrf_(const int *m,
             const int *n,
             float     *A,
             const int *IA,
             const int *JA,
             const int *DESCA,
             int       *ipiv,
             int       *INFO);
    void
    pzgetrf_(const int            *m,
             const int            *n,
             std::complex<double> *A,
             const int            *IA,
             const int            *JA,
             const int            *DESCA,
             int                  *ipiv,
             int                  *INFO);
    /**
     * Compute the inverse of a real symmetric positive definite
     * distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1) using the
     * Cholesky factorization sub( A ) = U**T*U or L*L**T computed by
     * PDPOTRF.
     *
     * http://www.netlib.org/scalapack/explore-html/d2/d44/pdpotri_8f_source.html
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lpotri.htm
     * https://software.intel.com/en-us/mkl-developer-reference-c-p-potri
     */
    void
    pdpotri_(const char *UPLO,
             const int  *N,
             double     *A,
             const int  *IA,
             const int  *JA,
             const int  *DESCA,
             int        *INFO);
    void
    pspotri_(const char *UPLO,
             const int  *N,
             float      *A,
             const int  *IA,
             const int  *JA,
             const int  *DESCA,
             int        *INFO);
    void
    pzpotri_(const char           *UPLO,
             const int            *N,
             std::complex<double> *A,
             const int            *IA,
             const int            *JA,
             const int            *DESCA,
             int                  *INFO);
 
    /**
     * PDGETRI computes the inverse of a distributed matrix using the LU
     * factorization computed by PDGETRF. This method inverts U and then
     * computes the inverse of sub( A ) = A(IA:IA+N-1,JA:JA+N-1) denoted
     * InvA by solving the system InvA*L = inv(U) for InvA.
     *
     * http://www.netlib.org/scalapack/explore-html/d3/df3/pdgetri_8f_source.html
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgetri.htm
     */
    void
    pdgetri_(const int *N,
             double    *A,
             const int *IA,
             const int *JA,
             const int *DESCA,
             const int *ipiv,
             double    *work,
             int       *lwork,
             int       *iwork,
             int       *liwork,
             int       *info);
    void
    psgetri_(const int *N,
             float     *A,
             const int *IA,
             const int *JA,
             const int *DESCA,
             const int *ipiv,
             float     *work,
             int       *lwork,
             int       *iwork,
             int       *liwork,
             int       *info);
    void
    pzgetri_(const int            *N,
             std::complex<double> *A,
             const int            *IA,
             const int            *JA,
             const int            *DESCA,
             const int            *ipiv,
             std::complex<double> *work,
             int                  *lwork,
             int                  *iwork,
             int                  *liwork,
             int                  *info);
 
    /**
     * PDTRTRI computes the inverse of a upper or lower triangular
     * distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1).
     *
     * http://www.netlib.org/scalapack/explore-html/d9/dc0/pdtrtri_8f_source.html
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lpdtri.htm
     * https://software.intel.com/en-us/mkl-developer-reference-c-p-trtri
     */
    void
    pdtrtri_(const char *UPLO,
             const char *DIAG,
             const int  *N,
             double     *A,
             const int  *IA,
             const int  *JA,
             const int  *DESCA,
             int        *INFO);
    void
    pstrtri_(const char *UPLO,
             const char *DIAG,
             const int  *N,
             float      *A,
             const int  *IA,
             const int  *JA,
             const int  *DESCA,
             int        *INFO);
 
    void
    pztrtri_(const char           *UPLO,
             const char           *DIAG,
             const int            *N,
             std::complex<double> *A,
             const int            *IA,
             const int            *JA,
             const int            *DESCA,
             int                  *INFO);
 
    /**
     * Estimate the reciprocal of the condition number (in the
     * l1-norm) of a real symmetric positive definite distributed matrix
     * using the Cholesky factorization.
     *
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lpocon.htm#lpocon
     * http://www.netlib.org/scalapack/explore-html/d4/df7/pdpocon_8f.html
     * https://software.intel.com/en-us/mkl-developer-reference-fortran-pocon
     */
    void
    pdpocon_(const char   *uplo,
             const int    *N,
             const double *A,
             const int    *IA,
             const int    *JA,
             const int    *DESCA,
             const double *ANORM,
             double       *RCOND,
             double       *WORK,
             const int    *LWORK,
             int          *IWORK,
             const int    *LIWORK,
             int          *INFO);
    void
    pspocon_(const char  *uplo,
             const int   *N,
             const float *A,
             const int   *IA,
             const int   *JA,
             const int   *DESCA,
             const float *ANORM,
             float       *RCOND,
             float       *WORK,
             const int   *LWORK,
             int         *IWORK,
             const int   *LIWORK,
             int         *INFO);
 
    /**
     * Norm of a real symmetric matrix
     *
     * http://www.netlib.org/scalapack/explore-html/dd/d12/pdlansy_8f_source.html
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_pdlansy.htm#pdlansy
     */
    double
    pdlansy_(const char   *norm,
             const char   *uplo,
             const int    *N,
             const double *A,
             const int    *IA,
             const int    *JA,
             const int    *DESCA,
             double       *work);
    float
    pslansy_(const char  *norm,
             const char  *uplo,
             const int   *N,
             const float *A,
             const int   *IA,
             const int   *JA,
             const int   *DESCA,
             float       *work);
 
    /**
     * Compute the Least Common Multiple (LCM) of two positive integers @p M and @p N.
     * In fact the routine Compute the greatest common divisor (GCD) and
     * use the fact that M*N = GCD*LCM.
     *
     * http://www.netlib.org/scalapack/explore-html/d0/d9b/ilcm_8f_source.html
     */
    int
    ilcm_(const int *M, const int *N);
 
    /**
     * Return the ceiling of the division of two integers.
     *
     * http://www.netlib.org/scalapack/explore-html/df/d07/iceil_8f_source.html
     */
    int
    iceil_(const int *i1, const int *i2);
 
    /**
     * Initialize the descriptor vector with the 8 input arguments
     */
    void
    descinit_(int       *desc,
              const int *m,
              const int *n,
              const int *mb,
              const int *nb,
              const int *irsrc,
              const int *icsrc,
              const int *ictxt,
              const int *lld,
              int       *info);
 
    /**
     * Compute the global index of a distributed matrix entry
     * pointed to by the local index @p indxloc of the process indicated by
     * @p iproc.
     *
     * @param indxloc The local index of the distributed matrix entry.
     * @param nb Block size, size of the blocks the distributed matrix is split
     * into.
     * @param iproc The coordinate of the process whose local array row or column
     * is to be determined
     * @param isrcproc  The coordinate of the process that possesses the first
     * row/column of the distributed matrix
     * @param nprocs The total number processes over which the distributed matrix
     * is distributed
     */
    int
    indxl2g_(const int *indxloc,
             const int *nb,
             const int *iproc,
             const int *isrcproc,
             const int *nprocs);
 
    /**
     * Compute the solution to a real system of linear equations
     */
    void
    pdgesv_(const int *n,
            const int *nrhs,
            double    *A,
            const int *ia,
            const int *ja,
            const int *desca,
            int       *ipiv,
            double    *B,
            const int *ib,
            const int *jb,
            const int *descb,
            int       *info);
    void
    psgesv_(const int *n,
            const int *nrhs,
            float     *A,
            const int *ia,
            const int *ja,
            const int *desca,
            int       *ipiv,
            float     *B,
            const int *ib,
            const int *jb,
            const int *descb,
            int       *info);
 
    /**
     * Perform one of the matrix-matrix operations:
     * @f{align*}{
     * \mathrm{sub}(C) &\dealcoloneq \alpha
     * op(\mathrm{sub}(A))op(\mathrm{sub}(B))
     *                            + \beta \mathrm{sub}(C), \\
     * \mathrm{sub}(C) &\dealcoloneq \alpha
     * op(\mathrm{sub}(A))op(\mathrm{sub}(B))
     *                            + beta sub(C),
     * @f}
     * where
     * $\mathrm{sub}(C)$ denotes C(IC:IC+M-1,JC:JC+N-1),  and, $op(X)$ is one of
     * $op(X) = X$ or $op(X) = X^T$.
     */
    void
    pdgemm_(const char   *transa,
            const char   *transb,
            const int    *m,
            const int    *n,
            const int    *k,
            const double *alpha,
            const double *A,
            const int    *IA,
            const int    *JA,
            const int    *DESCA,
            const double *B,
            const int    *IB,
            const int    *JB,
            const int    *DESCB,
            const double *beta,
            double       *C,
            const int    *IC,
            const int    *JC,
            const int    *DESCC);
    void
    psgemm_(const char  *transa,
            const char  *transb,
            const int   *m,
            const int   *n,
            const int   *k,
            const float *alpha,
            const float *A,
            const int   *IA,
            const int   *JA,
            const int   *DESCA,
            const float *B,
            const int   *IB,
            const int   *JB,
            const int   *DESCB,
            const float *beta,
            float       *C,
            const int   *IC,
            const int   *JC,
            const int   *DESCC);
 
    void
    pzgemm_(const char                 *transa,
            const char                 *transb,
            const int                  *m,
            const int                  *n,
            const int                  *k,
            const std::complex<double> *alpha,
            const std::complex<double> *A,
            const int                  *IA,
            const int                  *JA,
            const int                  *DESCA,
            const std::complex<double> *B,
            const int                  *IB,
            const int                  *JB,
            const int                  *DESCB,
            const std::complex<double> *beta,
            std::complex<double>       *C,
            const int                  *IC,
            const int                  *JC,
            const int                  *DESCC);
 
    /**
     * Return the value of the one norm, or the Frobenius norm, or the infinity
     * norm, or the element of largest absolute value of a distributed matrix
     */
    double
    pdlange_(char const   *norm,
             const int    *m,
             const int    *n,
             const double *A,
             const int    *ia,
             const int    *ja,
             const int    *desca,
             double       *work);
    float
    pslange_(const char  *norm,
             const int   *m,
             const int   *n,
             const float *A,
             const int   *ia,
             const int   *ja,
             const int   *desca,
             float       *work);
 
    /**
     * Compute the process coordinate which possesses the entry of a
     * distributed matrix specified by a global index
     */
    int
    indxg2p_(const int *glob,
             const int *nb,
             const int *iproc,
             const int *isproc,
             const int *nprocs);
 
    /**
     * Compute all eigenvalues and, optionally, eigenvectors of a real symmetric
     * matrix A by calling the recommended sequence of ScaLAPACK routines. In
     * its present form, the routine assumes a homogeneous system and makes no
     * checks for consistency of the eigenvalues or eigenvectors across the
     * different processes. Because of this, it is possible that a heterogeneous
     * system may return incorrect results without any error messages.
     *
     * http://www.netlib.org/scalapack/explore-html/d0/d1a/pdsyev_8f.html
     * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lsyev.htm#lsyev
     */
    void
    pdsyev_(const char *jobz,
            const char *uplo,
            const int  *m,
            double     *A,
            const int  *ia,
            const int  *ja,
            int        *desca,
            double     *w,
            double     *z,
            const int  *iz,
            const int  *jz,
            int        *descz,
            double     *work,
            const int  *lwork,
            int        *info);
    void
    pssyev_(const char *jobz,
            const char *uplo,
            const int  *m,
            float      *A,
            const int  *ia,
            const int  *ja,
            int        *desca,
            float      *w,
            float      *z,
            const int  *iz,
            const int  *jz,
            int        *descz,
            float      *work,
            const int  *lwork,
            int        *info);
    void
    pzheev_(const char           *jobz,
            const char           *uplo,
            const int            *m,
            std::complex<double> *A,
            const int            *ia,
            const int            *ja,
            int                  *desca,
            double               *w,
            std::complex<double> *z,
            const int            *iz,
            const int            *jz,
            int                  *descz,
            std::complex<double> *work,
            const int            *lwork,
            int                  *info);
 
    /**
     * Copy all or a part of a distributed matrix A to another distributed
     * matrix B. No communication is performed, pdlacpy performs a local copy
     * $\mathrm{sub}(A) \dealcoloneq \mathrm{sub}(B)$, where $\mathrm{sub}(A)$
     * denotes $A(ia:ia+m-1, ja:ja+n-1)$ and $\mathrm{sub}(B)$ denotes
     * $B(ib:ib+m-1, jb:jb+n-1)$.
     */
    void
    pdlacpy_(const char   *uplo,
             const int    *m,
             const int    *n,
             const double *A,
             const int    *ia,
             const int    *ja,
             const int    *desca,
             double       *B,
             const int    *ib,
             const int    *jb,
             const int    *descb);
    void
    pslacpy_(const char  *uplo,
             const int   *m,
             const int   *n,
             const float *A,
             const int   *ia,
             const int   *ja,
             const int   *desca,
             float       *B,
             const int   *ib,
             const int   *jb,
             const int   *descb);
 
    /**
     * Copies the content of a general rectangular distributed matrix @p A to another distributed matrix @p B
     * It is not required that the matrices A and B have the same process grid
     * or block size, e.g. copying a matrix from a one-dimensional to a
     * two-dimensional process grid
     * @p ictxt is a context which is at least a union of all processes in
     * context A and B
     */
    void
    pdgemr2d_(const int    *m,
              const int    *n,
              const double *A,
              const int    *ia,
              const int    *ja,
              const int    *desca,
              double       *B,
              const int    *ib,
              const int    *jb,
              const int    *descb,
              const int    *ictxt);
    void
    psgemr2d_(const int   *m,
              const int   *n,
              const float *A,
              const int   *ia,
              const int   *ja,
              const int   *desca,
              float       *B,
              const int   *ib,
              const int   *jb,
              const int   *descb,
              const int   *ictxt);
 
    /**
     * helper routines determining machine precision
     */
    double
    pdlamch_(const int *ictxt, const char *cmach);
    float
    pslamch_(const int *ictxt, const char *cmach);
 
 
    /**
     *  psyevx computes selected eigenvalues and, optionally, eigenvectors
     *  of a real symmetric matrix A. Eigenvalues/vectors can be selected by
     *  specifying a range of values or a range of indices for the desired
     *  eigenvalues.
     */
    void
    pdsyevx_(const char   *jobz,
             const char   *range,
             const char   *uplo,
             const int    *n,
             double       *A,
             const int    *ia,
             const int    *ja,
             const int    *desca,
             const double *VL,
             const double *VU,
             const int    *il,
             const int    *iu,
             const double *abstol,
             const int    *m,
             const int    *nz,
             double       *w,
             double       *orfac,
             double       *Z,
             const int    *iz,
             const int    *jz,
             const int    *descz,
             double       *work,
             int          *lwork,
             int          *iwork,
             int          *liwork,
             int          *ifail,
             int          *iclustr,
             double       *gap,
             int          *info);
    void
    pssyevx_(const char  *jobz,
             const char  *range,
             const char  *uplo,
             const int   *n,
             float       *A,
             const int   *ia,
             const int   *ja,
             const int   *desca,
             const float *VL,
             const float *VU,
             const int   *il,
             const int   *iu,
             const float *abstol,
             const int   *m,
             const int   *nz,
             float       *w,
             float       *orfac,
             float       *Z,
             const int   *iz,
             const int   *jz,
             const int   *descz,
             float       *work,
             int         *lwork,
             int         *iwork,
             int         *liwork,
             int         *ifail,
             int         *iclustr,
             float       *gap,
             int         *info);
    void
    pzheevx_(const char           *jobz,
             const char           *range,
             const char           *uplo,
             const int            *n,
             std::complex<double> *A,
             const int            *ia,
             const int            *ja,
             const int            *desca,
             const double         *VL,
             const double         *VU,
             const int            *il,
             const int            *iu,
             const double         *abstol,
             const int            *m,
             const int            *nz,
             double               *w,
             double               *orfac,
             std::complex<double> *Z,
             const int            *iz,
             const int            *jz,
             const int            *descz,
             std::complex<double> *work,
             int                  *lwork,
             int                  *iwork,
             int                  *liwork,
             int                  *ifail,
             int                  *iclustr,
             double               *gap,
             int                  *info);
 
    /*
     * PDGESVD computes the singular value decomposition (SVD) of an
     * M-by-N matrix A, optionally computing the left and/or right
     * singular vectors
     */
    void
    pdgesvd_(const char *jobu,
             const char *jobvt,
             const int  *m,
             const int  *n,
             double     *A,
             const int  *ia,
             const int  *ja,
             const int  *desca,
             double     *S,
             double     *U,
             const int  *iu,
             const int  *ju,
             const int  *descu,
             double     *VT,
             const int  *ivt,
             const int  *jvt,
             const int  *descvt,
             double     *work,
             int        *lwork,
             int        *info);
    void
    psgesvd_(const char *jobu,
             const char *jobvt,
             const int  *m,
             const int  *n,
             float      *A,
             const int  *ia,
             const int  *ja,
             const int  *desca,
             float      *S,
             float      *U,
             const int  *iu,
             const int  *ju,
             const int  *descu,
             float      *VT,
             const int  *ivt,
             const int  *jvt,
             const int  *descvt,
             float      *work,
             int        *lwork,
             int        *info);
 
    /*
     * P_GELS solves overdetermined or underdetermined real linear
     * systems involving an M-by-N matrix A, or its transpose,
     * using a QR or LQ factorization of A.  It is assumed that A has full rank.
     */
    void
    pdgels_(const char *trans,
            const int  *m,
            const int  *n,
            const int  *nrhs,
            double     *A,
            const int  *ia,
            const int  *ja,
            const int  *desca,
            double     *B,
            const int  *ib,
            const int  *jb,
            const int  *descb,
            double     *work,
            int        *lwork,
            int        *info);
    void
    psgels_(const char *trans,
            const int  *m,
            const int  *n,
            const int  *nrhs,
            float      *A,
            const int  *ia,
            const int  *ja,
            const int  *desca,
            float      *B,
            const int  *ib,
            const int  *jb,
            const int  *descb,
            float      *work,
            int        *lwork,
            int        *info);
 
    /*
     * Perform matrix sum:
     * @f{equation*}
     * C \dealcoloneq \beta C + \alpha op(A),
     * @f
     * where $op(A)$ denotes either $op(A) = A$ or $op(A)=A^T$.
     */
    void
    pdgeadd_(const char   *transa,
             const int    *m,
             const int    *n,
             const double *alpha,
             const double *A,
             const int    *IA,
             const int    *JA,
             const int    *DESCA,
             const double *beta,
             double       *C,
             const int    *IC,
             const int    *JC,
             const int    *DESCC);
    void
    psgeadd_(const char  *transa,
             const int   *m,
             const int   *n,
             const float *alpha,
             const float *A,
             const int   *IA,
             const int   *JA,
             const int   *DESCA,
             const float *beta,
             float       *C,
             const int   *IC,
             const int   *JC,
             const int   *DESCC);
    void
    pzgeadd_(const char                 *transa,
             const int                  *m,
             const int                  *n,
             const std::complex<double> *alpha,
             const std::complex<double> *A,
             const int                  *IA,
             const int                  *JA,
             const int                  *DESCA,
             const std::complex<double> *beta,
             std::complex<double>       *C,
             const int                  *IC,
             const int                  *JC,
             const int                  *DESCC);
 
    /**
     * Routine to transpose a matrix:
     * C = beta C + alpha A^T
     */
    void
    pdtran_(const int    *m,
            const int    *n,
            const double *alpha,
            const double *A,
            const int    *IA,
            const int    *JA,
            const int    *DESCA,
            const double *beta,
            double       *C,
            const int    *IC,
            const int    *JC,
            const int    *DESCC);
    void
    pstran_(const int   *m,
            const int   *n,
            const float *alpha,
            const float *A,
            const int   *IA,
            const int   *JA,
            const int   *DESCA,
            const float *beta,
            float       *C,
            const int   *IC,
            const int   *JC,
            const int   *DESCC);
 
    /**
     *  psyevr computes selected eigenvalues and, optionally, eigenvectors
     *  of a real symmetric matrix A using a parallel implementation of the MRR
     * algorithm. Eigenvalues/vectors can be selected by specifying a range of
     * values or a range of indices for the desired eigenvalues.
     */
    void
    pdsyevr_(const char   *jobz,
             const char   *range,
             const char   *uplo,
             const int    *n,
             double       *A,
             const int    *IA,
             const int    *JA,
             const int    *DESCA,
             const double *VL,
             const double *VU,
             const int    *IL,
             const int    *IU,
             int          *m,
             int          *nz,
             double       *w,
             double       *Z,
             const int    *IZ,
             const int    *JZ,
             const int    *DESCZ,
             double       *work,
             int          *lwork,
             int          *iwork,
             int          *liwork,
             int          *info);
    void
    pssyevr_(const char  *jobz,
             const char  *range,
             const char  *uplo,
             const int   *n,
             float       *A,
             const int   *IA,
             const int   *JA,
             const int   *DESCA,
             const float *VL,
             const float *VU,
             const int   *IL,
             const int   *IU,
             int         *m,
             int         *nz,
             float       *w,
             float       *Z,
             const int   *IZ,
             const int   *JZ,
             const int   *DESCZ,
             float       *work,
             int         *lwork,
             int         *iwork,
             int         *liwork,
             int         *info);
    void
    pzheevr_(const char           *jobz,
             const char           *range,
             const char           *uplo,
             const int            *n,
             std::complex<double> *A,
             const int            *IA,
             const int            *JA,
             const int            *DESCA,
             const double         *VL,
             const double         *VU,
             const int            *IL,
             const int            *IU,
             int                  *m,
             int                  *nz,
             double               *w,
             std::complex<double> *Z,
             const int            *IZ,
             const int            *JZ,
             const int            *DESCZ,
             std::complex<double> *work,
             int                  *lwork,
             double               *rwork,
             int                  *lrwork,
             int                  *iwork,
             int                  *liwork,
             int                  *info);
  }
 
 
 
  /*
   * In the following we have template wrappers for the ScaLAPACK routines
   * wrappers for other numeric types can be added in the future
   */
  template <typename number>
  inline void
  Cgerv2d(int /*context*/,
          int /*M*/,
          int /*N*/,
          number * /*A*/,
          int /*lda*/,
          int /*rsrc*/,
          int /*csrc*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  Cgerv2d(int context, int M, int N, double *A, int lda, int rsrc, int csrc)
  {
    Cdgerv2d(context, M, N, A, lda, rsrc, csrc);
  }
 
  inline void
  Cgerv2d(int context, int M, int N, float *A, int lda, int rsrc, int csrc)
  {
    Csgerv2d(context, M, N, A, lda, rsrc, csrc);
  }
 
 
  template <typename number>
  inline void
  Cgesd2d(int /*context*/,
          int /*M*/,
          int /*N*/,
          number * /*A*/,
          int /*lda*/,
          int /*rdest*/,
          int /*cdest*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  Cgesd2d(int context, int M, int N, double *A, int lda, int rdest, int cdest)
  {
    Cdgesd2d(context, M, N, A, lda, rdest, cdest);
  }
 
  inline void
  Cgesd2d(int context, int M, int N, float *A, int lda, int rdest, int cdest)
  {
    Csgesd2d(context, M, N, A, lda, rdest, cdest);
  }
 
  inline void
  pplacgv(const int *N,
          double    *A,
          const int *IA,
          const int *JA,
          const int *DESCA,
          const int *INCX)
  {}
 
  inline void
  pplacgv(const int            *N,
          std::complex<double> *A,
          const int            *IA,
          const int            *JA,
          const int            *DESCA,
          const int            *INCX)
  {
    pzlacgv_(N, A, IA, JA, DESCA, INCX);
  }
 
  template <typename number>
  inline void
  ppotrf(const char * /*UPLO*/,
         const int * /*N*/,
         number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         int * /*INFO*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  ppotrf(const char *UPLO,
         const int  *N,
         double     *A,
         const int  *IA,
         const int  *JA,
         const int  *DESCA,
         int        *INFO)
  {
    pdpotrf_(UPLO, N, A, IA, JA, DESCA, INFO);
  }
 
  inline void
  ppotrf(const char *UPLO,
         const int  *N,
         float      *A,
         const int  *IA,
         const int  *JA,
         const int  *DESCA,
         int        *INFO)
  {
    pspotrf_(UPLO, N, A, IA, JA, DESCA, INFO);
  }
 
  inline void
  ppotrf(const char           *UPLO,
         const int            *N,
         std::complex<double> *A,
         const int            *IA,
         const int            *JA,
         const int            *DESCA,
         int                  *INFO)
  {
    pzpotrf_(UPLO, N, A, IA, JA, DESCA, INFO);
  }
 
  template <typename number>
  inline void
  pgetrf(const int * /*m*/,
         const int * /*n*/,
         number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         int * /*ipiv*/,
         int * /*INFO*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgetrf(const int *m,
         const int *n,
         double    *A,
         const int *IA,
         const int *JA,
         const int *DESCA,
         int       *ipiv,
         int       *INFO)
  {
    pdgetrf_(m, n, A, IA, JA, DESCA, ipiv, INFO);
  }
 
  inline void
  pgetrf(const int *m,
         const int *n,
         float     *A,
         const int *IA,
         const int *JA,
         const int *DESCA,
         int       *ipiv,
         int       *INFO)
  {
    psgetrf_(m, n, A, IA, JA, DESCA, ipiv, INFO);
  }
 
  inline void
  pgetrf(const int            *m,
         const int            *n,
         std::complex<double> *A,
         const int            *IA,
         const int            *JA,
         const int            *DESCA,
         int                  *ipiv,
         int                  *INFO)
  {
    pzgetrf_(m, n, A, IA, JA, DESCA, ipiv, INFO);
  }
 
  template <typename number>
  inline void
  ppotri(const char * /*UPLO*/,
         const int * /*N*/,
         number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         int * /*INFO*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  ppotri(const char *UPLO,
         const int  *N,
         double     *A,
         const int  *IA,
         const int  *JA,
         const int  *DESCA,
         int        *INFO)
  {
    pdpotri_(UPLO, N, A, IA, JA, DESCA, INFO);
  }
 
  inline void
  ppotri(const char *UPLO,
         const int  *N,
         float      *A,
         const int  *IA,
         const int  *JA,
         const int  *DESCA,
         int        *INFO)
  {
    pspotri_(UPLO, N, A, IA, JA, DESCA, INFO);
  }
 
  inline void
  ppotri(const char           *UPLO,
         const int            *N,
         std::complex<double> *A,
         const int            *IA,
         const int            *JA,
         const int            *DESCA,
         int                  *INFO)
  {
    pzpotri_(UPLO, N, A, IA, JA, DESCA, INFO);
  }
 
  template <typename number>
  inline void
  pgetri(const int * /*N*/,
         number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         const int * /*ipiv*/,
         number * /*work*/,
         int * /*lwork*/,
         int * /*iwork*/,
         int * /*liwork*/,
         int * /*info*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgetri(const int *N,
         double    *A,
         const int *IA,
         const int *JA,
         const int *DESCA,
         const int *ipiv,
         double    *work,
         int       *lwork,
         int       *iwork,
         int       *liwork,
         int       *info)
  {
    pdgetri_(N, A, IA, JA, DESCA, ipiv, work, lwork, iwork, liwork, info);
  }
 
  inline void
  pgetri(const int *N,
         float     *A,
         const int *IA,
         const int *JA,
         const int *DESCA,
         const int *ipiv,
         float     *work,
         int       *lwork,
         int       *iwork,
         int       *liwork,
         int       *info)
  {
    psgetri_(N, A, IA, JA, DESCA, ipiv, work, lwork, iwork, liwork, info);
  }
 
  inline void
  pgetri(const int            *N,
         std::complex<double> *A,
         const int            *IA,
         const int            *JA,
         const int            *DESCA,
         const int            *ipiv,
         std::complex<double> *work,
         int                  *lwork,
         int                  *iwork,
         int                  *liwork,
         int                  *info)
  {
    pzgetri_(N, A, IA, JA, DESCA, ipiv, work, lwork, iwork, liwork, info);
  }
 
 
 
  template <typename number>
  inline void
  ptrtri(const char * /*UPLO*/,
         const char * /*DIAG*/,
         const int * /*N*/,
         number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         int * /*INFO*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  ptrtri(const char *UPLO,
         const char *DIAG,
         const int  *N,
         double     *A,
         const int  *IA,
         const int  *JA,
         const int  *DESCA,
         int        *INFO)
  {
    pdtrtri_(UPLO, DIAG, N, A, IA, JA, DESCA, INFO);
  }
 
  inline void
  ptrtri(const char *UPLO,
         const char *DIAG,
         const int  *N,
         float      *A,
         const int  *IA,
         const int  *JA,
         const int  *DESCA,
         int        *INFO)
  {
    pstrtri_(UPLO, DIAG, N, A, IA, JA, DESCA, INFO);
  }
 
  inline void
  ptrtri(const char           *UPLO,
         const char           *DIAG,
         const int            *N,
         std::complex<double> *A,
         const int            *IA,
         const int            *JA,
         const int            *DESCA,
         int                  *INFO)
  {
    pztrtri_(UPLO, DIAG, N, A, IA, JA, DESCA, INFO);
  }
 
  template <typename number>
  inline void
  ppocon(const char * /*uplo*/,
         const int * /*N*/,
         const number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         const number * /*ANORM*/,
         number * /*RCOND*/,
         number * /*WORK*/,
         const int * /*LWORK*/,
         int * /*IWORK*/,
         const int * /*LIWORK*/,
         int * /*INFO*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  ppocon(const char   *uplo,
         const int    *N,
         const double *A,
         const int    *IA,
         const int    *JA,
         const int    *DESCA,
         const double *ANORM,
         double       *RCOND,
         double       *WORK,
         const int    *LWORK,
         int          *IWORK,
         const int    *LIWORK,
         int          *INFO)
  {
    pdpocon_(uplo,
             N,
             A,
             IA,
             JA,
             DESCA,
             ANORM,
             RCOND,
             WORK,
             LWORK,
             IWORK,
             LIWORK,
             INFO);
  }
 
  inline void
  ppocon(const char  *uplo,
         const int   *N,
         const float *A,
         const int   *IA,
         const int   *JA,
         const int   *DESCA,
         const float *ANORM,
         float       *RCOND,
         float       *WORK,
         const int   *LWORK,
         int         *IWORK,
         const int   *LIWORK,
         int         *INFO)
  {
    pspocon_(uplo,
             N,
             A,
             IA,
             JA,
             DESCA,
             ANORM,
             RCOND,
             WORK,
             LWORK,
             IWORK,
             LIWORK,
             INFO);
  }
 
 
  template <typename number>
  inline number
  plansy(const char * /*norm*/,
         const char * /*uplo*/,
         const int * /*N*/,
         const number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         number * /*work*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline double
  plansy(const char   *norm,
         const char   *uplo,
         const int    *N,
         const double *A,
         const int    *IA,
         const int    *JA,
         const int    *DESCA,
         double       *work)
  {
    return pdlansy_(norm, uplo, N, A, IA, JA, DESCA, work);
  }
 
  inline float
  plansy(const char  *norm,
         const char  *uplo,
         const int   *N,
         const float *A,
         const int   *IA,
         const int   *JA,
         const int   *DESCA,
         float       *work)
  {
    return pslansy_(norm, uplo, N, A, IA, JA, DESCA, work);
  }
 
 
  template <typename number>
  inline void
  pgesv(const int * /*n*/,
        const int * /*nrhs*/,
        number * /*A*/,
        const int * /*ia*/,
        const int * /*ja*/,
        const int * /*desca*/,
        int * /*ipiv*/,
        number * /*B*/,
        const int * /*ib*/,
        const int * /*jb*/,
        const int * /*descb*/,
        int * /*info*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgesv(const int *n,
        const int *nrhs,
        double    *A,
        const int *ia,
        const int *ja,
        const int *desca,
        int       *ipiv,
        double    *B,
        const int *ib,
        const int *jb,
        const int *descb,
        int       *info)
  {
    pdgesv_(n, nrhs, A, ia, ja, desca, ipiv, B, ib, jb, descb, info);
  }
 
  inline void
  pgesv(const int *n,
        const int *nrhs,
        float     *A,
        const int *ia,
        const int *ja,
        const int *desca,
        int       *ipiv,
        float     *B,
        const int *ib,
        const int *jb,
        const int *descb,
        int       *info)
  {
    psgesv_(n, nrhs, A, ia, ja, desca, ipiv, B, ib, jb, descb, info);
  }
 
 
  template <typename number>
  inline void
  pgemm(const char * /*transa*/,
        const char * /*transb*/,
        const int * /*m*/,
        const int * /*n*/,
        const int * /*k*/,
        const number * /*alpha*/,
        const number * /*A*/,
        const int * /*IA*/,
        const int * /*JA*/,
        const int * /*DESCA*/,
        const number * /*B*/,
        const int * /*IB*/,
        const int * /*JB*/,
        const int * /*DESCB*/,
        const number * /*beta*/,
        number * /*C*/,
        const int * /*IC*/,
        const int * /*JC*/,
        const int * /*DESCC*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgemm(const char   *transa,
        const char   *transb,
        const int    *m,
        const int    *n,
        const int    *k,
        const double *alpha,
        const double *A,
        const int    *IA,
        const int    *JA,
        const int    *DESCA,
        const double *B,
        const int    *IB,
        const int    *JB,
        const int    *DESCB,
        const double *beta,
        double       *C,
        const int    *IC,
        const int    *JC,
        const int    *DESCC)
  {
    pdgemm_(transa,
            transb,
            m,
            n,
            k,
            alpha,
            A,
            IA,
            JA,
            DESCA,
            B,
            IB,
            JB,
            DESCB,
            beta,
            C,
            IC,
            JC,
            DESCC);
  }
 
  inline void
  pgemm(const char  *transa,
        const char  *transb,
        const int   *m,
        const int   *n,
        const int   *k,
        const float *alpha,
        const float *A,
        const int   *IA,
        const int   *JA,
        const int   *DESCA,
        const float *B,
        const int   *IB,
        const int   *JB,
        const int   *DESCB,
        const float *beta,
        float       *C,
        const int   *IC,
        const int   *JC,
        const int   *DESCC)
  {
    psgemm_(transa,
            transb,
            m,
            n,
            k,
            alpha,
            A,
            IA,
            JA,
            DESCA,
            B,
            IB,
            JB,
            DESCB,
            beta,
            C,
            IC,
            JC,
            DESCC);
  }
 
 
  inline void
  pgemm(const char                 *transa,
        const char                 *transb,
        const int                  *m,
        const int                  *n,
        const int                  *k,
        const std::complex<double> *alpha,
        const std::complex<double> *A,
        const int                  *IA,
        const int                  *JA,
        const int                  *DESCA,
        const std::complex<double> *B,
        const int                  *IB,
        const int                  *JB,
        const int                  *DESCB,
        const std::complex<double> *beta,
        std::complex<double>       *C,
        const int                  *IC,
        const int                  *JC,
        const int                  *DESCC)
  {
    pzgemm_(transa,
            transb,
            m,
            n,
            k,
            alpha,
            A,
            IA,
            JA,
            DESCA,
            B,
            IB,
            JB,
            DESCB,
            beta,
            C,
            IC,
            JC,
            DESCC);
  }
 
  template <typename number>
  inline number
  plange(const char * /*norm*/,
         const int * /*m*/,
         const int * /*n*/,
         const number * /*A*/,
         const int * /*ia*/,
         const int * /*ja*/,
         const int * /*desca*/,
         number * /*work*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline double
  plange(const char   *norm,
         const int    *m,
         const int    *n,
         const double *A,
         const int    *ia,
         const int    *ja,
         const int    *desca,
         double       *work)
  {
    return pdlange_(norm, m, n, A, ia, ja, desca, work);
  }
 
  inline float
  plange(const char  *norm,
         const int   *m,
         const int   *n,
         const float *A,
         const int   *ia,
         const int   *ja,
         const int   *desca,
         float       *work)
  {
    return pslange_(norm, m, n, A, ia, ja, desca, work);
  }
 
 
  template <typename number>
  inline void
  psyev(const char * /*jobz*/,
        const char * /*uplo*/,
        const int * /*m*/,
        number * /*A*/,
        const int * /*ia*/,
        const int * /*ja*/,
        int * /*desca*/,
        number * /*w*/,
        number * /*z*/,
        const int * /*iz*/,
        const int * /*jz*/,
        int * /*descz*/,
        number * /*work*/,
        const int * /*lwork*/,
        int * /*info*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  psyev(const char *jobz,
        const char *uplo,
        const int  *m,
        double     *A,
        const int  *ia,
        const int  *ja,
        int        *desca,
        double     *w,
        double     *z,
        const int  *iz,
        const int  *jz,
        int        *descz,
        double     *work,
        const int  *lwork,
        int        *info)
  {
    pdsyev_(
      jobz, uplo, m, A, ia, ja, desca, w, z, iz, jz, descz, work, lwork, info);
  }
 
  inline void
  psyev(const char *jobz,
        const char *uplo,
        const int  *m,
        float      *A,
        const int  *ia,
        const int  *ja,
        int        *desca,
        float      *w,
        float      *z,
        const int  *iz,
        const int  *jz,
        int        *descz,
        float      *work,
        const int  *lwork,
        int        *info)
  {
    pssyev_(
      jobz, uplo, m, A, ia, ja, desca, w, z, iz, jz, descz, work, lwork, info);
  }
 
  inline void
  psyev(const char           *jobz,
        const char           *uplo,
        const int            *m,
        std::complex<double> *A,
        const int            *ia,
        const int            *ja,
        int                  *desca,
        double               *w,
        std::complex<double> *z,
        const int            *iz,
        const int            *jz,
        int                  *descz,
        std::complex<double> *work,
        const int            *lwork,
        int                  *info)
  {
    pzheev_(
      jobz, uplo, m, A, ia, ja, desca, w, z, iz, jz, descz, work, lwork, info);
  }
 
  template <typename number>
  inline void
  placpy(const char * /*uplo*/,
         const int * /*m*/,
         const int * /*n*/,
         const number * /*A*/,
         const int * /*ia*/,
         const int * /*ja*/,
         const int * /*desca*/,
         number * /*B*/,
         const int * /*ib*/,
         const int * /*jb*/,
         const int * /*descb*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  placpy(const char   *uplo,
         const int    *m,
         const int    *n,
         const double *A,
         const int    *ia,
         const int    *ja,
         const int    *desca,
         double       *B,
         const int    *ib,
         const int    *jb,
         const int    *descb)
  {
    pdlacpy_(uplo, m, n, A, ia, ja, desca, B, ib, jb, descb);
  }
 
  inline void
  placpy(const char  *uplo,
         const int   *m,
         const int   *n,
         const float *A,
         const int   *ia,
         const int   *ja,
         const int   *desca,
         float       *B,
         const int   *ib,
         const int   *jb,
         const int   *descb)
  {
    pslacpy_(uplo, m, n, A, ia, ja, desca, B, ib, jb, descb);
  }
 
 
  template <typename number>
  inline void
  pgemr2d(const int * /*m*/,
          const int * /*n*/,
          const number * /*A*/,
          const int * /*ia*/,
          const int * /*ja*/,
          const int * /*desca*/,
          number * /*B*/,
          const int * /*ib*/,
          const int * /*jb*/,
          const int * /*descb*/,
          const int * /*ictxt*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgemr2d(const int    *m,
          const int    *n,
          const double *A,
          const int    *ia,
          const int    *ja,
          const int    *desca,
          double       *B,
          const int    *ib,
          const int    *jb,
          const int    *descb,
          const int    *ictxt)
  {
    pdgemr2d_(m, n, A, ia, ja, desca, B, ib, jb, descb, ictxt);
  }
 
  inline void
  pgemr2d(const int   *m,
          const int   *n,
          const float *A,
          const int   *ia,
          const int   *ja,
          const int   *desca,
          float       *B,
          const int   *ib,
          const int   *jb,
          const int   *descb,
          const int   *ictxt)
  {
    psgemr2d_(m, n, A, ia, ja, desca, B, ib, jb, descb, ictxt);
  }
 
 
  template <typename number>
  inline void
  plamch(const int * /*ictxt*/, const char * /*cmach*/, number & /*val*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  plamch(const int *ictxt, const char *cmach, double &val)
  {
    val = pdlamch_(ictxt, cmach);
  }
 
  inline void
  plamch(const int *ictxt, const char *cmach, float &val)
  {
    val = pslamch_(ictxt, cmach);
  }
 
 
  template <typename number>
  inline void
  psyevx(const char * /*jobz*/,
         const char * /*range*/,
         const char * /*uplo*/,
         const int * /*n*/,
         number * /*A*/,
         const int * /*ia*/,
         const int * /*ja*/,
         const int * /*desca*/,
         number * /*VL*/,
         number * /*VU*/,
         const int * /*il*/,
         const int * /*iu*/,
         number * /*abstol*/,
         const int * /*m*/,
         const int * /*nz*/,
         number * /*w*/,
         number * /*orfac*/,
         number * /*Z*/,
         const int * /*iz*/,
         const int * /*jz*/,
         const int * /*descz*/,
         number * /*work*/,
         int * /*lwork*/,
         int * /*iwork*/,
         int * /*liwork*/,
         int * /*ifail*/,
         int * /*iclustr*/,
         number * /*gap*/,
         int * /*info*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  psyevx(const char *jobz,
         const char *range,
         const char *uplo,
         const int  *n,
         double     *A,
         const int  *ia,
         const int  *ja,
         const int  *desca,
         double     *VL,
         double     *VU,
         const int  *il,
         const int  *iu,
         double     *abstol,
         const int  *m,
         const int  *nz,
         double     *w,
         double     *orfac,
         double     *Z,
         const int  *iz,
         const int  *jz,
         const int  *descz,
         double     *work,
         int        *lwork,
         int        *iwork,
         int        *liwork,
         int        *ifail,
         int        *iclustr,
         double     *gap,
         int        *info)
  {
    pdsyevx_(jobz,
             range,
             uplo,
             n,
             A,
             ia,
             ja,
             desca,
             VL,
             VU,
             il,
             iu,
             abstol,
             m,
             nz,
             w,
             orfac,
             Z,
             iz,
             jz,
             descz,
             work,
             lwork,
             iwork,
             liwork,
             ifail,
             iclustr,
             gap,
             info);
  }
 
  inline void
  psyevx(const char *jobz,
         const char *range,
         const char *uplo,
         const int  *n,
         float      *A,
         const int  *ia,
         const int  *ja,
         const int  *desca,
         float      *VL,
         float      *VU,
         const int  *il,
         const int  *iu,
         float      *abstol,
         const int  *m,
         const int  *nz,
         float      *w,
         float      *orfac,
         float      *Z,
         const int  *iz,
         const int  *jz,
         const int  *descz,
         float      *work,
         int        *lwork,
         int        *iwork,
         int        *liwork,
         int        *ifail,
         int        *iclustr,
         float      *gap,
         int        *info)
  {
    pssyevx_(jobz,
             range,
             uplo,
             n,
             A,
             ia,
             ja,
             desca,
             VL,
             VU,
             il,
             iu,
             abstol,
             m,
             nz,
             w,
             orfac,
             Z,
             iz,
             jz,
             descz,
             work,
             lwork,
             iwork,
             liwork,
             ifail,
             iclustr,
             gap,
             info);
  }
 
  inline void
  psyevx(const char           *jobz,
         const char           *range,
         const char           *uplo,
         const int            *n,
         std::complex<double> *A,
         const int            *ia,
         const int            *ja,
         const int            *desca,
         double               *VL,
         double               *VU,
         const int            *il,
         const int            *iu,
         double               *abstol,
         const int            *m,
         const int            *nz,
         double               *w,
         double               *orfac,
         std::complex<double> *Z,
         const int            *iz,
         const int            *jz,
         const int            *descz,
         std::complex<double> *work,
         int                  *lwork,
         int                  *iwork,
         int                  *liwork,
         int                  *ifail,
         int                  *iclustr,
         double               *gap,
         int                  *info)
  {
    pzheevx_(jobz,
             range,
             uplo,
             n,
             A,
             ia,
             ja,
             desca,
             VL,
             VU,
             il,
             iu,
             abstol,
             m,
             nz,
             w,
             orfac,
             Z,
             iz,
             jz,
             descz,
             work,
             lwork,
             iwork,
             liwork,
             ifail,
             iclustr,
             gap,
             info);
  }
 
  template <typename number>
  inline void
  pgesvd(const char * /*jobu*/,
         const char * /*jobvt*/,
         const int * /*m*/,
         const int * /*n*/,
         number * /*A*/,
         const int * /*ia*/,
         const int * /*ja*/,
         const int * /*desca*/,
         number * /*S*/,
         number * /*U*/,
         const int * /*iu*/,
         const int * /*ju*/,
         const int * /*descu*/,
         number * /*VT*/,
         const int * /*ivt*/,
         const int * /*jvt*/,
         const int * /*descvt*/,
         number * /*work*/,
         int * /*lwork*/,
         int * /*info*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgesvd(const char *jobu,
         const char *jobvt,
         const int  *m,
         const int  *n,
         double     *A,
         const int  *ia,
         const int  *ja,
         const int  *desca,
         double     *S,
         double     *U,
         const int  *iu,
         const int  *ju,
         const int  *descu,
         double     *VT,
         const int  *ivt,
         const int  *jvt,
         const int  *descvt,
         double     *work,
         int        *lwork,
         int        *info)
  {
    pdgesvd_(jobu,
             jobvt,
             m,
             n,
             A,
             ia,
             ja,
             desca,
             S,
             U,
             iu,
             ju,
             descu,
             VT,
             ivt,
             jvt,
             descvt,
             work,
             lwork,
             info);
  }
 
  inline void
  pgesvd(const char *jobu,
         const char *jobvt,
         const int  *m,
         const int  *n,
         float      *A,
         const int  *ia,
         const int  *ja,
         const int  *desca,
         float      *S,
         float      *U,
         const int  *iu,
         const int  *ju,
         const int  *descu,
         float      *VT,
         const int  *ivt,
         const int  *jvt,
         const int  *descvt,
         float      *work,
         int        *lwork,
         int        *info)
  {
    psgesvd_(jobu,
             jobvt,
             m,
             n,
             A,
             ia,
             ja,
             desca,
             S,
             U,
             iu,
             ju,
             descu,
             VT,
             ivt,
             jvt,
             descvt,
             work,
             lwork,
             info);
  }
 
 
  template <typename number>
  inline void
  pgels(const char * /*trans*/,
        const int * /*m*/,
        const int * /*n*/,
        const int * /*nrhs*/,
        number * /*A*/,
        const int * /*ia*/,
        const int * /*ja*/,
        const int * /*desca*/,
        number * /*B*/,
        const int * /*ib*/,
        const int * /*jb*/,
        const int * /*descb*/,
        number * /*work*/,
        int * /*lwork*/,
        int * /*info*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgels(const char *trans,
        const int  *m,
        const int  *n,
        const int  *nrhs,
        double     *A,
        const int  *ia,
        const int  *ja,
        const int  *desca,
        double     *B,
        const int  *ib,
        const int  *jb,
        const int  *descb,
        double     *work,
        int        *lwork,
        int        *info)
  {
    pdgels_(
      trans, m, n, nrhs, A, ia, ja, desca, B, ib, jb, descb, work, lwork, info);
  }
 
  inline void
  pgels(const char *trans,
        const int  *m,
        const int  *n,
        const int  *nrhs,
        float      *A,
        const int  *ia,
        const int  *ja,
        const int  *desca,
        float      *B,
        const int  *ib,
        const int  *jb,
        const int  *descb,
        float      *work,
        int        *lwork,
        int        *info)
  {
    psgels_(
      trans, m, n, nrhs, A, ia, ja, desca, B, ib, jb, descb, work, lwork, info);
  }
 
 
  template <typename number>
  inline void
  pgeadd(const char * /*transa*/,
         const int * /*m*/,
         const int * /*n*/,
         const number * /*alpha*/,
         const number * /*A*/,
         const int * /*IA*/,
         const int * /*JA*/,
         const int * /*DESCA*/,
         const number * /*beta*/,
         number * /*C*/,
         const int * /*IC*/,
         const int * /*JC*/,
         const int * /*DESCC*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  pgeadd(const char   *transa,
         const int    *m,
         const int    *n,
         const double *alpha,
         const double *A,
         const int    *IA,
         const int    *JA,
         const int    *DESCA,
         const double *beta,
         double       *C,
         const int    *IC,
         const int    *JC,
         const int    *DESCC)
  {
    pdgeadd_(transa, m, n, alpha, A, IA, JA, DESCA, beta, C, IC, JC, DESCC);
  }
 
  inline void
  pgeadd(const char  *transa,
         const int   *m,
         const int   *n,
         const float *alpha,
         const float *A,
         const int   *IA,
         const int   *JA,
         const int   *DESCA,
         const float *beta,
         float       *C,
         const int   *IC,
         const int   *JC,
         const int   *DESCC)
  {
    psgeadd_(transa, m, n, alpha, A, IA, JA, DESCA, beta, C, IC, JC, DESCC);
  }
 
  inline void
  pgeadd(const char                 *transa,
         const int                  *m,
         const int                  *n,
         const std::complex<double> *alpha,
         const std::complex<double> *A,
         const int                  *IA,
         const int                  *JA,
         const int                  *DESCA,
         const std::complex<double> *beta,
         std::complex<double>       *C,
         const int                  *IC,
         const int                  *JC,
         const int                  *DESCC)
  {
    pzgeadd_(transa, m, n, alpha, A, IA, JA, DESCA, beta, C, IC, JC, DESCC);
  }
 
  template <typename number>
  inline void
  ptran(const int * /*m*/,
        const int * /*n*/,
        const number * /*alpha*/,
        const number * /*A*/,
        const int * /*IA*/,
        const int * /*JA*/,
        const int * /*DESCA*/,
        const number * /*beta*/,
        number * /*C*/,
        const int * /*IC*/,
        const int * /*JC*/,
        const int * /*DESCC*/)
  {
    AssertThrow(false, dealii::ExcNotImplemented());
  }
 
  inline void
  ptran(const int    *m,
        const int    *n,
        const double *alpha,
        const double *A,
        const int    *IA,
        const int    *JA,
        const int    *DESCA,
        const double *beta,
        double       *C,
        const int    *IC,
        const int    *JC,
        const int    *DESCC)
  {
    pdtran_(m, n, alpha, A, IA, JA, DESCA, beta, C, IC, JC, DESCC);
  }
 
  inline void
  ptran(const int   *m,
        const int   *n,
        const float *alpha,
        const float *A,
        const int   *IA,
        const int   *JA,
        const int   *DESCA,
        const float *beta,
        float       *C,
        const int   *IC,
        const int   *JC,
        const int   *DESCC)
  {
    pstran_(m, n, alpha, A, IA, JA, DESCA, beta, C, IC, JC, DESCC);
  }
 
 
 
  inline void
  psyevr(const char   *jobz,
         const char   *range,
         const char   *uplo,
         const int    *n,
         double       *A,
         const int    *IA,
         const int    *JA,
         const int    *DESCA,
         const double *VL,
         const double *VU,
         const int    *IL,
         const int    *IU,
         int          *m,
         int          *nz,
         double       *w,
         double       *Z,
         const int    *IZ,
         const int    *JZ,
         const int    *DESCZ,
         double       *work,
         int          *lwork,
         int          *iwork,
         int          *liwork,
         int          *info)
  {
    /*
     * Netlib ScaLAPACK performs floating point tests (e.g. divide-by-zero)
     * within the call to pdsyevr causing floating point exceptions to be thrown
     * (at least in debug mode). Therefore, we wrap the calls to pdsyevr into
     * the following code to suppress the exception.
     */
#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
    fenv_t fp_exceptions;
    feholdexcept(&fp_exceptions);
#endif
 
    pdsyevr_(jobz,
             range,
             uplo,
             n,
             A,
             IA,
             JA,
             DESCA,
             VL,
             VU,
             IL,
             IU,
             m,
             nz,
             w,
             Z,
             IZ,
             JZ,
             DESCZ,
             work,
             lwork,
             iwork,
             liwork,
             info);
 
#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
    fesetenv(&fp_exceptions);
#endif
  }
 
  inline void
  psyevr(const char  *jobz,
         const char  *range,
         const char  *uplo,
         const int   *n,
         float       *A,
         const int   *IA,
         const int   *JA,
         const int   *DESCA,
         const float *VL,
         const float *VU,
         const int   *IL,
         const int   *IU,
         int         *m,
         int         *nz,
         float       *w,
         float       *Z,
         const int   *IZ,
         const int   *JZ,
         const int   *DESCZ,
         float       *work,
         int         *lwork,
         int         *iwork,
         int         *liwork,
         int         *info)
  {
    /*
     * Netlib ScaLAPACK performs floating point tests (e.g. divide-by-zero)
     * within the call to pssyevr causing floating point exceptions to be thrown
     * (at least in debug mode). Therefore, we wrap the calls to pssyevr into
     * the following code to suppress the exception.
     */
#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
    fenv_t fp_exceptions;
    feholdexcept(&fp_exceptions);
#endif
 
    pssyevr_(jobz,
             range,
             uplo,
             n,
             A,
             IA,
             JA,
             DESCA,
             VL,
             VU,
             IL,
             IU,
             m,
             nz,
             w,
             Z,
             IZ,
             JZ,
             DESCZ,
             work,
             lwork,
             iwork,
             liwork,
             info);
 
#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
    fesetenv(&fp_exceptions);
#endif
  }
 
  inline void
  psyevr(const char           *jobz,
         const char           *range,
         const char           *uplo,
         const int            *n,
         std::complex<double> *A,
         const int            *IA,
         const int            *JA,
         const int            *DESCA,
         const double         *VL,
         const double         *VU,
         const int            *IL,
         const int            *IU,
         int                  *m,
         int                  *nz,
         double               *w,
         std::complex<double> *Z,
         const int            *IZ,
         const int            *JZ,
         const int            *DESCZ,
         std::complex<double> *work,
         int                  *lwork,
         double               *rwork,
         int                  *lrwork,
         int                  *iwork,
         int                  *liwork,
         int                  *info)
  {
    /*
     * Netlib ScaLAPACK performs floating point tests (e.g. divide-by-zero)
     * within the call to pdsyevr causing floating point exceptions to be thrown
     * (at least in debug mode). Therefore, we wrap the calls to pdsyevr into
     * the following code to suppress the exception.
     */
#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
    fenv_t fp_exceptions;
    feholdexcept(&fp_exceptions);
#endif
 
    pzheevr_(jobz,
             range,
             uplo,
             n,
             A,
             IA,
             JA,
             DESCA,
             VL,
             VU,
             IL,
             IU,
             m,
             nz,
             w,
             Z,
             IZ,
             JZ,
             DESCZ,
             work,
             lwork,
             rwork,
             lrwork,
             iwork,
             liwork,
             info);
 
#ifdef DEAL_II_HAVE_FP_EXCEPTIONS
    fesetenv(&fp_exceptions);
#endif
  }
 
  inline int
  lworkFromWork(std::vector<double> &work)
  {
    return static_cast<int>(work[0]);
  }
 
  inline int
  lworkFromWork(std::vector<float> &work)
  {
    return static_cast<int>(work[0]);
  }
 
  inline int
  lworkFromWork(std::vector<std::complex<double>> &work)
  {
    return static_cast<int>(work[0].real());
  }
} // namespace dftfe
#endif // dftfe_scalapack_templates_h