exponential_linear_algebra.h

Go to the documentation of this file.

#ifndef EXPONENTIAL_LINEAR_ALGEBRA_H
#define EXPONENTIAL_LINEAR_ALGEBRA_H

#include "header.h"
#include "solver_options.h"
#include "solver_props.h"



static inline
void matvec_m_by_m (const int m, const double * __restrict__ A,
                    const double * __restrict__ V, double * __restrict__ Av) {
    //for each row

    for (int i = 0; i < m; ++i) {
        Av[i] = A[i] * V[0];

        //go across a row of A, multiplying by a column
        for (int j = 1; j < m; ++j) {
            Av[i] += A[j * STRIDE + i] * V[j];
        }
    }
}


static inline void matvec_m_by_m_plusequal (const int m, const double * __restrict__ A,
                                            const double * __restrict__ V, double * __restrict__ Av)
{
    //for each row

    for (int i = 0; i < m; ++i) {
        Av[i] = A[i] * V[0];

        //go across a row of A, multiplying by a column of phiHm

        for (int j = 1; j < m; ++j) {
            Av[i] += A[j * STRIDE + i] * V[j];
        }

        Av[i] += V[i];
    }
}

static inline
void matvec_n_by_m_scale (const int m, const double scale,
                          const double * __restrict__ A, const double * __restrict__ V,
                          double * __restrict__ Av) {
    //for each row

    for (int i = 0; i < NSP; ++i) {
        Av[i] = A[i] * V[0];

        //go across a row of A, multiplying by a column of phiHm

        for (int j = 1; j < m; ++j) {
            Av[i] += A[j * NSP + i] * V[j];
        }

        Av[i] *= scale;
    }
}

static inline
void matvec_n_by_m_scale_special (const int m, const double* __restrict__ scale,
                                  const double * __restrict__ A, const double** __restrict__ V,
                                  double** __restrict__ Av) {
    //for each row
    for (int i = 0; i < NSP; ++i) {
        Av[0][i] = A[i] * V[0][0];
        Av[1][i] = A[i] * V[1][0];
        Av[2][i] = A[i] * V[2][0];

        //go across a row of A, multiplying by a column of phiHm

        for (int j = 1; j < m; ++j) {
            Av[0][i] += A[j * NSP + i] * V[0][j];
            Av[1][i] += A[j * NSP + i] * V[1][j];
            Av[2][i] += A[j * NSP + i] * V[2][j];
        }

        Av[0][i] *= scale[0];
        Av[1][i] *= scale[1];
        Av[2][i] = Av[2][i] * scale[2] + V[3][i] + V[4][i];
    }
}

static inline
void matvec_n_by_m_scale_special2 (const int m, const double* __restrict__ scale,
                                   const double* __restrict__ A, const double** __restrict__ V,
                                   double** __restrict__ Av) {
    //for each row
    for (int i = 0; i < NSP; ++i) {
        Av[0][i] = A[i] * V[0][0];
        Av[1][i] = A[i] * V[1][0];

        //go across a row of A, multiplying by a column of phiHm

        for (int j = 1; j < m; ++j) {
            Av[0][i] += A[j * NSP + i] * V[0][j];
            Av[1][i] += A[j * NSP + i] * V[1][j];
        }

        Av[0][i] *= scale[0];
        Av[1][i] *= scale[1];
    }
}


static inline
void matvec_n_by_m_scale_add (const int m, const double scale, const double* __restrict__ A,
                              const double* __restrict__ V, double* __restrict__ Av,
                              const double* __restrict__ add) {
    //for each row
    for (int i = 0; i < NSP; ++i) {
        Av[i] = A[i] * V[0];

        //go across a row of A, multiplying by a column of phiHm

        for (int j = 1; j < m; ++j) {
            Av[i] += A[j * NSP + i] * V[j];
        }

        Av[i] = Av[i] * scale + add[i];
    }
}


static inline
void matvec_n_by_m_scale_add_subtract (const int m, const double scale,
                                       const double* __restrict__ A, const double* __restrict__ V,
                                       double* __restrict__ Av, const double* __restrict__ add,
                                       const double* __restrict__ sub) {
    //for each row

    for (int i = 0; i < NSP; ++i) {
        Av[i] = A[i] * V[0];

        //go across a row of A, multiplying by a column of phiHm

        for (int j = 1; j < m; ++j) {
            Av[i] += A[j * NSP + i] * V[j];
        }

        Av[i] = Av[i] * scale + 2.0 * add[i] - sub[i];
    }
}


static inline
void scale (const double* __restrict__ y0, const double* __restrict__ y1, double* __restrict__ sc) {

    for (int i = 0; i < NSP; ++i) {
        sc[i] = 1.0 / (ATOL + fmax(fabs(y0[i]), fabs(y1[i])) * RTOL);
    }
}


static inline
void scale_init (const double* __restrict__ y0, double* __restrict__ sc) {

    for (int i = 0; i < NSP; ++i) {
        sc[i] = 1.0 / (ATOL + fabs(y0[i]) * RTOL);
    }
}


static inline
double sc_norm (const double* __restrict__ nums, const double* __restrict__ sc) {
    double norm = 0.0;

    for (int i = 0; i < NSP; ++i) {
        norm += nums[i] * nums[i] * sc[i] * sc[i];
    }

    return sqrt(norm / NSP);
}

static inline
double two_norm(const double* __restrict__ v)
{
    double norm = 0.0;

    for (int i = 0; i < NSP; ++i) {
        norm += v[i] * v[i];
    }
    return sqrt(norm);
}

static inline
double normalize (const double* __restrict__ v, double* __restrict__ v_out) {

    double norm = two_norm(v);

    if (norm == 0)
        norm = 1;

    double m_norm = 1.0 / norm;


    for (int i = 0; i < NSP; ++i) {
        v_out[i] = v[i] * m_norm;
    }
    return norm;
}

static inline
double dotproduct(const double* __restrict__ w, const double* __restrict__ Vm)
{
    double sum = 0;

    for(int i = 0; i < NSP; i++)
    {
        sum += w[i] * Vm[i];
    }
    return sum;
}

static inline void scale_subtract(const double s, const double* __restrict__ Vm, double* __restrict__ w)
{

    for (int i = 0; i < NSP; i++)
    {
        w[i] -= s * Vm[i];
    }
}


static inline void scale_mult(const double s, const double* __restrict__ w, double* __restrict__ Vm)
{

    for (int i = 0; i < NSP; i++)
    {
        Vm[i] = w[i] * s;
    }
}
#endif