AuxVector.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       --*/
/*---------------------------------*/
 
SHORT fasp_dvec_isnan (const dvector *u)
{
    INT i;
    
    for ( i = 0; i < u->row; i++ ) {
        if ( isnan(u->val[i]) ) return TRUE;
    }
 
    return FALSE;
}
 
dvector fasp_dvec_create (const INT m)
{
    dvector u;
    
    u.row = m;
    u.val = (REAL *)fasp_mem_calloc(m,sizeof(REAL));
    
    return u;
}
 
ivector fasp_ivec_create (const INT m)
{    
    ivector u;
    
    u.row = m;
    u.val = (INT *)fasp_mem_calloc(m,sizeof(INT)); 
    
    return u;
}
 
void fasp_dvec_alloc (const INT  m,
                      dvector   *u)
{    
    u->row = m;
    u->val = (REAL*)fasp_mem_calloc(m,sizeof(REAL)); 
    
    return;
}
 
void fasp_ivec_alloc (const INT  m,
                      ivector   *u)
{    
    
    u->row = m;
    u->val = (INT*)fasp_mem_calloc(m,sizeof(INT));
    
    return;
}
 
void fasp_dvec_free (dvector *u)
{    
    if ( u == NULL ) return;
    
    fasp_mem_free(u->val); u->val = NULL; u->row = 0;
}
 
void fasp_ivec_free (ivector *u)
{    
    if ( u == NULL ) return;
    
    fasp_mem_free(u->val); u->val = NULL; u->row = 0;
}
 
void fasp_dvec_rand (const INT  n,
                     dvector   *x)
{
    const INT va = 0;
    const INT vb = n;
    
    INT s=1, i,j;
    
    srand(s);
    for ( i = 0; i < n; ++i ) {
        j = 1 + (INT) (((REAL)n)*rand()/(RAND_MAX+1.0));
        x->val[i] = (((REAL)j)-va)/(vb-va);
    }
    x->row = n;
}
 
void fasp_dvec_set (INT         n,
                    dvector    *x,
                    const REAL  val)
{
    INT   i;
    REAL *xpt = x->val;
    
    if ( n > 0 ) x->row = n;
    else n = x->row;
    
#ifdef _OPENMP 
    // variables for OpenMP 
    INT myid, mybegin, myend;
    INT nthreads = fasp_get_num_threads();
#endif
   
    if (val == 0.0) {
        
#ifdef _OPENMP 
        if (n > OPENMP_HOLDS) {
#pragma omp parallel for private(myid, mybegin, myend) 
            for (myid = 0; myid < nthreads; myid++ ) {
                fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
                memset(&xpt[mybegin], 0x0, sizeof(REAL)*(myend-mybegin));
            }
        }
        else {
#endif
            memset(xpt, 0x0, sizeof(REAL)*n);
#ifdef _OPENMP
        }
#endif
        
    }
    
    else {
        
#ifdef _OPENMP 
        if (n > OPENMP_HOLDS) {
#pragma omp parallel for private(myid, mybegin, myend) 
            for (myid = 0; myid < nthreads; myid++ ) {
                fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
                for (i=mybegin; i<myend; ++i) xpt[i]=val;
            }
        }
        else {
#endif
            for (i=0; i<n; ++i) xpt[i]=val;
#ifdef _OPENMP
        }
#endif
        
    }
}
 
void fasp_ivec_set (INT        n,
                    ivector   *u,
                    const INT  m)
{    
    SHORT nthreads = 1, use_openmp = FALSE;
    INT   i;
    
    if ( n > 0 ) u->row = n;
    else n = u->row;
 
#ifdef _OPENMP
    if ( n > OPENMP_HOLDS ) {
        use_openmp = TRUE;
        nthreads = fasp_get_num_threads();
    }
#endif
    
    if (use_openmp) {
        INT mybegin, myend, myid;
#ifdef _OPENMP 
#pragma omp parallel for private(myid, mybegin, myend, i) 
#endif
        for (myid = 0; myid < nthreads; myid++ ) {
            fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
            for (i=mybegin; i<myend; ++i) u->val[i] = m;
        }        
    }    
    else {
        for (i=0; i<n; ++i) u->val[i] = m;
    }
}
 
void fasp_dvec_cp (const dvector  *x,
                   dvector        *y)
{
    y->row = x->row;
    memcpy(y->val,x->val,x->row*sizeof(REAL));
}
 
REAL fasp_dvec_maxdiff (const dvector *x,
                        const dvector *y)
{
    const INT length = x->row;
    const REAL *xpt = x->val, *ypt = y->val;
    
    SHORT use_openmp = FALSE;
    INT   i;
    REAL  Linf = 0.0, diffi = 0.0;
 
#ifdef _OPENMP 
    INT myid, mybegin, myend, nthreads;
    if ( length > OPENMP_HOLDS ) {
        use_openmp = TRUE;
        nthreads = fasp_get_num_threads();
    }
#endif
    
    if(use_openmp) {
#ifdef _OPENMP 
        REAL temp = 0.;
#pragma omp parallel firstprivate(temp) private(myid, mybegin, myend, i, diffi) 
        {
            myid = omp_get_thread_num();
            fasp_get_start_end(myid, nthreads, length, &mybegin, &myend);
            for(i=mybegin; i<myend; i++) {
                if ((diffi = ABS(xpt[i]-ypt[i])) > temp) temp = diffi;
            }
#pragma omp critical
            if(temp > Linf) Linf = temp;
        }
#endif
    }
    else {
        for (i=0; i<length; ++i) {
            if ((diffi = ABS(xpt[i]-ypt[i])) > Linf) Linf = diffi;
        }
    }
    
    return Linf;
}
 
void fasp_dvec_symdiagscale (dvector        *b,
                             const dvector  *diag)
{
    // information about dvector
    const INT    n = b->row;
    REAL      *val = b->val;
    
    // local variables
    SHORT use_openmp = FALSE;
    INT   i;
    
    if ( diag->row != n ) {
        printf("### ERROR: Sizes of diag = %d != dvector = %d!", diag->row, n);
        fasp_chkerr(ERROR_MISC, __FUNCTION__);
    }
    
#ifdef _OPENMP 
    INT mybegin, myend, myid, nthreads;
    if ( n > OPENMP_HOLDS ){
        use_openmp = TRUE;
        nthreads = fasp_get_num_threads();
    }
#endif
 
    if (use_openmp) {
#ifdef _OPENMP 
#pragma omp parallel for private(myid, mybegin,myend) 
        for (myid = 0; myid < nthreads; myid++ ) {
            fasp_get_start_end(myid, nthreads, n, &mybegin, &myend);
            for (i=mybegin; i<myend; ++i) val[i] = val[i]/sqrt(diag->val[i]);
        }        
#endif
    }    
    else {
        for (i=0; i<n; ++i) val[i] = val[i]/sqrt(diag->val[i]);
    }
    
    return;
}
 
/*---------------------------------*/
/*--        End of File          --*/
/*---------------------------------*/