BlaVector.c

Go to the documentation of this file.

 
#include <math.h>
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include "fasp.h"
#include "fasp_functs.h"
 
/*---------------------------------*/
/*--      Public Functions       --*/
/*---------------------------------*/
 
void fasp_blas_dvec_axpy (const REAL     a,
                          const dvector *x,
                          dvector       *y)
{
    const INT   m   = x->row;
    const REAL *xpt = x->val;
    REAL       *ypt = y->val;
 
    SHORT use_openmp = FALSE;
    INT   i;
    
#ifdef _OPENMP
    INT myid, mybegin, myend, nthreads;
    if ( m > OPENMP_HOLDS ) {
        use_openmp = TRUE;
        nthreads = fasp_get_num_threads();
    }
#endif
    
    if ( y->row != m ) {
        printf("### ERROR: Vectors have different dimensions!\n");
        fasp_chkerr(ERROR_DATA_STRUCTURE, __FUNCTION__);
    }
    
    if (use_openmp) {
#ifdef _OPENMP
#pragma omp parallel private(myid,mybegin,myend,i) num_threads(nthreads)
        {
            myid = omp_get_thread_num();
            fasp_get_start_end (myid, nthreads, m, &mybegin, &myend);
            for ( i = mybegin; i < myend; ++i ) ypt[i] += a*xpt[i];
        }
#endif
    }
    else {
        for ( i = 0; i < m; ++i ) ypt[i] += a*xpt[i];
    }
}
 
void fasp_blas_dvec_axpyz (const REAL     a,
                           const dvector *x,
                           const dvector *y,
                           dvector       *z)
{
    const INT     m = x->row;
    const REAL *xpt = x->val, *ypt = y->val;
    REAL       *zpt = z->val;
    
    if ( y->row != m ) {
        printf("### ERROR: Vectors have different dimensions!\n");
        fasp_chkerr(ERROR_DATA_STRUCTURE, __FUNCTION__);
    }
    
    z->row = m;
    
    memcpy (zpt, ypt, m*sizeof(REAL));
    fasp_blas_darray_axpy (m, a, xpt, zpt);
}
 
REAL fasp_blas_dvec_norm1 (const dvector *x)
{
    const INT   length = x->row;
    const REAL *xpt = x->val;
    
    register REAL onenorm = 0.0;
    SHORT use_openmp = FALSE;
    INT   i;
    
#ifdef _OPENMP
    if ( length > OPENMP_HOLDS ) {
        use_openmp = TRUE;
    }
#endif
    
    if ( use_openmp ) {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:onenorm) private(i)
#endif
        for ( i = 0; i < length; ++i ) onenorm += ABS(xpt[i]);
    }
    else {
        for ( i = 0; i < length; ++i ) onenorm += ABS(xpt[i]);
    }
    
    return onenorm;
}
 
REAL fasp_blas_dvec_norm2 (const dvector *x)
{
    const INT   length = x->row;
    const REAL *xpt = x->val;
    
    register REAL twonorm = 0.0;
    SHORT use_openmp = FALSE;
    INT   i;
    
#ifdef _OPENMP
    if ( length > OPENMP_HOLDS ) use_openmp = TRUE;
#endif
    
    if ( use_openmp ) {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:twonorm) private(i)
#endif
        for ( i = 0; i < length; ++i ) twonorm += xpt[i]*xpt[i];
    }
    else {
        for ( i = 0; i < length; ++i ) twonorm += xpt[i]*xpt[i];
    }
    
    return sqrt(twonorm);
}
 
REAL fasp_blas_dvec_norminf (const dvector *x)
{
    const INT      length=x->row;
    const REAL    *xpt=x->val;
    
    register REAL  infnorm = 0.0;
    register INT   i;
    
    for ( i = 0; i < length; ++i ) infnorm = MAX(infnorm, ABS(xpt[i]));
    
    return infnorm;
}
 
REAL fasp_blas_dvec_dotprod (const dvector *x,
                             const dvector *y)
{
    const INT   length = x->row;
    const REAL *xpt = x->val, *ypt = y->val;
    
    register REAL value = 0.0;
    SHORT use_openmp = FALSE;
    INT   i;
    
#ifdef _OPENMP
    if ( length > OPENMP_HOLDS ) use_openmp = TRUE;
#endif
    
    if (use_openmp) {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:value) private(i)
#endif
        for ( i = 0; i < length; ++i ) value += xpt[i] * ypt[i];
    }
    else {
        for ( i = 0; i < length; ++i ) value += xpt[i] * ypt[i];
    }
    
    return value;
}
 
REAL fasp_blas_dvec_relerr (const dvector *x,
                            const dvector *y)
{
    const INT   length = x->row;
    const REAL *xpt = x->val, *ypt = y->val;
    
    SHORT use_openmp = FALSE;
    REAL diff = 0.0, temp = 0.0;
    INT   i;
    
    if ( length != y->row ) {
        printf("### ERROR: Vectors have different dimensions!\n");
        fasp_chkerr(ERROR_DATA_STRUCTURE, __FUNCTION__);
    }
    
#ifdef _OPENMP
    if ( length > OPENMP_HOLDS ) use_openmp = TRUE;
#endif
    
    if ( use_openmp ) {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:temp,diff) private(i)
#endif
        for ( i = 0; i < length; ++i ) {
            temp += xpt[i]*xpt[i];
            diff += pow(xpt[i]-ypt[i],2);
        }
    }
    else {
        for ( i = 0; i < length; ++i ) {
            temp += xpt[i]*xpt[i];
            diff += pow(xpt[i]-ypt[i],2);
        }
    }
    
    return sqrt(diff/temp);
}
 
/*---------------------------------*/
/*--        End of File          --*/
/*---------------------------------*/