Actual source code: bvec1.c

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

     7:  #include vecimpl.h
     8:  #include src/vec/impls/dvecimpl.h
     9:  #include petscblaslapack.h

     13: PetscErrorCode VecDot_Seq(Vec xin,Vec yin,PetscScalar *z)
     14: {
     15:   PetscScalar *ya,*xa;
     17: #if !defined(PETSC_USE_COMPLEX)
     18:   PetscBLASInt bn = (PetscBLASInt)xin->n, one = 1;
     19: #endif

     22:   VecGetArray(xin,&xa);
     23:   if (xin != yin) {VecGetArray(yin,&ya);}
     24:   else ya = xa;
     25: #if defined(PETSC_USE_COMPLEX)
     26:   /* cannot use BLAS dot for complex because compiler/linker is
     27:      not happy about returning a double complex */
     28:   {
     29:     PetscInt    i;
     30:     PetscScalar sum = 0.0;
     31:     for (i=0; i<xin->n; i++) {
     32:       sum += xa[i]*PetscConj(ya[i]);
     33:     }
     34:     *z = sum;
     35:   }
     36: #else
     37:   *z = BLASdot_(&bn,xa,&one,ya,&one);
     38: #endif
     39:   VecRestoreArray(xin,&xa);
     40:   if (xin != yin) {VecRestoreArray(yin,&ya);}
     41:   PetscLogFlops(2*xin->n-1);
     42:   return(0);
     43: }

     47: PetscErrorCode VecTDot_Seq(Vec xin,Vec yin,PetscScalar *z)
     48: {
     49:   PetscScalar *ya,*xa;
     51: #if !defined(PETSC_USE_COMPLEX)
     52:  PetscBLASInt bn = (PetscBLASInt)xin->n, one = 1;
     53: #endif

     56:   VecGetArray(xin,&xa);
     57:   if (xin != yin) {VecGetArray(yin,&ya);}
     58:   else ya = xa;
     59: #if defined(PETSC_USE_COMPLEX)
     60:   /* cannot use BLAS dot for complex because compiler/linker is
     61:      not happy about returning a double complex */
     62:   PetscInt    i;
     63:   PetscScalar sum = 0.0;
     64:   for (i=0; i<xin->n; i++) {
     65:     sum += xa[i]*ya[i];
     66:   }
     67:   *z = sum;
     68: #else
     69:   *z = BLASdot_(&bn,xa,&one,ya,&one);
     70: #endif
     71:   VecRestoreArray(xin,&xa);
     72:   if (xin != yin) {VecRestoreArray(yin,&ya);}
     73:   PetscLogFlops(2*xin->n-1);
     74:   return(0);
     75: }

     79: PetscErrorCode VecScale_Seq(Vec xin, PetscScalar alpha)
     80: {
     81:   Vec_Seq        *x = (Vec_Seq*)xin->data;
     82:   PetscBLASInt   bn = (PetscBLASInt)xin->n, one = 1;

     86:   if (alpha == 0.0) {
     87:     VecSet_Seq(xin,alpha);
     88:   } else if (alpha != 1.0) {
     89:     PetscScalar a = alpha;
     90:     BLASscal_(&bn,&a,x->array,&one);
     91:     PetscLogFlops(xin->n);
     92:   }
     93:   return(0);
     94: }

     98: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
     99: {
    100:   Vec_Seq     *x = (Vec_Seq *)xin->data;
    101:   PetscScalar *ya;

    105:   if (xin != yin) {
    106:     VecGetArray(yin,&ya);
    107:     PetscMemcpy(ya,x->array,xin->n*sizeof(PetscScalar));
    108:     VecRestoreArray(yin,&ya);
    109:   }
    110:   return(0);
    111: }

    115: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
    116: {
    117:   Vec_Seq        *x = (Vec_Seq *)xin->data;
    118:   PetscScalar    *ya;
    120:   PetscBLASInt   bn = (PetscBLASInt)xin->n, one = 1;

    123:   if (xin != yin) {
    124:     VecGetArray(yin,&ya);
    125:     BLASswap_(&bn,x->array,&one,ya,&one);
    126:     VecRestoreArray(yin,&ya);
    127:   }
    128:   return(0);
    129: }

    133: PetscErrorCode VecAXPY_Seq(Vec yin,PetscScalar alpha,Vec xin)
    134: {
    135:   Vec_Seq        *y = (Vec_Seq *)yin->data;
    137:   PetscBLASInt   bn = (PetscBLASInt)yin->n, one = 1;
    138:   PetscScalar    *xarray;

    141:   /* assume that the BLAS handles alpha == 1.0 efficiently since we have no fast code for it */
    142:   if (alpha != 0.0) {
    143:     PetscScalar oalpha = alpha;
    144:     VecGetArray(xin,&xarray);
    145:     BLASaxpy_(&bn,&oalpha,xarray,&one,y->array,&one);
    146:     VecRestoreArray(xin,&xarray);
    147:     PetscLogFlops(2*yin->n);
    148:   }
    149:   return(0);
    150: }

    154: PetscErrorCode VecAXPBY_Seq(Vec yin,PetscScalar alpha,PetscScalar beta,Vec xin)
    155: {
    156:   Vec_Seq        *y = (Vec_Seq *)yin->data;
    158:   PetscInt       n = yin->n,i;
    159:   PetscScalar    *yy = y->array,*xx ,a = alpha,b = beta;

    162:   if (a == 0.0) {
    163:     VecScale_Seq(yin,beta);
    164:   } else if (b == 1.0) {
    165:     VecAXPY_Seq(yin,alpha,xin);
    166:   } else if (a == 1.0) {
    167:     VecAYPX_Seq(yin,beta,xin);
    168:   } else if (b == 0.0) {
    169:     VecGetArray(xin,&xx);
    170:     for (i=0; i<n; i++) {
    171:       yy[i] = a*xx[i];
    172:     }
    173:     VecRestoreArray(xin,&xx);
    174:     PetscLogFlops(xin->n);
    175:   } else {
    176:     VecGetArray(xin,&xx);
    177:     for (i=0; i<n; i++) {
    178:       yy[i] = a*xx[i] + b*yy[i];
    179:     }
    180:     VecRestoreArray(xin,&xx);
    181:     PetscLogFlops(3*xin->n);
    182:   }
    183:   return(0);
    184: }