#include <math.h>
double
ddot(int n, const double *x, int incX, const double *y, int incY)
{
    int     i;
    double  alpha = 0;
    for (i=0; i<n; ++i) {
        alpha += x[i*incX]*y[i*incY];
    }
    return alpha;
}
void
daxpy(int n, double alpha, const double *x, int incX, double *y, int incY)
{
    int i;
    if (alpha==0) {
        return;
    }
    for (i=0; i<n; ++i) {
        y[i*incY] += alpha*x[i*incX];
    }
}
void
dscal(int n, double alpha, double *x, int incX)
{
    int i;
    if (alpha==1) {
        return;
    }
    for (i=0; i<n; ++i) {
        x[i*incX] *= alpha;
    }
}
void
dcopy(int n, const double *x, int incX, double *y, int incY)
{
    int i;
    for (i=0; i<n; ++i) {
        y[i*incY] = x[i*incX];
    }
}
void
dswap(int n, double *x, int incX, double *y, int incY)
{
    int i;
    for (i=0; i<n; ++i) {
        double tmp;
        tmp = x[i*incX];
        x[i*incX] = y[i*incY];
        y[i*incY] = tmp;
    }
}
double
damax(int n, const double *x, int incX)
{
    int     i;
    double  result = 0;
    for (i=0; i<n; ++i) {
        if (fabs(x[i*incX])>result) {
            result = fabs(x[i*incX]);
        }
    }
    return result;
}
void
dgecopy(int m, int n,
        const double *X, int incRowX, int incColX,
        double *Y, int incRowY, int incColY)
{
    int i, j;
    for (j=0; j<n; ++j) {
        for (i=0; i<m; ++i) {
            Y[i*incRowY+j*incColY] = X[i*incRowX+j*incColX];
        }
    }
}
void
dgescal(int m, int n, double alpha,
        double *X, int incRowX, int incColX)
{
    int i, j;
    for (j=0; j<n; ++j) {
        for (i=0; i<m; ++i) {
            X[i*incRowX+j*incColX] *= alpha;
        }
    }
}
void
dgeaxpy(int m, int n, double alpha,
        const double *X, int incRowX, int incColX,
        double *Y, int incRowY, int incColY)
{
    int i, j;
    if (alpha==0 || m==0 || n==0) {
        return;
    }
    for (j=0; j<n; ++j) {
        for (i=0; i<m; ++i) {
            Y[i*incRowY+j*incColY] += alpha*X[i*incRowX+j*incColX];
        }
    }
}
double
dgenrm1(int m, int n, const double *A, int incRowA, int incColA)
{
    int     i, j;
    double  result = 0;
    for (j=0; j<n; ++j) {
        double sum = 0;
        for (i=0; i<m; ++i) {
            sum += fabs(A[i*incRowA+j*incColA]);
        }
        if (sum>result) {
            result = sum;
        }
    }
    return result;
}
double
dtrnrm1(int m, int n,  int unitDiag, int lower,
        const double *A, int incRowA, int incColA)
{
    int     i, j;
    double  result = 0;
    for (j=0; j<n; ++j) {
        double sum = 0;
        for (i=0; i<m; ++i) {
            if (unitDiag && (i==j)) {
                sum += 1.0;
            } else if ((lower && (i>=j)) || (!lower && (i<=j))) {
                sum += fabs(A[i*incRowA+j*incColA]);
            }
        }
        if (sum>result) {
            result = sum;
        }
    }
    return result;
}