#include <ulmaux.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
#include <sys/times.h>
#include <unistd.h>
double
walltime()
{
struct tms ts;
static double ClockTick=0.0;
if (ClockTick==0.0) {
ClockTick = 1.0 / ((double) sysconf(_SC_CLK_TCK));
}
return ((double) times(&ts)) * ClockTick;
}
void
randDGeMatrix(size_t m, size_t n, bool withNan,
double *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
A[i*incRowA + j*incColA] = withNan
? nan("")
: 2.*(rand()-RAND_MAX/2)/RAND_MAX;
}
}
}
void
printfDGeMatrix(const char * fmt, size_t m, size_t n,
const double *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
printf(fmt, A[i*incRowA + j*incColA]);
}
printf("\n");
}
printf("\n");
}
void
printDGeMatrix(size_t m, size_t n,
const double *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
printfDGeMatrix("%9.2lf ", m, n, A, incRowA, incColA);
}
void
printIGeMatrix(size_t m, size_t n,
const size_t *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
printf("%6zu ", A[i*incRowA + j*incColA]);
}
printf("\n");
}
printf("\n");
}
extern inline size_t
ptrdiff_abs(ptrdiff_t x);
void
dfill_nan(size_t bufsize, double *buf)
{
for (size_t i=0; i<bufsize; ++i) {
buf[i] = nan("");
}
}
void
ifill_rand(size_t bufsize, size_t *buf)
{
for (size_t i=0; i<bufsize; ++i) {
buf[i] = rand();
}
}
void
dgecopy(size_t m, size_t n,
const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
double *B, ptrdiff_t incRowB, ptrdiff_t incColB)
{
if (ptrdiff_abs(incRowB)>ptrdiff_abs(incColB)
&& ptrdiff_abs(incRowA)>ptrdiff_abs(incColA))
{
dgecopy(n, m, A, incColA, incRowA, B, incColB, incRowB);
return;
}
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
B[i*incRowB + j*incColB] = A[i*incRowA + j*incColA];
}
}
}
void
dgeaxpy(size_t m, size_t n, double alpha,
const double *X, ptrdiff_t incRowX, ptrdiff_t incColX,
double *Y, ptrdiff_t incRowY, ptrdiff_t incColY)
{
if (ptrdiff_abs(incRowY)>ptrdiff_abs(incColY)
&& ptrdiff_abs(incRowX)>ptrdiff_abs(incColX))
{
dgeaxpy(n, m, alpha, X, incColX, incRowX, Y, incColY, incRowY);
return;
}
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
Y[i*incRowY + j*incColY] += alpha * X[i*incRowX + j*incColX];
}
}
}
double
dgenrm_inf(size_t m, size_t n,
const double *A,
ptrdiff_t incRowA, ptrdiff_t incColA)
{
double result = 0;
for (size_t i=0; i<m; ++i) {
double sum = 0;
for (size_t j=0; j<n; ++j) {
sum += fabs(A[i*incRowA+j*incColA]);
}
if (isnan(sum)) {
return sum;
}
if (sum>result) {
result = sum;
}
}
return result;
}