1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
#include <cstdio>
#include <cstdlib>
struct DGeMatrix
{
double *data;
long numRows, numCols;
};
struct DGeMatrixView
{
double *data;
long numRows, numCols;
long leadingDimension;
};
// Constructor for DGeMatrix
void
DGeMatrix_alloc(DGeMatrix *A, long m, long n)
{
A->data = (double *)malloc(m*n*sizeof(double));
A->numRows = m;
A->numCols = n;
}
// Destructor for DGeMatrix
void
DGeMatrix_free(DGeMatrix *A)
{
free(A->data);
A->data = 0;
}
// View from DGematrix
void
DGeMatrix_view(const DGeMatrix *A, DGeMatrixView *B,
long fromRow, long toRow,
long fromCol, long toCol)
{
long ldA = A->numRows;
B->data = A->data + fromRow + fromCol*ldA;
B->numRows = toRow-fromRow+1;
B->numCols = toCol-fromCol+1;
B->leadingDimension = ldA;
}
// Change entry in DGeMatrix
void
DGeMatrix_set(const DGeMatrix *A, long row, long col, double value)
{
A->data[row+col*A->numRows] = value;
}
// Get entry from DGeMatrix
double
DGeMatrix_get(const DGeMatrix *A, long row, long col)
{
return A->data[row+col*A->numRows];
}
// Change entry in DGeMatrixView
void
DGeMatrixView_set(const DGeMatrixView *A, long row, long col, double value)
{
A->data[row+col*A->leadingDimension] = value;
}
// Get entry from DGeMatrixView
double
DGeMatrixView_get(const DGeMatrixView *A, long row, long col)
{
return A->data[row+col*A->leadingDimension];
}
int
main()
{
// DGeMatrix A(3,4);
DGeMatrix A;
DGeMatrix_alloc(&A, 3, 4);
// A = 1, 2, 3, 4,
// 5, 6, 7, 8,
// 9, 10, 11, 12;
//
// The list initializer always fills the matrix row wise. Even if col wise
// would be faster. But anyway, the sole purpose of the list initializer is
// for writing compact examples in the tutorial
//
long counter = 1;
for (long j=0; j<A.numCols; ++j) {
for (long i=0; i<A.numRows; ++i) {
DGeMatrix_set(&A, i, j, counter);
++counter;
}
}
// cout << "A = " << A << endl;
printf("A =\n");
for (long i=0; i<A.numRows; ++i) {
for (long j=0; j<A.numCols; ++j) {
printf("%8.3lf", DGeMatrix_get(&A, i, j));
}
printf("\n");
}
printf("\n");
// DGeMatrixView B = A(_(2,3),_);
DGeMatrixView B;
DGeMatrix_view(&A, &B, 1, 2, 0, 3);
// cout << "B = " << B << endl;
printf("B =\n");
for (long i=0; i<B.numRows; ++i) {
for (long j=0; j<B.numCols; ++j) {
printf("%8.3lf", DGeMatrixView_get(&B, i, j));
}
printf("\n");
}
printf("\n");
// B(1,2) = 42;
DGeMatrixView_set(&B, 0, 1, 42);
// cout << "Changed B(1,2)." << std::endl;
printf("Changed B(1,2).\n");
// cout << "A = " << A << endl;
printf("A =\n");
for (long i=0; i<A.numRows; ++i) {
for (long j=0; j<A.numCols; ++j) {
printf("%8.3lf", DGeMatrix_get(&A, i, j));
}
printf("\n");
}
printf("\n");
// cout << "B = " << B << endl;
printf("B =\n");
for (long i=0; i<B.numRows; ++i) {
for (long j=0; j<B.numCols; ++j) {
printf("%8.3lf", DGeMatrixView_get(&B, i, j));
}
printf("\n");
}
printf("\n");
// Destructor for A
DGeMatrix_free(&A);
}
|