# Possible solution

#### Content

• Check out how functions like dgescal handle cases where the matrix is not stored in col major format.

• Observe that function dgenorm_inf accesses elements always row-wise. Make sure that you are aware that traversing elements column-wise would require an additional vector of length $$m$$.

• The test was modified such that each possible storage order for $$A$$ and $$B$$ gets tested!

## Source code

#include <math.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>

//------------------------------------------------------------------------------

#define MY_ABS(x)   ((x)<0 ? -(x) : (x))

//------------------------------------------------------------------------------

void
initMatrix(size_t m, size_t n, double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
bool withNan)
{
// if A is row major initialize A^T
if (MY_ABS(incRowA) > MY_ABS(incColA)) {
initMatrix(n, m, A, incColA, incRowA, withNan);
return;
}
// if A is col major
if (withNan) {
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
A[i*incRowA+j*incColA] = nan("");
}
}
} else {
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
double rValue = ((double)rand() - RAND_MAX/2)*2/RAND_MAX;
A[i*incRowA+j*incColA] = rValue;
}
}
}
}

void
printMatrix(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("%10.3lf ", A[i*incRowA+j*incColA]);
}
printf("\n");
}
printf("\n");
}

//------------------------------------------------------------------------------

void
dgeaxpy(size_t m, size_t n, double alpha,
const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
double *B, ptrdiff_t incRowB, ptrdiff_t incColB)
{
if (m==0 || n==0) {
return;
}
// if B is row major:   B^T <- alpha*A^T + B^T
if (MY_ABS(incRowB) > MY_ABS(incColB)) {
dgeaxpy(n, m, alpha, A, incColA, incRowA, B, incColB, incRowB);
return;
}
// B is col major:
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
B[i*incRowB+j*incColB] += alpha*A[i*incRowA+j*incColA];
}
}
}

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 (m==0 || n==0) {
return;
}
// if B is row major:   B^T <- A^T
if (MY_ABS(incRowB) > MY_ABS(incColB)) {
dgecopy(n, m, A, incColA, incRowA, B, incColB, incRowB);
return;
}
// B is col major:
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
dgescal(size_t m, size_t n, double alpha,
double *A, ptrdiff_t incRowA, ptrdiff_t incColA)
{
if (alpha==1 || m==0 || n==0) {
return;
}
// if A is row major: A^T <- alpha*A^T
if (MY_ABS(incRowA) > MY_ABS(incColA)) {
dgescal(n, m, alpha, A, incColA, incRowA);
return;
}
// A is col major:
if (alpha!=0) {
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
A[i*incRowA+j*incColA] *= alpha;
}
}
} else {
for (size_t j=0; j<n; ++j) {
for (size_t i=0; i<m; ++i) {
A[i*incRowA+j*incColA] = 0;
}
}
}
}

// This operation is not cache friendly!
double
dgenorm_inf(size_t m, size_t n,
const double *A, ptrdiff_t incRowA, ptrdiff_t incColA)
{
double res = 0;
for (size_t i=0; i<m; ++i) {
double asum = 0;
for (size_t j=0; j<n; ++j) {
asum += fabs(A[i*incRowA+j*incColA]);
}
if (asum>res) {
res = asum;
}
}
return res;
}

//------------------------------------------------------------------------------

#ifndef DIM_M
#define DIM_M 4
#endif

#ifndef DIM_N
#define DIM_N 5
#endif

double A[DIM_M*DIM_N];
double B[DIM_M*DIM_N];

const int rowMajorA[] = {0, 1, 0, 1};
const int rowMajorB[] = {0, 0, 1, 1};
const int numTests = sizeof(rowMajorA)/sizeof(int);

int
main()
{
for (int test=0; test<numTests; ++test) {
ptrdiff_t incRowA = rowMajorA[test] ? DIM_N : 1;
ptrdiff_t incColA = rowMajorA[test] ? 1 : DIM_M;

ptrdiff_t incRowB = rowMajorB[test] ? DIM_N : 1;
ptrdiff_t incColB = rowMajorB[test] ? 1 : DIM_M;

printf("Matrix A is stored %s\n", rowMajorA[test] ? "row major"
: "col major");
printf("Matrix B is stored %s\n", rowMajorB[test] ? "row major"
: "col major");
srand(0);
initMatrix(DIM_M, DIM_N, A, incRowA, incColA, false);
initMatrix(DIM_M, DIM_N, B, incRowB, incColB, false);

printf("A =\n");
printMatrix(DIM_M, DIM_N, A, incRowA, incColA);
printf("B =\n");
printMatrix(DIM_M, DIM_N, B, incRowB, incColB);

printf("B = 2*B\n");
dgescal(DIM_M, DIM_N, 2, B, incRowB, incColB);
printf("B =\n");
printMatrix(DIM_M, DIM_N, B, incRowB, incColB);

printf("||A||_inf = %lf\n",
dgenorm_inf(DIM_M, DIM_N, A, incRowA, incColA));
printf("||B||_inf = %lf\n",
dgenorm_inf(DIM_M, DIM_N, B, incRowB, incColB));

dgeaxpy(DIM_M, DIM_N, -1, A, incRowA, incColA, B, incRowB, incColB);
printf("||B-A||_inf = %lf\n",
dgenorm_inf(DIM_M, DIM_N, B, incRowB, incColB));

dgecopy(DIM_M, DIM_N, A, incRowA, incColA, B, incRowB, incColB);
printf("B = A\n");
printf("B =\n");
printMatrix(DIM_M, DIM_N, B, incRowB, incColB);

dgeaxpy(DIM_M, DIM_N, -1, A, incRowA, incColA, B, incRowB, incColB);
printf("||B-A||_inf = %lf\n",
dgenorm_inf(DIM_M, DIM_N, B, incRowB, incColB));
}
}


## Test run

theon$gcc -Wall -std=c99 -O3 norm.c theon$ ./a.out
Matrix A is stored col major
Matrix B is stored col major
A =
-1.000     -0.786     -0.260     -0.655      0.976
0.310      0.033     -0.487     -0.404      0.601
-0.390     -0.021     -0.252      0.287     -0.071
0.350      0.205      0.651      0.579      0.078

B =
0.251      0.959      0.480     -0.918      0.030
-0.500     -0.342     -0.654      0.195     -0.205
0.408     -0.109     -0.969     -0.509     -0.635
0.433      0.417      0.565      0.113      0.292

B = 2*B
B =
0.502      1.918      0.960     -1.835      0.061
-1.000     -0.684     -1.309      0.390     -0.410
0.816     -0.218     -1.937     -1.017     -1.270
0.865      0.834      1.130      0.226      0.584

||A||_inf = 3.676748
||B||_inf = 5.275918
||B-A||_inf = 7.521714
B = A
B =
-1.000     -0.786     -0.260     -0.655      0.976
0.310      0.033     -0.487     -0.404      0.601
-0.390     -0.021     -0.252      0.287     -0.071
0.350      0.205      0.651      0.579      0.078

||B-A||_inf = 0.000000
Matrix A is stored row major
Matrix B is stored col major
A =
-1.000      0.310     -0.390      0.350     -0.786
0.033     -0.021      0.205     -0.260     -0.487
-0.252      0.651     -0.655     -0.404      0.287
0.579      0.976      0.601     -0.071      0.078

B =
0.251      0.959      0.480     -0.918      0.030
-0.500     -0.342     -0.654      0.195     -0.205
0.408     -0.109     -0.969     -0.509     -0.635
0.433      0.417      0.565      0.113      0.292

B = 2*B
B =
0.502      1.918      0.960     -1.835      0.061
-1.000     -0.684     -1.309      0.390     -0.410
0.816     -0.218     -1.937     -1.017     -1.270
0.865      0.834      1.130      0.226      0.584

||A||_inf = 2.837123
||B||_inf = 5.275918
||B-A||_inf = 7.492294
B = A
B =
-1.000      0.310     -0.390      0.350     -0.786
0.033     -0.021      0.205     -0.260     -0.487
-0.252      0.651     -0.655     -0.404      0.287
0.579      0.976      0.601     -0.071      0.078

||B-A||_inf = 0.000000
Matrix A is stored col major
Matrix B is stored row major
A =
-1.000     -0.786     -0.260     -0.655      0.976
0.310      0.033     -0.487     -0.404      0.601
-0.390     -0.021     -0.252      0.287     -0.071
0.350      0.205      0.651      0.579      0.078

B =
0.251     -0.500      0.408      0.433      0.959
-0.342     -0.109      0.417      0.480     -0.654
-0.969      0.565     -0.918      0.195     -0.509
0.113      0.030     -0.205     -0.635      0.292

B = 2*B
B =
0.502     -1.000      0.816      0.865      1.918
-0.684     -0.218      0.834      0.960     -1.309
-1.937      1.130     -1.835      0.390     -1.017
0.226      0.061     -0.410     -1.270      0.584

||A||_inf = 3.676748
||B||_inf = 6.309641
||B-A||_inf = 5.839961
B = A
B =
-1.000     -0.786     -0.260     -0.655      0.976
0.310      0.033     -0.487     -0.404      0.601
-0.390     -0.021     -0.252      0.287     -0.071
0.350      0.205      0.651      0.579      0.078

||B-A||_inf = 0.000000
Matrix A is stored row major
Matrix B is stored row major
A =
-1.000      0.310     -0.390      0.350     -0.786
0.033     -0.021      0.205     -0.260     -0.487
-0.252      0.651     -0.655     -0.404      0.287
0.579      0.976      0.601     -0.071      0.078

B =
0.251     -0.500      0.408      0.433      0.959
-0.342     -0.109      0.417      0.480     -0.654
-0.969      0.565     -0.918      0.195     -0.509
0.113      0.030     -0.205     -0.635      0.292

B = 2*B
B =
0.502     -1.000      0.816      0.865      1.918
-0.684     -0.218      0.834      0.960     -1.309
-1.937      1.130     -1.835      0.390     -1.017
0.226      0.061     -0.410     -1.270      0.584

||A||_inf = 2.837123
||B||_inf = 6.309641
||B-A||_inf = 7.238380
B = A
B =
-1.000      0.310     -0.390      0.350     -0.786
0.033     -0.021      0.205     -0.260     -0.487
-0.252      0.651     -0.655     -0.404      0.287
0.579      0.976      0.601     -0.071      0.078

||B-A||_inf = 0.000000
theon\$