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             :    BV private kernels that use the BLAS
      12             : */
      13             : 
      14             : #include <slepc/private/bvimpl.h>
      15             : #include <slepcblaslapack.h>
      16             : 
      17             : #define BLOCKSIZE 64
      18             : 
      19             : /*
      20             :     C := alpha*A*B + beta*C
      21             : 
      22             :     A is mxk (ld=lda), B is kxn (ld=ldb), C is mxn (ld=ldc)
      23             : */
      24        5104 : PetscErrorCode BVMult_BLAS_Private(BV bv,PetscInt m_,PetscInt n_,PetscInt k_,PetscScalar alpha,const PetscScalar *A,PetscInt lda_,const PetscScalar *B,PetscInt ldb_,PetscScalar beta,PetscScalar *C,PetscInt ldc_)
      25             : {
      26        5104 :   PetscBLASInt   m,n,k,lda,ldb,ldc;
      27             : #if defined(PETSC_HAVE_FBLASLAPACK) || defined(PETSC_HAVE_F2CBLASLAPACK)
      28             :   PetscBLASInt   l,bs=BLOCKSIZE;
      29             : #endif
      30             : 
      31        5104 :   PetscFunctionBegin;
      32        5104 :   PetscCall(PetscBLASIntCast(m_,&m));
      33        5104 :   PetscCall(PetscBLASIntCast(n_,&n));
      34        5104 :   PetscCall(PetscBLASIntCast(k_,&k));
      35        5104 :   PetscCall(PetscBLASIntCast(lda_,&lda));
      36        5104 :   PetscCall(PetscBLASIntCast(ldb_,&ldb));
      37        5104 :   PetscCall(PetscBLASIntCast(ldc_,&ldc));
      38             : #if defined(PETSC_HAVE_FBLASLAPACK) || defined(PETSC_HAVE_F2CBLASLAPACK)
      39             :   l = m % bs;
      40             :   if (l) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&k,&alpha,(PetscScalar*)A,&lda,(PetscScalar*)B,&ldb,&beta,C,&ldc));
      41             :   for (;l<m;l+=bs) {
      42             :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&bs,&n,&k,&alpha,(PetscScalar*)A+l,&lda,(PetscScalar*)B,&ldb,&beta,C+l,&ldc));
      43             :   }
      44             : #else
      45        5104 :   if (m) PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&alpha,(PetscScalar*)A,&lda,(PetscScalar*)B,&ldb,&beta,C,&ldc));
      46             : #endif
      47        5104 :   PetscCall(PetscLogFlops(2.0*m*n*k));
      48        5104 :   PetscFunctionReturn(PETSC_SUCCESS);
      49             : }
      50             : 
      51             : /*
      52             :     y := alpha*A*x + beta*y
      53             : 
      54             :     A is nxk (ld=lda)
      55             : */
      56      719898 : PetscErrorCode BVMultVec_BLAS_Private(BV bv,PetscInt n_,PetscInt k_,PetscScalar alpha,const PetscScalar *A,PetscInt lda_,const PetscScalar *x,PetscScalar beta,PetscScalar *y)
      57             : {
      58      719898 :   PetscBLASInt   n,k,lda,one=1;
      59             : 
      60      719898 :   PetscFunctionBegin;
      61      719898 :   PetscCall(PetscBLASIntCast(n_,&n));
      62      719898 :   PetscCall(PetscBLASIntCast(k_,&k));
      63      719898 :   PetscCall(PetscBLASIntCast(lda_,&lda));
      64      719898 :   if (n) PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&k,&alpha,A,&lda,x,&one,&beta,y,&one));
      65      719898 :   PetscCall(PetscLogFlops(2.0*n*k));
      66      719898 :   PetscFunctionReturn(PETSC_SUCCESS);
      67             : }
      68             : 
      69             : /*
      70             :     A(:,s:e-1) := A*B(:,s:e-1)
      71             : 
      72             :     A is mxk (ld=lda), B is kxn (ld=ldb), n=e-s
      73             : */
      74       23743 : PetscErrorCode BVMultInPlace_BLAS_Private(BV bv,PetscInt m_,PetscInt k_,PetscInt s,PetscInt e,PetscScalar *A,PetscInt lda_,const PetscScalar *B,PetscInt ldb_,PetscBool btrans)
      75             : {
      76       23743 :   PetscScalar    *pb,zero=0.0,one=1.0;
      77       23743 :   PetscBLASInt   m,n,k,l,lda,ldb,bs=BLOCKSIZE;
      78       23743 :   PetscInt       j,n_=e-s;
      79       23743 :   const char     *bt;
      80             : 
      81       23743 :   PetscFunctionBegin;
      82       23743 :   PetscCall(PetscBLASIntCast(m_,&m));
      83       23743 :   PetscCall(PetscBLASIntCast(n_,&n));
      84       23743 :   PetscCall(PetscBLASIntCast(k_,&k));
      85       23743 :   PetscCall(PetscBLASIntCast(lda_,&lda));
      86       23743 :   PetscCall(PetscBLASIntCast(ldb_,&ldb));
      87       23743 :   PetscCall(BVAllocateWork_Private(bv,BLOCKSIZE*n_));
      88       23743 :   if (PetscUnlikely(btrans)) {
      89           3 :     pb = (PetscScalar*)B+s;
      90           3 :     bt = "C";
      91             :   } else {
      92       23740 :     pb = (PetscScalar*)B+s*ldb;
      93       23740 :     bt = "N";
      94             :   }
      95       23743 :   l = m % bs;
      96       23743 :   if (l) {
      97       22497 :     PetscCallBLAS("BLASgemm",BLASgemm_("N",bt,&l,&n,&k,&one,A,&lda,pb,&ldb,&zero,bv->work,&l));
      98      152513 :     for (j=0;j<n;j++) PetscCall(PetscArraycpy(A+(s+j)*lda,bv->work+j*l,l));
      99             :   }
     100      104823 :   for (;l<m;l+=bs) {
     101       81080 :     PetscCallBLAS("BLASgemm",BLASgemm_("N",bt,&bs,&n,&k,&one,A+l,&lda,pb,&ldb,&zero,bv->work,&bs));
     102      738264 :     for (j=0;j<n;j++) PetscCall(PetscArraycpy(A+(s+j)*lda+l,bv->work+j*bs,bs));
     103             :   }
     104       23743 :   PetscCall(PetscLogFlops(2.0*m*n*k));
     105       23743 :   PetscFunctionReturn(PETSC_SUCCESS);
     106             : }
     107             : 
     108             : /*
     109             :     V := V*B
     110             : 
     111             :     V is mxn (ld=m), B is nxn (ld=k)
     112             : */
     113         115 : PetscErrorCode BVMultInPlace_Vecs_Private(BV bv,PetscInt m_,PetscInt n_,PetscInt k_,Vec *V,const PetscScalar *B,PetscBool btrans)
     114             : {
     115         115 :   PetscScalar       zero=0.0,one=1.0,*out,*pout;
     116         115 :   const PetscScalar *pin;
     117         115 :   PetscBLASInt      m = 0,n,k,l,bs=BLOCKSIZE;
     118         115 :   PetscInt          j;
     119         115 :   const char        *bt;
     120             : 
     121         115 :   PetscFunctionBegin;
     122         115 :   PetscCall(PetscBLASIntCast(m_,&m));
     123         115 :   PetscCall(PetscBLASIntCast(n_,&n));
     124         115 :   PetscCall(PetscBLASIntCast(k_,&k));
     125         115 :   PetscCall(BVAllocateWork_Private(bv,2*BLOCKSIZE*n_));
     126         115 :   out = bv->work+BLOCKSIZE*n_;
     127         115 :   if (btrans) bt = "C";
     128         114 :   else bt = "N";
     129         115 :   l = m % bs;
     130         115 :   if (l) {
     131         563 :     for (j=0;j<n;j++) {
     132         448 :       PetscCall(VecGetArrayRead(V[j],&pin));
     133         448 :       PetscCall(PetscArraycpy(bv->work+j*l,pin,l));
     134         448 :       PetscCall(VecRestoreArrayRead(V[j],&pin));
     135             :     }
     136         115 :     PetscCallBLAS("BLASgemm",BLASgemm_("N",bt,&l,&n,&n,&one,bv->work,&l,(PetscScalar*)B,&k,&zero,out,&l));
     137         563 :     for (j=0;j<n;j++) {
     138         448 :       PetscCall(VecGetArray(V[j],&pout));
     139         448 :       PetscCall(PetscArraycpy(pout,out+j*l,l));
     140         448 :       PetscCall(VecRestoreArray(V[j],&pout));
     141             :     }
     142             :   }
     143         523 :   for (;l<m;l+=bs) {
     144        1956 :     for (j=0;j<n;j++) {
     145        1548 :       PetscCall(VecGetArrayRead(V[j],&pin));
     146        1548 :       PetscCall(PetscArraycpy(bv->work+j*bs,pin+l,bs));
     147        1548 :       PetscCall(VecRestoreArrayRead(V[j],&pin));
     148             :     }
     149         408 :     PetscCallBLAS("BLASgemm",BLASgemm_("N",bt,&bs,&n,&n,&one,bv->work,&bs,(PetscScalar*)B,&k,&zero,out,&bs));
     150        1956 :     for (j=0;j<n;j++) {
     151        1548 :       PetscCall(VecGetArray(V[j],&pout));
     152        1548 :       PetscCall(PetscArraycpy(pout+l,out+j*bs,bs));
     153        1548 :       PetscCall(VecRestoreArray(V[j],&pout));
     154             :     }
     155             :   }
     156         115 :   PetscCall(PetscLogFlops(2.0*n*n*k));
     157         115 :   PetscFunctionReturn(PETSC_SUCCESS);
     158             : }
     159             : 
     160             : /*
     161             :     B := alpha*A + beta*B
     162             : 
     163             :     A,B are nxk
     164             : */
     165        1492 : PetscErrorCode BVAXPY_BLAS_Private(BV bv,PetscInt n_,PetscInt k_,PetscScalar alpha,const PetscScalar *A,PetscInt lda_,PetscScalar beta,PetscScalar *B,PetscInt ldb_)
     166             : {
     167        1492 :   PetscBLASInt   m,one=1;
     168        1492 :   PetscInt       i,j;
     169             : 
     170        1492 :   PetscFunctionBegin;
     171        1492 :   if (lda_==n_ && ldb_==n_) {
     172        1492 :     PetscCall(PetscBLASIntCast(n_*k_,&m));
     173        1492 :     if (beta!=(PetscScalar)1.0) PetscCallBLAS("BLASscal",BLASscal_(&m,&beta,B,&one));
     174        1492 :     PetscCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,A,&one,B,&one));
     175             :   } else {
     176           0 :     if (beta!=(PetscScalar)1.0) {
     177           0 :       for (j=0;j<k_;j++) {
     178           0 :         for (i=0;i<n_;i++) {
     179           0 :           B[i+j*ldb_] = alpha*A[i+j*lda_] + beta*B[i+j*ldb_];
     180             :         }
     181             :       }
     182             :     } else {
     183           0 :       for (j=0;j<k_;j++) {
     184           0 :         for (i=0;i<n_;i++) {
     185           0 :           B[i+j*ldb_] += alpha*A[i+j*lda_];
     186             :         }
     187             :       }
     188             :     }
     189             :   }
     190        1492 :   PetscCall(PetscLogFlops((beta==(PetscScalar)1.0)?2.0*n_*k_:3.0*n_*k_));
     191        1492 :   PetscFunctionReturn(PETSC_SUCCESS);
     192             : }
     193             : 
     194             : /*
     195             :     C := A'*B
     196             : 
     197             :     A' is mxk (ld=lda), B is kxn (ld=ldb), C is mxn (ld=ldc)
     198             : */
     199       30088 : PetscErrorCode BVDot_BLAS_Private(BV bv,PetscInt m_,PetscInt n_,PetscInt k_,const PetscScalar *A,PetscInt lda_,const PetscScalar *B,PetscInt ldb_,PetscScalar *C,PetscInt ldc_,PetscBool mpi)
     200             : {
     201       30088 :   PetscScalar    zero=0.0,one=1.0,*CC;
     202       30088 :   PetscBLASInt   m,n,k,lda,ldb,ldc,j;
     203       30088 :   PetscMPIInt    len;
     204             : 
     205       30088 :   PetscFunctionBegin;
     206       30088 :   PetscCall(PetscBLASIntCast(m_,&m));
     207       30088 :   PetscCall(PetscBLASIntCast(n_,&n));
     208       30088 :   PetscCall(PetscBLASIntCast(k_,&k));
     209       30088 :   PetscCall(PetscBLASIntCast(lda_,&lda));
     210       30088 :   PetscCall(PetscBLASIntCast(ldb_,&ldb));
     211       30088 :   PetscCall(PetscBLASIntCast(ldc_,&ldc));
     212       30088 :   if (mpi) {
     213        3941 :     if (ldc==m) {
     214        2562 :       PetscCall(BVAllocateWork_Private(bv,m*n));
     215        2562 :       if (k) PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&m,&n,&k,&one,(PetscScalar*)A,&lda,(PetscScalar*)B,&ldb,&zero,bv->work,&ldc));
     216           0 :       else PetscCall(PetscArrayzero(bv->work,m*n));
     217        2562 :       PetscCall(PetscMPIIntCast(m*n,&len));
     218        7686 :       PetscCallMPI(MPIU_Allreduce(bv->work,C,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)bv)));
     219             :     } else {
     220        1379 :       PetscCall(BVAllocateWork_Private(bv,2*m*n));
     221        1379 :       CC = bv->work+m*n;
     222        1379 :       if (k) PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&m,&n,&k,&one,(PetscScalar*)A,&lda,(PetscScalar*)B,&ldb,&zero,bv->work,&m));
     223           0 :       else PetscCall(PetscArrayzero(bv->work,m*n));
     224        1379 :       PetscCall(PetscMPIIntCast(m*n,&len));
     225        4137 :       PetscCallMPI(MPIU_Allreduce(bv->work,CC,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)bv)));
     226       16641 :       for (j=0;j<n;j++) PetscCall(PetscArraycpy(C+j*ldc,CC+j*m,m));
     227             :     }
     228             :   } else {
     229       26147 :     if (k) PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&m,&n,&k,&one,(PetscScalar*)A,&lda,(PetscScalar*)B,&ldb,&zero,C,&ldc));
     230             :   }
     231       30088 :   PetscCall(PetscLogFlops(2.0*m*n*k));
     232       30088 :   PetscFunctionReturn(PETSC_SUCCESS);
     233             : }
     234             : 
     235             : /*
     236             :     y := A'*x
     237             : 
     238             :     A is nxk (ld=lda)
     239             : */
     240      420118 : PetscErrorCode BVDotVec_BLAS_Private(BV bv,PetscInt n_,PetscInt k_,const PetscScalar *A,PetscInt lda_,const PetscScalar *x,PetscScalar *y,PetscBool mpi)
     241             : {
     242      420118 :   PetscScalar    zero=0.0,done=1.0;
     243      420118 :   PetscBLASInt   n,k,lda,one=1;
     244      420118 :   PetscMPIInt    len;
     245             : 
     246      420118 :   PetscFunctionBegin;
     247      420118 :   PetscCall(PetscBLASIntCast(n_,&n));
     248      420118 :   PetscCall(PetscBLASIntCast(k_,&k));
     249      420118 :   PetscCall(PetscBLASIntCast(lda_,&lda));
     250      420118 :   if (mpi) {
     251       36244 :     PetscCall(BVAllocateWork_Private(bv,k));
     252       36244 :     if (n) PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&k,&done,A,&lda,x,&one,&zero,bv->work,&one));
     253           0 :     else PetscCall(PetscArrayzero(bv->work,k));
     254       36244 :     PetscCall(PetscMPIIntCast(k,&len));
     255      108732 :     PetscCallMPI(MPIU_Allreduce(bv->work,y,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)bv)));
     256             :   } else {
     257      383874 :     if (n) PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&k,&done,A,&lda,x,&one,&zero,y,&one));
     258             :   }
     259      420118 :   PetscCall(PetscLogFlops(2.0*n*k));
     260      420118 :   PetscFunctionReturn(PETSC_SUCCESS);
     261             : }
     262             : 
     263             : /*
     264             :     Scale n scalars
     265             : */
     266      194914 : PetscErrorCode BVScale_BLAS_Private(BV bv,PetscInt n_,PetscScalar *A,PetscScalar alpha)
     267             : {
     268      194914 :   PetscBLASInt   n,one=1;
     269             : 
     270      194914 :   PetscFunctionBegin;
     271      194914 :   if (PetscUnlikely(alpha == (PetscScalar)0.0)) PetscCall(PetscArrayzero(A,n_));
     272      194558 :   else if (alpha!=(PetscScalar)1.0) {
     273      194555 :     PetscCall(PetscBLASIntCast(n_,&n));
     274      194555 :     PetscCallBLAS("BLASscal",BLASscal_(&n,&alpha,A,&one));
     275      194555 :     PetscCall(PetscLogFlops(1.0*n));
     276             :   }
     277      194914 :   PetscFunctionReturn(PETSC_SUCCESS);
     278             : }