#include <stdio.h>
#include <math.h>
void
initGeMatrix(int m, int n, double *A, int incRowA, int incColA)
{
int i, j;
for (i=0; i<m; ++i) {
for (j=0; j<n; ++j) {
A[i*incRowA+j*incColA] = i*n + j + 1;
}
}
}
void
printGeMatrix(int m, int n, const double *A, int incRowA, int incColA)
{
int i, j;
for (i=0; i<m; ++i) {
for (j=0; j<n; ++j) {
printf("%10.4lf ", A[i*incRowA+j*incColA]);
}
printf("\n");
}
printf("\n\n");
}
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;
}
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;
}
#ifndef DIM_M
#define DIM_M 5
#endif
#ifndef DIM_N
#define DIM_N 5
#endif
#ifndef ROWMAJOR
#define ROWMAJOR 0
#endif
#if (ROWMAJOR==1)
# define INCROW_A DIM_N
# define INCCOL_A 1
#else
# define INCROW_A 1
# define INCCOL_A DIM_M
#endif
double A[DIM_M*DIM_N];
int
main()
{
initGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);
printf("A =\n");
printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);
printf("A(:,1) * A(:,2) = %lf\n",
ddot(DIM_M, &A[0*INCCOL_A], INCROW_A, &A[1*INCCOL_A], INCROW_A));
daxpy(DIM_M, 2.0, &A[1*INCCOL_A], INCROW_A, &A[2*INCCOL_A], INCROW_A);
printf("A(:,3) = A(:3) + 2*A(:,2)\n");
printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);
dscal(DIM_M, 0.5, &A[3*INCCOL_A], INCROW_A);
printf("A(:,4) = 0.5*A(:4)\n");
printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);
printf("inf-nrm(A(:,1)) = %lf\n", damax(DIM_M, &A[0*INCCOL_A], INCROW_A));
printf("inf-nrm(A(2,:)) = %lf\n", damax(DIM_N, &A[1*INCROW_A], INCCOL_A));
printf("nrm1(A(1:4,1:5)) = %lf\n",
dgenrm1(4, 5, A, INCROW_A, INCCOL_A));
printf("nrm1(A(2:4,3:5)) = %lf\n",
dgenrm1(3, 3, &A[1*INCROW_A+2*INCCOL_A], INCROW_A, INCCOL_A));
dswap(DIM_M, &A[1*INCCOL_A], INCROW_A, &A[2*INCCOL_A], INCROW_A);
printf("swap(A(:,2), A(:,3))\n");
printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);
dcopy(DIM_M, &A[2*INCCOL_A], INCROW_A, &A[3*INCCOL_A], INCROW_A);
printf("A(:,4) = A(:,3)\n");
printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);
return 0;
}