Example Code (Matrix Operations)#

Use a simple C array for vector (pod01_vector.cpp)#
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <iostream>
#include <iomanip>

int main(int argc, char ** argv)
{
    constexpr size_t width = 5;

    double vector[width];

    // Populate a vector.
    for (size_t i=0; i<width; ++i)
    {
        vector[i] = i;
    }

    std::cout << "vector elements in memory:" << std::endl << " ";
    for (size_t i=0; i<width; ++i)
    {
        std::cout << " " << vector[i];
    }
    std::cout << std::endl;

    return 0;
}
Build “Use a simple C array for vector (pod01_vector.cpp)”#
$ g++ pod01_vector.cpp -o pod01_vector -std=c++17 -O3 -g -m64
Use a 2D array on stack for a square matrix (pod02_matrix_auto.cpp)#
 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
#include <iostream>
#include <iomanip>

int main(int argc, char ** argv)
{
    constexpr size_t width = 5;

    double amatrix[width][width];

    // Populate the matrix on stack (row-major 2D array).
    for (size_t i=0; i<width; ++i) // the i-th row
    {
        for (size_t j=0; j<width; ++j) // the j-th column
        {
            amatrix[i][j] = i*10 + j;
        }
    }

    std::cout << "2D array elements:";
    for (size_t i=0; i<width; ++i) // the i-th row
    {
        std::cout << std::endl << " ";
        for (size_t j=0; j<width; ++j) // the j-th column
        {
            std::cout << " " << std::setfill('0') << std::setw(2)
                      << amatrix[i][j];
        }
    }
    std::cout << std::endl;

    return 0;
}
Build “Use a 2D array on stack for a square matrix (pod02_matrix_auto.cpp)”#
$ g++ pod02_matrix_auto.cpp -o pod02_matrix_auto -std=c++17 -O3 -g -m64
C++ does not support variable-length arrays (pod_bad_matrix.cpp)#
 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
#include <iostream>
#include <iomanip>

void work(double * buffer, size_t width)
{
    // This should not work since width is unknown in compile time.
    double (*matrix)[width] = reinterpret_cast<double (*)[width]>(buffer);

    for (size_t i=0; i<width; ++i) // the i-th row
    {
        for (size_t j=0; j<width; ++j) // the j-th column
        {
            matrix[i][j] = i*10 + j;
        }
    }

    std::cout << "matrix:";
    for (size_t i=0; i<width; ++i)
    {
        std::cout << std::endl << " ";
        for (size_t j=0; j<width; ++j)
        {
            std::cout << " " << std::setfill('0') << std::setw(2)
                      << matrix[i][j];
        }
    }
    std::cout << std::endl;
}

int main(int argc, char ** argv)
{
    size_t width = 5;

    double * buffer = new double[width*width];

    work(buffer, width);

    delete[] buffer;

    return 0;
}
Row-majored 2D array (pod03_matrix_rowmajor.cpp)#
 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
#include <iostream>
#include <iomanip>

int main(int argc, char ** argv)
{
    constexpr size_t width = 5;

    double * buffer = new double[width*width];
    std::cout << "buffer address: " << buffer << std::endl;

    // Populate a buffer (row-major 2D array).
    for (size_t i=0; i<width; ++i) // the i-th row
    {
        for (size_t j=0; j<width; ++j) // the j-th column
        {
            buffer[i*width + j] = i*10 + j;
        }
    }

    // Make a pointer to double[width].  Note width is a constexpr.
    double (*matrix)[width] = reinterpret_cast<double (*)[width]>(buffer);
    std::cout << "matrix address: " << matrix << std::endl;

    std::cout << "matrix (row-major) elements as 2D array:";
    for (size_t i=0; i<width; ++i) // the i-th row
    {
        std::cout << std::endl << " ";
        for (size_t j=0; j<width; ++j) // the j-th column
        {
            std::cout << " " << std::setfill('0') << std::setw(2)
                      << matrix[i][j];
        }
    }
    std::cout << std::endl;

    std::cout << "matrix (row-major) elements in memory:" << std::endl << " ";
    for (size_t i=0; i<width*width; ++i)
    {
        std::cout << " " << std::setfill('0') << std::setw(2) << buffer[i];
    }
    std::cout << std::endl;
    std::cout << "row majoring: "
              << "the fastest moving index is the trailing index"
              << std::endl;

    delete[] buffer;

    return 0;
}
Build “Row-majored 2D array (pod03_matrix_rowmajor.cpp)”#
$ g++ pod03_matrix_rowmajor.cpp -o pod03_matrix_rowmajor -std=c++17 -O3 -g -m64
Column-majored 2D array (pod04_matrix_colmajor.cpp)#
 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
#include <iostream>
#include <iomanip>

int main(int argc, char ** argv)
{
    constexpr size_t width = 5;

    double * buffer = new double[width*width];
    std::cout << "buffer address: " << buffer << std::endl;

    // Populate a buffer (column-major 2D array).
    for (size_t i=0; i<width; ++i) // the i-th row
    {
        for (size_t j=0; j<width; ++j) // the j-th column
        {
            buffer[j*width + i] = i*10 + j;
        }
    }

    // Make a pointer to double[width].  Note width is a constexpr.
    double (*matrix)[width] = reinterpret_cast<double (*)[width]>(buffer);
    std::cout << "matrix address: " << matrix << std::endl;

    std::cout << "matrix (column-major) elements as 2D array:";
    for (size_t i=0; i<width; ++i) // the i-th row
    {
        std::cout << std::endl << " ";
        for (size_t j=0; j<width; ++j) // the j-th column
        {
            std::cout << " " << std::setfill('0') << std::setw(2)
                      << matrix[j][i];
        }
    }
    std::cout << std::endl;

    std::cout << "matrix (column-major) elements in memory:" << std::endl << " ";
    for (size_t i=0; i<width*width; ++i)
    {
        std::cout << " " << std::setfill('0') << std::setw(2) << buffer[i];
    }
    std::cout << std::endl;
    std::cout << "column majoring: "
              << "the fastest moving index is the leading index"
              << std::endl;

    delete[] buffer;

    return 0;
}
Build “Column-majored 2D array (pod04_matrix_colmajor.cpp)”#
$ g++ pod04_matrix_colmajor.cpp -o pod04_matrix_colmajor -std=c++17 -O3 -g -m64
Skeleton implementation for a C++ matrix class (ma01_matrix_class.cpp)#
 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
#include <iostream>
#include <iomanip>

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol)
      : m_nrow(nrow), m_ncol(ncol)
    {
        size_t nelement = nrow * ncol;
        m_buffer = new double[nelement];
    }

    // TODO: copy and move constructors and assignment operators.

    ~Matrix()
    {
        delete[] m_buffer;
    }

    // No bound check.
    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[row*m_ncol + col];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[row*m_ncol + col];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

private:

    size_t m_nrow;
    size_t m_ncol;
    double * m_buffer;

};

/**
 * Populate the matrix object.
 */
void populate(Matrix & matrix)
{
    for (size_t i=0; i<matrix.nrow(); ++i) // the i-th row
    {
        for (size_t j=0; j<matrix.ncol(); ++j) // the j-th column
        {
            matrix(i, j) = i*10 + j;
        }
    }
}

int main(int argc, char ** argv)
{
    size_t width = 5;

    Matrix matrix(width, width);

    populate(matrix);

    std::cout << "matrix:";
    for (size_t i=0; i<matrix.nrow(); ++i) // the i-th row
    {
        std::cout << std::endl << " ";
        for (size_t j=0; j<matrix.ncol(); ++j) // the j-th column
        {
            std::cout << " " << std::setfill('0') << std::setw(2)
                      << matrix(i, j);
        }
    }
    std::cout << std::endl;

    return 0;
}
Build “Skeleton implementation for a C++ matrix class (ma01_matrix_class.cpp)”#
$ g++ ma01_matrix_class.cpp -o ma01_matrix_class -std=c++17 -O3 -g -m64
Example code for matrix-vector multiplication (ma02_matrix_vector.cpp)#
  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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol)
      : m_nrow(nrow), m_ncol(ncol)
    {
        reset_buffer(nrow, ncol);
    }

    Matrix(Matrix const & other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
    {
        reset_buffer(other.m_nrow, other.m_ncol);
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
    }

    Matrix & operator=(Matrix const & other)
    {
        if (this == &other) { return *this; }
        if (m_nrow != other.m_nrow || m_ncol != other.m_ncol)
        {
            reset_buffer(other.m_nrow, other.m_ncol);
        }
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
        return *this;
    }

    // TODO: move constructors and assignment operators.

    ~Matrix()
    {
        reset_buffer(0, 0);
    }

    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[index(row, col)];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[index(row, col)];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

    size_t size() const { return m_nrow * m_ncol; }
    double buffer(size_t i) const { return m_buffer[i]; }
    std::vector<double> buffer_vector() const
    {
        return std::vector<double>(m_buffer, m_buffer+size());
    }

    bool is_transposed() const { return m_transpose; }

    Matrix & transpose()
    {
        m_transpose = !m_transpose;
        std::swap(m_nrow, m_ncol);
        return *this;
    }

private:

    size_t index(size_t row, size_t col) const
    {
        if (m_transpose) { return row          + col * m_nrow; }
        else             { return row * m_ncol + col         ; }
    }

    void reset_buffer(size_t nrow, size_t ncol)
    {
        if (m_buffer) { delete[] m_buffer; }
        const size_t nelement = nrow * ncol;
        if (nelement) { m_buffer = new double[nelement]; }
        else          { m_buffer = nullptr; }
        m_nrow = nrow;
        m_ncol = ncol;
    }

    size_t m_nrow = 0;
    size_t m_ncol = 0;
    bool m_transpose = false;
    double * m_buffer = nullptr;

};

/*
 * Naive matrix vector multiplication.
 */
std::vector<double> operator*(Matrix const & mat, std::vector<double> const & vec)
{
    if (mat.ncol() != vec.size())
    {
        throw std::out_of_range("matrix column differs from vector size");
    }

    std::vector<double> ret(mat.nrow());

    for (size_t i=0; i<mat.nrow(); ++i) // the i-th row
    {
        double v = 0;
        for (size_t j=0; j<mat.ncol(); ++j) // the j-th column
        {
            v += mat(i,j) * vec[j];
        }
        ret[i] = v;
    }

    return ret;
}

std::ostream & operator << (std::ostream & ostr, Matrix const & mat)
{
    for (size_t i=0; i<mat.nrow(); ++i) // the i-th row
    {
        ostr << std::endl << " ";
        for (size_t j=0; j<mat.ncol(); ++j) // the j-th column
        {
            ostr << " " << std::setw(2) << mat(i, j);
        }
    }

    return ostr;
}

std::ostream & operator << (std::ostream & ostr, std::vector<double> const & vec)
{
    for (size_t i=0; i<vec.size(); ++i)
    {
        std::cout << " " << vec[i];
    }

    return ostr;
}

int main(int argc, char ** argv)
{
    size_t width = 5;

    std::cout << ">>> square matrix-vector multiplication:" << std::endl;
    Matrix mat(width, width);

    for (size_t i=0; i<mat.nrow(); ++i) // the i-th row
    {
        for (size_t j=0; j<mat.ncol(); ++j) // the j-th column
        {
            mat(i, j) = i == j ? 1 : 0;
        }
    }

    std::vector<double> vec{1, 0, 0, 0, 0};
    std::vector<double> res = mat * vec;

    std::cout << "matrix A:" << mat << std::endl;
    std::cout << "vector b:" << vec << std::endl;
    std::cout << "A*b =" << res << std::endl;

    std::cout << ">>> m*n matrix-vector multiplication:" << std::endl;
    Matrix mat2(2, 3);

    double v = 1;
    for (size_t i=0; i<mat2.nrow(); ++i) // the i-th row
    {
        for (size_t j=0; j<mat2.ncol(); ++j) // the j-th column
        {
            mat2(i, j) = v;
            v += 1;
        }
    }

    std::vector<double> vec2{1, 2, 3};
    std::vector<double> res2 = mat2 * vec2;

    std::cout << "matrix A:" << mat2 << std::endl;
    std::cout << "vector b:" << vec2 << std::endl;
    std::cout << "A*b =" << res2 << std::endl;

    std::cout << ">>> transposed matrix-vector multiplication:" << std::endl;
    mat2.transpose();
    std::vector<double> vec3{1, 2};
    std::vector<double> res3 = mat2 * vec3;

    std::cout << "matrix A:" << mat2 << std::endl;
    std::cout << "matrix A buffer:" << mat2.buffer_vector() << std::endl;
    std::cout << "vector b:" << vec3 << std::endl;
    std::cout << "A*b =" << res3 << std::endl;

    std::cout << ">>> copied transposed matrix-vector multiplication:" << std::endl;
    Matrix mat3 = mat2;
    res3 = mat3 * vec3;

    std::cout << "matrix A:" << mat3 << std::endl;
    std::cout << "matrix A buffer:" << mat3.buffer_vector() << std::endl;
    std::cout << "vector b:" << vec3 << std::endl;
    std::cout << "A*b =" << res3 << std::endl;

    return 0;
}
Build “Example code for matrix-vector multiplication (ma02_matrix_vector.cpp)”#
$ g++ ma02_matrix_vector.cpp -o ma02_matrix_vector -std=c++17 -O3 -g -m64
Example code for matrix-matrix multiplication (ma03_matrix_matrix.cpp)#
  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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol)
      : m_nrow(nrow), m_ncol(ncol)
    {
        reset_buffer(nrow, ncol);
    }

    Matrix(size_t nrow, size_t ncol, std::vector<double> const & vec)
      : m_nrow(nrow), m_ncol(ncol)
    {
        reset_buffer(nrow, ncol);
        (*this) = vec;
    }

    Matrix & operator=(std::vector<double> const & vec)
    {
        if (size() != vec.size())
        {
            throw std::out_of_range("number of elements mismatch");
        }

        size_t k = 0;
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = vec[k];
                ++k;
            }
        }

        return *this;
    }

    Matrix(Matrix const & other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
    {
        reset_buffer(other.m_nrow, other.m_ncol);
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
    }

    Matrix & operator=(Matrix const & other)
    {
        if (this == &other) { return *this; }
        if (m_nrow != other.m_nrow || m_ncol != other.m_ncol)
        {
            reset_buffer(other.m_nrow, other.m_ncol);
        }
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
        return *this;
    }

    Matrix(Matrix && other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
    {
        reset_buffer(0, 0);
        std::swap(m_nrow, other.m_nrow);
        std::swap(m_ncol, other.m_ncol);
        std::swap(m_buffer, other.m_buffer);
    }

    Matrix & operator=(Matrix && other)
    {
        if (this == &other) { return *this; }
        reset_buffer(0, 0);
        std::swap(m_nrow, other.m_nrow);
        std::swap(m_ncol, other.m_ncol);
        std::swap(m_buffer, other.m_buffer);
        return *this;
    }

    ~Matrix()
    {
        reset_buffer(0, 0);
    }

    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[index(row, col)];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[index(row, col)];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

    size_t size() const { return m_nrow * m_ncol; }
    double buffer(size_t i) const { return m_buffer[i]; }
    std::vector<double> buffer_vector() const
    {
        return std::vector<double>(m_buffer, m_buffer+size());
    }

private:

    size_t index(size_t row, size_t col) const
    {
        return row + col * m_nrow;
    }

    void reset_buffer(size_t nrow, size_t ncol)
    {
        if (m_buffer) { delete[] m_buffer; }
        const size_t nelement = nrow * ncol;
        if (nelement) { m_buffer = new double[nelement]; }
        else          { m_buffer = nullptr; }
        m_nrow = nrow;
        m_ncol = ncol;
    }

    size_t m_nrow = 0;
    size_t m_ncol = 0;
    double * m_buffer = nullptr;

};

/*
 * Naive matrix matrix multiplication.
 */
Matrix operator*(Matrix const & mat1, Matrix const & mat2)
{
    if (mat1.ncol() != mat2.nrow())
    {
        throw std::out_of_range(
            "the number of first matrix column "
            "differs from that of second matrix row");
    }

    Matrix ret(mat1.nrow(), mat2.ncol());

    for (size_t i=0; i<ret.nrow(); ++i)
    {
        for (size_t k=0; k<ret.ncol(); ++k)
        {
            double v = 0;
            for (size_t j=0; j<mat1.ncol(); ++j)
            {
                v += mat1(i,j) * mat2(j,k);
            }
            ret(i,k) = v;
        }
    }

    return ret;
}

std::ostream & operator << (std::ostream & ostr, Matrix const & mat)
{
    for (size_t i=0; i<mat.nrow(); ++i)
    {
        ostr << std::endl << " ";
        for (size_t j=0; j<mat.ncol(); ++j)
        {
            ostr << " " << std::setw(2) << mat(i, j);
        }
    }

    return ostr;
}

int main(int argc, char ** argv)
{
    std::cout << ">>> A(2x3) times B(3x2):" << std::endl;
    Matrix mat1(2, 3, std::vector<double>{1, 2, 3, 4, 5, 6});
    Matrix mat2(3, 2, std::vector<double>{1, 2, 3, 4, 5, 6});

    Matrix mat3 = mat1 * mat2;

    std::cout << "matrix A (2x3):" << mat1 << std::endl;
    std::cout << "matrix B (3x2):" << mat2 << std::endl;
    std::cout << "result matrix C (2x2) = AB:" << mat3 << std::endl;

    std::cout << ">>> B(3x2) times A(2x3):" << std::endl;
    Matrix mat4 = mat2 * mat1;
    std::cout << "matrix B (3x2):" << mat2 << std::endl;
    std::cout << "matrix A (2x3):" << mat1 << std::endl;
    std::cout << "result matrix D (3x3) = BA:" << mat4 << std::endl;

    return 0;
}
Build “Example code for matrix-matrix multiplication (ma03_matrix_matrix.cpp)”#
$ g++ ma03_matrix_matrix.cpp -o ma03_matrix_matrix -std=c++17 -O3 -g -m64
Example code for linear solver (la01_gesv.cpp)#
  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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>

#ifdef HASMKL
#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#else // HASMKL
#ifdef __MACH__
#include <clapack.h>
#include <Accelerate/Accelerate.h>
#endif // __MACH__
#endif // HASMKL

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol, bool column_major)
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
    }

    Matrix
    (
        size_t nrow, size_t ncol, bool column_major
      , std::vector<double> const & vec
    )
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
        (*this) = vec;
    }

    Matrix & operator=(std::vector<double> const & vec)
    {
        if (size() != vec.size())
        {
            throw std::out_of_range("number of elements mismatch");
        }

        size_t k = 0;
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = vec[k];
                ++k;
            }
        }

        return *this;
    }

    Matrix(Matrix const & other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(other.m_nrow, other.m_ncol);
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
    }

    Matrix & operator=(Matrix const & other)
    {
        if (this == &other) { return *this; }
        if (m_nrow != other.m_nrow || m_ncol != other.m_ncol)
        {
            reset_buffer(other.m_nrow, other.m_ncol);
        }
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
        return *this;
    }

    Matrix(Matrix && other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(0, 0);
        std::swap(m_buffer, other.m_buffer);
    }

    Matrix & operator=(Matrix && other)
    {
        if (this == &other) { return *this; }
        reset_buffer(0, 0);
        std::swap(m_nrow, other.m_nrow);
        std::swap(m_ncol, other.m_ncol);
        std::swap(m_buffer, other.m_buffer);
        m_column_major = other.m_column_major;
        return *this;
    }

    ~Matrix()
    {
        reset_buffer(0, 0);
    }

    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[index(row, col)];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[index(row, col)];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

    size_t size() const { return m_nrow * m_ncol; }
    double buffer(size_t i) const { return m_buffer[i]; }
    std::vector<double> buffer_vector() const
    {
        return std::vector<double>(m_buffer, m_buffer+size());
    }

    double * data() const { return m_buffer; }

private:

    size_t index(size_t row, size_t col) const
    {
        if (m_column_major) { return row          + col * m_nrow; }
        else                { return row * m_ncol + col         ; }
    }

    void reset_buffer(size_t nrow, size_t ncol)
    {
        if (m_buffer) { delete[] m_buffer; }
        const size_t nelement = nrow * ncol;
        if (nelement) { m_buffer = new double[nelement]; }
        else          { m_buffer = nullptr; }
        m_nrow = nrow;
        m_ncol = ncol;
    }

    size_t m_nrow = 0;
    size_t m_ncol = 0;
    bool m_column_major = false;
    double * m_buffer = nullptr;

};

std::ostream & operator << (std::ostream & ostr, Matrix const & mat)
{
    for (size_t i=0; i<mat.nrow(); ++i)
    {
        ostr << std::endl << " ";
        for (size_t j=0; j<mat.ncol(); ++j)
        {
            ostr << " " << std::setw(2) << mat(i, j);
        }
    }

    ostr << std::endl << " data: ";
    for (size_t i=0; i<mat.size(); ++i)
    {
        ostr << " " << std::setw(2) << mat.data()[i];
    }

    return ostr;
}

std::ostream & operator << (std::ostream & ostr, std::vector<double> const & vec)
{
    for (size_t i=0; i<vec.size(); ++i)
    {
        std::cout << " " << vec[i];
    }

    return ostr;
}

int main(int argc, char ** argv)
{
    const size_t n = 3;
    int status;

    std::cout << ">>> Solve Ax=b (row major)" << std::endl;
    Matrix mat(n, n, /* column_major */ false);
    mat(0,0) = 3; mat(0,1) = 5; mat(0,2) = 2;
    mat(1,0) = 2; mat(1,1) = 1; mat(1,2) = 3;
    mat(2,0) = 4; mat(2,1) = 3; mat(2,2) = 2;
    Matrix b(n, 2, false);
    b(0,0) = 57; b(0,1) = 23;
    b(1,0) = 22; b(1,1) = 12;
    b(2,0) = 41; b(2,1) = 84;
    std::vector<int> ipiv(n);

    std::cout << "A:" << mat << std::endl;
    std::cout << "b:" << b << std::endl;

#if !defined(HASMKL) || defined(NOMKL)
    {
        int nn = n;
        int bncol = b.ncol();
        int bnrow = b.nrow();
        int matnrow = mat.nrow();

        dgesv_( // column major.
            &nn // int * n: number of linear equation
          , &bncol // int * nrhs: number of RHS
          , mat.data() // double * a: array (lda, n)
          , &matnrow // int * lda: leading dimension of array a
          , ipiv.data() // int * ipiv: pivot indices
          , b.data() // double * b: array (ldb, nrhs)
          , &bnrow // int * ldb: leading dimension of array b
          , &status
          // for column major matrix, ldb remains the leading dimension.
        );
    }
#else // HASMKL NOMKL
    /*
     * "Note that using row-major ordering may require more memory and time
     * than column-major ordering, because the routine must transpose the
     * row-major order to the column-major order required by the underlying
     * LAPACK routine."  See:
     *
     * - https://www.netlib.org/lapack/lapacke.html#_array_arguments
     * - https://github.com/Reference-LAPACK/lapack/blob/2a39774316821436989cb755a59255cfa1ae9d99/LAPACKE/src/lapacke_dgesv_work.c#L63
     */
    status = LAPACKE_dgesv(
        LAPACK_ROW_MAJOR // int matrix_layout
      , n // lapack_int n
      , b.ncol() // lapack_int nrhs
      , mat.data() // double * a
      , mat.ncol() // lapack_int lda
      , ipiv.data() // lapack_int * ipiv
      , b.data() // double * b
      , b.ncol() // lapack_int ldb
      // for row major matrix, ldb becomes the trailing dimension.
    );
#endif // HASMKL NOMKL

    std::cout << "solution x:" << b << std::endl;
    std::cout << "dgesv status: " << status << std::endl;

    std::cout << ">>> Solve Ax=b (column major)" << std::endl;
    Matrix mat2 = Matrix(n, n, /* column_major */ true);
    mat2(0,0) = 3; mat2(0,1) = 5; mat2(0,2) = 2;
    mat2(1,0) = 2; mat2(1,1) = 1; mat2(1,2) = 3;
    mat2(2,0) = 4; mat2(2,1) = 3; mat2(2,2) = 2;
    Matrix b2(n, 2, true);
    b2(0,0) = 57; b2(0,1) = 23;
    b2(1,0) = 22; b2(1,1) = 12;
    b2(2,0) = 41; b2(2,1) = 84;

    std::cout << "A:" << mat2 << std::endl;
    std::cout << "b:" << b2 << std::endl;

#if !defined(HASMKL) || defined(NOMKL)
    {
        int nn = n;
        int b2ncol = b2.ncol();
        int b2nrow = b2.nrow();
        int mat2nrow = mat2.nrow();

        dgesv_( // column major.
            &nn // int * n: number of linear equation
          , &b2ncol // int * nrhs: number of RHS
          , mat2.data() // double * a: array (lda, n)
          , &mat2nrow // int * lda: leading dimension of array a
          , ipiv.data() // int * ipiv: pivot indices
          , b2.data() // double * b: array (ldb, nrhs)
          , &b2nrow // int * ldb: leading dimension of array b
          , &status
          // for column major matrix, ldb remains the leading dimension.
        );
    }
#else // HASMKL NOMKL
    status = LAPACKE_dgesv(
        LAPACK_COL_MAJOR // int matrix_layout
      , n // lapack_int n
      , b2.ncol() // lapack_int nrhs
      , mat2.data() // double * a
      , mat2.nrow() // lapack_int lda
      , ipiv.data() // lapack_int * ipiv
      , b2.data() // double * b
      , b2.nrow() // lapack_int ldb
      // for column major matrix, ldb remains the leading dimension.
    );
#endif // HASMKL NOMKL

    std::cout << "solution x:" << b2 << std::endl;
    std::cout << "dgesv status: " << status << std::endl;

    return 0;
}
Example code for eigenvalue problems (la02_geev.cpp)#
  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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>

#ifdef HASMKL
#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#else // HASMKL
#ifdef __MACH__
#include <clapack.h>
#include <Accelerate/Accelerate.h>
#endif // __MACH__
#endif // HASMKL

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol, bool column_major)
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
    }

    Matrix
    (
        size_t nrow, size_t ncol, bool column_major
      , std::vector<double> const & vec
    )
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
        (*this) = vec;
    }

    Matrix & operator=(std::vector<double> const & vec)
    {
        if (size() != vec.size())
        {
            throw std::out_of_range("number of elements mismatch");
        }

        size_t k = 0;
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = vec[k];
                ++k;
            }
        }

        return *this;
    }

    Matrix(Matrix const & other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(other.m_nrow, other.m_ncol);
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
    }

    Matrix & operator=(Matrix const & other)
    {
        if (this == &other) { return *this; }
        if (m_nrow != other.m_nrow || m_ncol != other.m_ncol)
        {
            reset_buffer(other.m_nrow, other.m_ncol);
        }
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
        return *this;
    }

    Matrix(Matrix && other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(0, 0);
        std::swap(m_buffer, other.m_buffer);
    }

    Matrix & operator=(Matrix && other)
    {
        if (this == &other) { return *this; }
        reset_buffer(0, 0);
        std::swap(m_nrow, other.m_nrow);
        std::swap(m_ncol, other.m_ncol);
        std::swap(m_buffer, other.m_buffer);
        m_column_major = other.m_column_major;
        return *this;
    }

    ~Matrix()
    {
        reset_buffer(0, 0);
    }

    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[index(row, col)];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[index(row, col)];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

    size_t size() const { return m_nrow * m_ncol; }
    double buffer(size_t i) const { return m_buffer[i]; }
    std::vector<double> buffer_vector() const
    {
        return std::vector<double>(m_buffer, m_buffer+size());
    }

    double * data() const { return m_buffer; }

private:

    size_t index(size_t row, size_t col) const
    {
        if (m_column_major) { return row          + col * m_nrow; }
        else                { return row * m_ncol + col         ; }
    }

    void reset_buffer(size_t nrow, size_t ncol)
    {
        if (m_buffer) { delete[] m_buffer; }
        const size_t nelement = nrow * ncol;
        if (nelement) { m_buffer = new double[nelement]; }
        else          { m_buffer = nullptr; }
        m_nrow = nrow;
        m_ncol = ncol;
    }

    size_t m_nrow = 0;
    size_t m_ncol = 0;
    bool m_column_major = false;
    double * m_buffer = nullptr;

};

std::ostream & operator << (std::ostream & ostr, Matrix const & mat)
{
    for (size_t i=0; i<mat.nrow(); ++i)
    {
        ostr << std::endl << " ";
        for (size_t j=0; j<mat.ncol(); ++j)
        {
            ostr << " " << std::setw(2) << mat(i, j);
        }
    }

    ostr << std::endl << " data: ";
    for (size_t i=0; i<mat.size(); ++i)
    {
        ostr << " " << std::setw(2) << mat.data()[i];
    }

    return ostr;
}

std::ostream & operator << (std::ostream & ostr, std::vector<double> const & vec)
{
    for (size_t i=0; i<vec.size(); ++i)
    {
        std::cout << " " << vec[i];
    }

    return ostr;
}

/*
 * See references:
 * * https://software.intel.com/en-us/mkl-developer-reference-c-geev
 * * https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/lapacke_dgeev_row.c.htm
 */
int main(int argc, char ** argv)
{
    const size_t n = 3;
    int status;

    std::cout << ">>> Solve Ax=lx (row major)" << std::endl;
    Matrix mat(n, n, /* column_major */ true);
    mat(0,0) = 3; mat(0,1) = 5; mat(0,2) = 2;
    mat(1,0) = 2; mat(1,1) = 1; mat(1,2) = 3;
    mat(2,0) = 4; mat(2,1) = 3; mat(2,2) = 2;
    std::vector<double> wr(n), wi(n);
    Matrix vl(n, n, /* column_major */ true), vr(n, n, /* column_major */ true);

    std::cout << "A:" << mat << std::endl;

#if !defined(HASMKL) || defined(NOMKL)
    {
        char jobvl = 'V';
        char jobvr = 'V';
        int nn = n;
        int vlnrow = vl.nrow();
        int vrnrow = vr.nrow();
        int matnrow = mat.nrow();
        int lwork = 4*n;
        std::vector<double> work(lwork);

        dgeev_( // column major.
            &jobvl
          , &jobvr
          , &nn // int * n: number of linear equation
          , mat.data() // double *: a
          , &matnrow // int *: lda
          , wr.data() // double *: wr
          , wi.data() // double *: wi
          , vl.data() // double *: vl
          , &vlnrow // int *: ldvl
          , vr.data() // double *: vr
          , &vrnrow // int *: ldvr
          , work.data() // double *: work array
          , &lwork // int *: lwork
          , &status
          // for column major matrix, ldb remains the leading dimension.
        );
    }
#else // HASMKL NOMKL
    status = LAPACKE_dgeev(
        LAPACK_COL_MAJOR // int matrix_layout
      , 'V' // char jobvl; 'V' to compute left eigenvectors, 'N' to not compute them
      , 'V' // char jobvr; 'V' to compute right eigenvectors, 'N' to not compute them
      , n // lapack_int n
      , mat.data() // double * a
      , mat.nrow() // lapack_int lda
      , wr.data() // double * wr
      , wi.data() // double * wi
      , vl.data() // double * vl
      , vl.nrow() // lapack_int ldvl
      , vr.data() // double * vr
      , vr.nrow() // lapack_int ldvr
    );
#endif // HASMKL NOMKL

    std::cout << "dgeev status: " << status << std::endl;

    std::cout << "eigenvalues:" << std::endl;
    std::cout << "      (real)      (imag)" << std::endl;
    std::cout << std::fixed;
    for (size_t i = 0 ; i < n ; ++i )
    {
        std::cout
            << "(  " << std::setw(9) << wr[i]
            << ",  " << std::setw(9) << wi[i]
            << ")" << std::endl;
    }
    std::cout << std::resetiosflags(std::ios_base::basefield);

    auto print_eigenvectors = [&wi](std::ostream & ostr, Matrix const & mat)
    {
        for (size_t i = 0 ; i < mat.nrow() ; ++i)
        {
            size_t j = 0;
            while (j < mat.ncol())
            {
                if ( wi[j] == (double)0.0 )
                {
                    ostr << "   " << std::setw(9) << mat(i, j);
                    ++j;
                }
                else
                {
                    auto eimag = mat(i, j+1);
                    ostr
                        << " ( " << std::setw(9) << mat(i, j)
                        << ", "  << std::setw(9) << eimag << ")";
                    if (eimag != 0) { eimag = -eimag; }
                    ostr
                        << " ( " << std::setw(9) << mat(i, j)
                        << ", "  << std::setw(9) << eimag << ")";
                    j += 2;
                }
            }
            ostr << std::endl;
        }
    };

    std::cout << "left eigenvectors:" << std::endl;
    print_eigenvectors(std::cout, vl);

    std::cout << "right eigenvectors:" << std::endl;
    print_eigenvectors(std::cout, vr);

    return 0;
}
Example code for eigenvalue problems for symmetric matrices (la03_syev.cpp)#
  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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>

#ifdef HASMKL
#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#else // HASMKL
#ifdef __MACH__
#include <clapack.h>
#include <Accelerate/Accelerate.h>
#endif // __MACH__
#endif // HASMKL

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol, bool column_major)
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
    }

    Matrix
    (
        size_t nrow, size_t ncol, bool column_major
      , std::vector<double> const & vec
    )
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
        (*this) = vec;
    }

    Matrix & operator=(std::vector<double> const & vec)
    {
        if (size() != vec.size())
        {
            throw std::out_of_range("number of elements mismatch");
        }

        size_t k = 0;
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = vec[k];
                ++k;
            }
        }

        return *this;
    }

    Matrix(Matrix const & other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(other.m_nrow, other.m_ncol);
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
    }

    Matrix & operator=(Matrix const & other)
    {
        if (this == &other) { return *this; }
        if (m_nrow != other.m_nrow || m_ncol != other.m_ncol)
        {
            reset_buffer(other.m_nrow, other.m_ncol);
        }
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
        return *this;
    }

    Matrix(Matrix && other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(0, 0);
        std::swap(m_buffer, other.m_buffer);
    }

    Matrix & operator=(Matrix && other)
    {
        if (this == &other) { return *this; }
        reset_buffer(0, 0);
        std::swap(m_nrow, other.m_nrow);
        std::swap(m_ncol, other.m_ncol);
        std::swap(m_buffer, other.m_buffer);
        m_column_major = other.m_column_major;
        return *this;
    }

    ~Matrix()
    {
        reset_buffer(0, 0);
    }

    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[index(row, col)];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[index(row, col)];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

    size_t size() const { return m_nrow * m_ncol; }
    double buffer(size_t i) const { return m_buffer[i]; }
    std::vector<double> buffer_vector() const
    {
        return std::vector<double>(m_buffer, m_buffer+size());
    }

    double * data() const { return m_buffer; }

private:

    size_t index(size_t row, size_t col) const
    {
        if (m_column_major) { return row          + col * m_nrow; }
        else                { return row * m_ncol + col         ; }
    }

    void reset_buffer(size_t nrow, size_t ncol)
    {
        if (m_buffer) { delete[] m_buffer; }
        const size_t nelement = nrow * ncol;
        if (nelement) { m_buffer = new double[nelement]; }
        else          { m_buffer = nullptr; }
        m_nrow = nrow;
        m_ncol = ncol;
    }

    size_t m_nrow = 0;
    size_t m_ncol = 0;
    bool m_column_major = false;
    double * m_buffer = nullptr;

};

std::ostream & operator << (std::ostream & ostr, Matrix const & mat)
{
    for (size_t i=0; i<mat.nrow(); ++i)
    {
        ostr << std::endl << " ";
        for (size_t j=0; j<mat.ncol(); ++j)
        {
            ostr << " " << std::setw(2) << mat(i, j);
        }
    }

    ostr << std::endl << " data: ";
    for (size_t i=0; i<mat.size(); ++i)
    {
        ostr << " " << std::setw(2) << mat.data()[i];
    }

    return ostr;
}

std::ostream & operator << (std::ostream & ostr, std::vector<double> const & vec)
{
    for (size_t i=0; i<vec.size(); ++i)
    {
        std::cout << " " << vec[i];
    }

    return ostr;
}

/*
 * See references:
 * * https://software.intel.com/en-us/mkl-developer-reference-c-syev
 * * https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/lapacke_dsyev_row.c.htm
 */
int main(int argc, char ** argv)
{
    const size_t n = 3;
    int status;

    std::cout << ">>> Solve Ax=lx (row major, A symmetric)" << std::endl;
    Matrix mat(n, n, /* column_major */ true);
    mat(0,0) = 3; mat(0,1) = 5; mat(0,2) = 2;
    mat(1,0) = 5; mat(1,1) = 1; mat(1,2) = 3;
    mat(2,0) = 2; mat(2,1) = 3; mat(2,2) = 2;
    std::vector<double> w(n);

    std::cout << "A:" << mat << std::endl;

#if !defined(HASMKL) || defined(NOMKL)
    {
        char jobz = 'V';
        char uplo = 'U';
        int nn = n;
        int matnrow = mat.nrow();
        int lwork = 3*n;
        std::vector<double> work(lwork);

        dsyev_(
            &jobz
          , &uplo
          , &nn // int * n: number of linear equation
          , mat.data() // double *: a
          , &matnrow // int *: lda
          , w.data() // double *: w
          , work.data() // double *: work array
          , &lwork // int *: lwork
          , &status
          // for column major matrix, ldb remains the leading dimension.
        );
    }
#else // HASMKL NOMKL
    status = LAPACKE_dsyev(
        LAPACK_COL_MAJOR // int matrix_layout
      , 'V' // char jobz;
            // 'V' to compute both eigenvalues and eigenvectors,
            // 'N' only eigenvalues
      , 'U' // char uplo;
            // 'U' use the upper triangular of input a,
            // 'L' use the lower
      , n // lapack_int n
      , mat.data() // double * a
      , mat.nrow() // lapack_int lda
      , w.data() // double * w
    );
#endif // HASMKL NOMKL

    std::cout << "dsyev status: " << status << std::endl;
    std::cout << "eigenvalues: " << w << std::endl;
    std::cout << "eigenvectors:" << mat << std::endl;

    return 0;
}
Example code for singular-value decomposition (la04_gesvd.cpp)#
  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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>
#include <algorithm>

#ifdef HASMKL
#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#else // HASMKL
#ifdef __MACH__
#include <clapack.h>
#include <Accelerate/Accelerate.h>
#endif // __MACH__
#endif // HASMKL

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol, bool column_major)
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
    }

    Matrix
    (
        size_t nrow, size_t ncol, bool column_major
      , std::vector<double> const & vec
    )
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
        (*this) = vec;
    }

    Matrix & operator=(std::vector<double> const & vec)
    {
        if (size() != vec.size())
        {
            throw std::out_of_range("number of elements mismatch");
        }

        size_t k = 0;
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = vec[k];
                ++k;
            }
        }

        return *this;
    }

    Matrix(Matrix const & other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(other.m_nrow, other.m_ncol);
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
    }

    Matrix & operator=(Matrix const & other)
    {
        if (this == &other) { return *this; }
        if (m_nrow != other.m_nrow || m_ncol != other.m_ncol)
        {
            reset_buffer(other.m_nrow, other.m_ncol);
        }
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
        return *this;
    }

    Matrix(Matrix && other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(0, 0);
        std::swap(m_buffer, other.m_buffer);
    }

    Matrix & operator=(Matrix && other)
    {
        if (this == &other) { return *this; }
        reset_buffer(0, 0);
        std::swap(m_nrow, other.m_nrow);
        std::swap(m_ncol, other.m_ncol);
        std::swap(m_buffer, other.m_buffer);
        m_column_major = other.m_column_major;
        return *this;
    }

    ~Matrix()
    {
        reset_buffer(0, 0);
    }

    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[index(row, col)];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[index(row, col)];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

    size_t size() const { return m_nrow * m_ncol; }
    double buffer(size_t i) const { return m_buffer[i]; }
    std::vector<double> buffer_vector() const
    {
        return std::vector<double>(m_buffer, m_buffer+size());
    }

    double * data() const { return m_buffer; }

private:

    size_t index(size_t row, size_t col) const
    {
        if (m_column_major) { return row          + col * m_nrow; }
        else                { return row * m_ncol + col         ; }
    }

    void reset_buffer(size_t nrow, size_t ncol)
    {
        if (m_buffer) { delete[] m_buffer; }
        const size_t nelement = nrow * ncol;
        if (nelement) { m_buffer = new double[nelement]; }
        else          { m_buffer = nullptr; }
        m_nrow = nrow;
        m_ncol = ncol;
    }

    size_t m_nrow = 0;
    size_t m_ncol = 0;
    bool m_column_major = false;
    double * m_buffer = nullptr;

};

std::ostream & operator << (std::ostream & ostr, Matrix const & mat)
{
    for (size_t i=0; i<mat.nrow(); ++i)
    {
        ostr << std::endl << " ";
        for (size_t j=0; j<mat.ncol(); ++j)
        {
            ostr << " " << std::setw(10) << mat(i, j);
        }
    }

    ostr << std::endl << " data: ";
    for (size_t i=0; i<mat.size(); ++i)
    {
        ostr << " " << mat.data()[i];
    }

    return ostr;
}

std::ostream & operator << (std::ostream & ostr, std::vector<double> const & vec)
{
    for (size_t i=0; i<vec.size(); ++i)
    {
        std::cout << " " << vec[i];
    }

    return ostr;
}

/*
 * See references:
 * * https://software.intel.com/en-us/mkl-developer-reference-c-gesvd
 * * https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/lapacke_dgesvd_row.c.htm
 */
int main(int argc, char ** argv)
{
    const size_t m = 3, n = 4;
    int status;

    std::cout << ">>> SVD" << std::endl;
    Matrix mat(m, n, /* column_major */ true);
    mat(0,0) = 3; mat(0,1) = 5; mat(0,2) = 2; mat(0, 3) = 6;
    mat(1,0) = 2; mat(1,1) = 1; mat(1,2) = 3; mat(1, 3) = 2;
    mat(2,0) = 4; mat(2,1) = 3; mat(2,2) = 2; mat(2, 3) = 4;
    std::vector<double> s(m), superb(m);
    Matrix u(m, m, /* column_major */ true);
    Matrix vt(n, n, /* column_major */ true);

    std::cout << "A:" << mat << std::endl;

#if !defined(HASMKL) || defined(NOMKL)
    {
        char jobu = 'A';
        char jobv = 'A';
        int mm = m;
        int nn = n;
        int matnrow = mat.nrow();
        int lwork = 5 * std::max(m, n);
        std::vector<double> work(lwork);

        dgesvd_( // column major.
            &jobu
          , &jobv
          , &mm // int *: m
          , &nn // int *: n
          , mat.data() // double *: a
          , &matnrow // int *: lda
          , s.data() // double *: s
          , u.data() // double *: u
          , &mm
          , vt.data() // double *: vt
          , &nn // int *: ldvt
          , work.data() // double *: work
          , &lwork // int *: lwork
          , &status
          // for column major matrix, ldb remains the leading dimension.
        );
    }
#else // HASMKL NOMKL
    status = LAPACKE_dgesvd(
        LAPACK_COL_MAJOR // int matrix_layout;
      , 'A' // char jobu;
      , 'A' // char jobvt;
      , m // lapack_int m
      , n // lapack_int n
      , mat.data() // double * a
      , mat.nrow() // lapack_int lda
      , s.data() // double * s
      , u.data() // double * u
      , u.nrow() // lapack_int ldu
      , vt.data() // double * vt
      , vt.nrow() // lapack_int ldvt
      , superb.data() // double * superb
    );
#endif // HASMKL NOMKL

    std::cout << "dgesvd status: " << status << std::endl;
    std::cout << "singular values: " << s << std::endl;
    std::cout << "u: " << u << std::endl;
    std::cout << "vt: " << vt << std::endl;

    return 0;
}
Example code for linear least-square problems (la05_gels.cpp)#
  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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>
#include <algorithm>

#ifdef HASMKL
#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#else // HASMKL
#ifdef __MACH__
#include <clapack.h>
#include <Accelerate/Accelerate.h>
#endif // __MACH__
#endif // HASMKL

class Matrix {

public:

    Matrix(size_t nrow, size_t ncol, bool column_major)
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
    }

    Matrix
    (
        size_t nrow, size_t ncol, bool column_major
      , std::vector<double> const & vec
    )
      : m_nrow(nrow), m_ncol(ncol), m_column_major(column_major)
    {
        reset_buffer(nrow, ncol);
        (*this) = vec;
    }

    Matrix & operator=(std::vector<double> const & vec)
    {
        if (size() != vec.size())
        {
            throw std::out_of_range("number of elements mismatch");
        }

        size_t k = 0;
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = vec[k];
                ++k;
            }
        }

        return *this;
    }

    Matrix(Matrix const & other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(other.m_nrow, other.m_ncol);
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
    }

    Matrix & operator=(Matrix const & other)
    {
        if (this == &other) { return *this; }
        if (m_nrow != other.m_nrow || m_ncol != other.m_ncol)
        {
            reset_buffer(other.m_nrow, other.m_ncol);
        }
        for (size_t i=0; i<m_nrow; ++i)
        {
            for (size_t j=0; j<m_ncol; ++j)
            {
                (*this)(i,j) = other(i,j);
            }
        }
        return *this;
    }

    Matrix(Matrix && other)
      : m_nrow(other.m_nrow), m_ncol(other.m_ncol)
      , m_column_major(other.m_column_major)
    {
        reset_buffer(0, 0);
        std::swap(m_buffer, other.m_buffer);
    }

    Matrix & operator=(Matrix && other)
    {
        if (this == &other) { return *this; }
        reset_buffer(0, 0);
        std::swap(m_nrow, other.m_nrow);
        std::swap(m_ncol, other.m_ncol);
        std::swap(m_buffer, other.m_buffer);
        m_column_major = other.m_column_major;
        return *this;
    }

    ~Matrix()
    {
        reset_buffer(0, 0);
    }

    double   operator() (size_t row, size_t col) const
    {
        return m_buffer[index(row, col)];
    }
    double & operator() (size_t row, size_t col)
    {
        return m_buffer[index(row, col)];
    }

    size_t nrow() const { return m_nrow; }
    size_t ncol() const { return m_ncol; }

    size_t size() const { return m_nrow * m_ncol; }
    double buffer(size_t i) const { return m_buffer[i]; }
    std::vector<double> buffer_vector() const
    {
        return std::vector<double>(m_buffer, m_buffer+size());
    }

    double * data() const { return m_buffer; }

private:

    size_t index(size_t row, size_t col) const
    {
        if (m_column_major) { return row          + col * m_nrow; }
        else                { return row * m_ncol + col         ; }
    }

    void reset_buffer(size_t nrow, size_t ncol)
    {
        if (m_buffer) { delete[] m_buffer; }
        const size_t nelement = nrow * ncol;
        if (nelement) { m_buffer = new double[nelement]; }
        else          { m_buffer = nullptr; }
        m_nrow = nrow;
        m_ncol = ncol;
    }

    size_t m_nrow = 0;
    size_t m_ncol = 0;
    bool m_column_major = false;
    double * m_buffer = nullptr;

};

std::ostream & operator << (std::ostream & ostr, Matrix const & mat)
{
    for (size_t i=0; i<mat.nrow(); ++i)
    {
        ostr << std::endl << " ";
        for (size_t j=0; j<mat.ncol(); ++j)
        {
            ostr << " " << std::setw(10) << mat(i, j);
        }
    }

    ostr << std::endl << " data: ";
    for (size_t i=0; i<mat.size(); ++i)
    {
        ostr << " " << mat.data()[i];
    }

    return ostr;
}

std::ostream & operator << (std::ostream & ostr, std::vector<double> const & vec)
{
    for (size_t i=0; i<vec.size(); ++i)
    {
        std::cout << " " << vec[i];
    }

    return ostr;
}

/*
 * See references:
 * * https://software.intel.com/en-us/mkl-developer-reference-c-gels
 * * https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/lapacke_dgels_row.c.htm
 */
int main(int argc, char ** argv)
{
    const size_t m = 4, n = 3;
    int status;

    std::cout << ">>> least square" << std::endl;
    // Use least-square to fit the data of (x, y) tuple:
    // (1, 17), (2, 58), (3, 165), (4, 360) to
    // the equation: a_1 x^3 + a_2 x^2 + a_3 x.
    Matrix mat(m, n, /* column_major */ true);
    mat(0,0) = 1; mat(0,1) = 1; mat(0,2) = 1;
    mat(1,0) = 8; mat(1,1) = 4; mat(1,2) = 2;
    mat(2,0) = 27; mat(2,1) = 9; mat(2,2) = 3;
    mat(3,0) = 64; mat(3,1) = 16; mat(3,2) = 4;
    std::vector<double> y{17, 58, 165, 360};
    // The equation f(x) = 3x^3 + 7^2x + 8x can perfectly fit the following
    // RHS:
    // std::vector<double> y{18, 68, 168, 336};

    std::cout << "J:" << mat << std::endl;
    std::cout << "y:" << y << std::endl;

#if !defined(HASMKL) || defined(NOMKL)
    {
        char trans = 'N';
        int mm = m;
        int nn = n;
        int nrhs = 1;
        int lwork = mm*nn + std::max(mm*nn, nrhs);
        std::vector<double> work(lwork);

        dgels_( // column major.
            &trans
          , &mm // int *: m
          , &nn // int *: n
          , &nrhs // int *: nrhs
          , mat.data() // double *: a for the 'J' matrix
          , &mm // int *: lda
          , y.data() // double *: the 'b' RHS buffer
          , &mm // int *: ldb
          , work.data() // double *: working buffer
          , &lwork // int *: size of working buffer
          , &status
          // for column major matrix, ldb remains the leading dimension.
        );
    }
#else // HASMKL NOMKL
    status = LAPACKE_dgels(
        LAPACK_COL_MAJOR // int matrix_layout
      , 'N' // transpose;
            // 'N' is no transpose,
            // 'T' is transpose,
            // 'C' conjugate transpose
      , m // number of rows of matrix
      , n // number of columns of matrix
      , 1 // nrhs; number of columns of RHS
      , mat.data() // a; the 'J' matrix
      , m // lda; leading dimension of matrix
      , y.data() // b; RHS
      , m // ldb; leading dimension of RHS
    );
#endif // HASMKL NOMKL

    std::cout << "dgels status: " << status << std::endl;
    std::cout << "a: " << y << std::endl;

    return 0;
}
Script to plot the least-square sample problem (la05_gels_plot.py)#
 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
#!/usr/bin/env python3

# [begin example]
import numpy as np
import matplotlib.pyplot as plt

poly = np.poly1d(np.array([5.35749, -2.04348, 12.5266, 0], dtype='float64'))
xin = np.array([1, 2, 3, 4], dtype='float64')
yin = np.array([18, 68, 168, 336], dtype='float64')
xp = np.linspace(1, 4, 100)

plt.rc('figure', figsize=(12, 8))
plt.plot(xin, yin, '.', xp, poly(xp), '-')
plt.xlim(0, 5)
plt.ylim(0, 400)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
# [end example]

import os
imagedir = os.path.join(os.path.dirname(__file__), '..', 'image')
imagebase = os.path.splitext(os.path.basename(__file__))[0] + '.png'
imagepath = os.path.join(imagedir, imagebase)
print('write to {}'.format(imagepath))
plt.savefig(imagepath, dpi=150)