#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;
}