1: /*
  2:    Defines the BLAS based vector operations. Code shared by parallel
  3:   and sequential vectors.
  4: */

  6: #include <../src/vec/vec/impls/dvecimpl.h>
  7: #include <petscblaslapack.h>

  9: #if defined(PETSC_USE_REAL_SINGLE) && defined(PETSC_BLASLAPACK_SNRM2_RETURNS_DOUBLE) && !defined(PETSC_USE_COMPLEX)
 10: static PetscErrorCode VecXDot_Seq_Private(Vec xin, Vec yin, PetscScalar *z, double (*const BLASfn)(const PetscBLASInt *, const PetscScalar *, const PetscBLASInt *, const PetscScalar *, const PetscBLASInt *))
 11: #else
 12: static PetscErrorCode VecXDot_Seq_Private(Vec xin, Vec yin, PetscScalar *z, PetscScalar (*const BLASfn)(const PetscBLASInt *, const PetscScalar *, const PetscBLASInt *, const PetscScalar *, const PetscBLASInt *))
 13: #endif
 14: {
 15:   const PetscInt     n   = xin->map->n;
 16:   const PetscBLASInt one = 1;
 17:   const PetscScalar *ya, *xa;
 18:   PetscBLASInt       bn;

 20:   PetscFunctionBegin;
 21:   PetscCall(PetscBLASIntCast(n, &bn));
 22:   if (n > 0) PetscCall(PetscLogFlops(2.0 * n - 1));
 23:   PetscCall(VecGetArrayRead(xin, &xa));
 24:   PetscCall(VecGetArrayRead(yin, &ya));
 25:   /* arguments ya, xa are reversed because BLAS complex conjugates the first argument, PETSc
 26:      the second */
 27:   PetscCallBLAS("BLASdot", *z = BLASfn(&bn, ya, &one, xa, &one));
 28:   PetscCall(VecRestoreArrayRead(xin, &xa));
 29:   PetscCall(VecRestoreArrayRead(yin, &ya));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: PetscErrorCode VecDot_Seq(Vec xin, Vec yin, PetscScalar *z)
 34: {
 35:   PetscFunctionBegin;
 36:   PetscCall(VecXDot_Seq_Private(xin, yin, z, BLASdot_));
 37:   PetscFunctionReturn(PETSC_SUCCESS);
 38: }

 40: PetscErrorCode VecTDot_Seq(Vec xin, Vec yin, PetscScalar *z)
 41: {
 42:   PetscFunctionBegin;
 43:   /*
 44:     pay close attention!!! xin and yin are SWAPPED here so that the eventual BLAS call is
 45:     dot(&bn, xa, &one, ya, &one)
 46:   */
 47:   PetscCall(VecXDot_Seq_Private(yin, xin, z, BLASdotu_));
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

 51: PetscErrorCode VecScale_Seq(Vec xin, PetscScalar alpha)
 52: {
 53:   PetscFunctionBegin;
 54:   if (alpha == (PetscScalar)0.0) {
 55:     PetscCall(VecSet_Seq(xin, alpha));
 56:   } else if (alpha != (PetscScalar)1.0) {
 57:     const PetscBLASInt one = 1;
 58:     PetscBLASInt       bn;
 59:     PetscScalar       *xarray;

 61:     PetscCall(PetscBLASIntCast(xin->map->n, &bn));
 62:     PetscCall(PetscLogFlops(bn));
 63:     PetscCall(VecGetArray(xin, &xarray));
 64:     PetscCallBLAS("BLASscal", BLASscal_(&bn, &alpha, xarray, &one));
 65:     PetscCall(VecRestoreArray(xin, &xarray));
 66:   }
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

 70: PetscErrorCode VecAXPY_Seq(Vec yin, PetscScalar alpha, Vec xin)
 71: {
 72:   PetscFunctionBegin;
 73:   /* assume that the BLAS handles alpha == 1.0 efficiently since we have no fast code for it */
 74:   if (alpha != (PetscScalar)0.0) {
 75:     const PetscScalar *xarray;
 76:     PetscScalar       *yarray;
 77:     const PetscBLASInt one = 1;
 78:     PetscBLASInt       bn;

 80:     PetscCall(PetscBLASIntCast(yin->map->n, &bn));
 81:     PetscCall(PetscLogFlops(2.0 * bn));
 82:     PetscCall(VecGetArrayRead(xin, &xarray));
 83:     PetscCall(VecGetArray(yin, &yarray));
 84:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bn, &alpha, xarray, &one, yarray, &one));
 85:     PetscCall(VecRestoreArrayRead(xin, &xarray));
 86:     PetscCall(VecRestoreArray(yin, &yarray));
 87:   }
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: PetscErrorCode VecAXPBY_Seq(Vec yin, PetscScalar a, PetscScalar b, Vec xin)
 92: {
 93:   PetscFunctionBegin;
 94:   if (a == (PetscScalar)0.0) {
 95:     PetscCall(VecScale_Seq(yin, b));
 96:   } else if (b == (PetscScalar)1.0) {
 97:     PetscCall(VecAXPY_Seq(yin, a, xin));
 98:   } else if (a == (PetscScalar)1.0) {
 99:     PetscCall(VecAYPX_Seq(yin, b, xin));
100:   } else {
101:     const PetscInt     n = yin->map->n;
102:     const PetscScalar *xx;
103:     PetscScalar       *yy;

105:     PetscCall(VecGetArrayRead(xin, &xx));
106:     PetscCall(VecGetArray(yin, &yy));
107:     if (b == (PetscScalar)0.0) {
108:       for (PetscInt i = 0; i < n; ++i) yy[i] = a * xx[i];
109:       PetscCall(PetscLogFlops(n));
110:     } else {
111:       for (PetscInt i = 0; i < n; ++i) yy[i] = a * xx[i] + b * yy[i];
112:       PetscCall(PetscLogFlops(3.0 * n));
113:     }
114:     PetscCall(VecRestoreArrayRead(xin, &xx));
115:     PetscCall(VecRestoreArray(yin, &yy));
116:   }
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: PetscErrorCode VecAXPBYPCZ_Seq(Vec zin, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec xin, Vec yin)
121: {
122:   const PetscInt     n = zin->map->n;
123:   const PetscScalar *yy, *xx;
124:   PetscInt           flops = 4 * n; // common case
125:   PetscScalar       *zz;

127:   PetscFunctionBegin;
128:   PetscCall(VecGetArrayRead(xin, &xx));
129:   PetscCall(VecGetArrayRead(yin, &yy));
130:   PetscCall(VecGetArray(zin, &zz));
131:   if (alpha == (PetscScalar)1.0) {
132:     for (PetscInt i = 0; i < n; ++i) zz[i] = xx[i] + beta * yy[i] + gamma * zz[i];
133:   } else if (gamma == (PetscScalar)1.0) {
134:     for (PetscInt i = 0; i < n; ++i) zz[i] = alpha * xx[i] + beta * yy[i] + zz[i];
135:   } else if (gamma == (PetscScalar)0.0) {
136:     for (PetscInt i = 0; i < n; ++i) zz[i] = alpha * xx[i] + beta * yy[i];
137:     flops -= n;
138:   } else {
139:     for (PetscInt i = 0; i < n; ++i) zz[i] = alpha * xx[i] + beta * yy[i] + gamma * zz[i];
140:     flops += n;
141:   }
142:   PetscCall(VecRestoreArrayRead(xin, &xx));
143:   PetscCall(VecRestoreArrayRead(yin, &yy));
144:   PetscCall(VecRestoreArray(zin, &zz));
145:   PetscCall(PetscLogFlops(flops));
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }