{ "cells": [ { "cell_type": "markdown", "id": "61dead2f-10d1-4141-8892-898359e344bb", "metadata": {}, "source": [ "https://www.cs.usfca.edu/~galles/visualization/RotateScale2D.html\n", "# 1D Arrays for 2D Matrices\n", "\n", "\n", "## Representing a matrix\n", "\n", "A matrix is a 2d array of objects (numbers in our case) with m rows and n\n", "columns, with a total of $S=m\\times n$ elements. See\n", "\n", "- \n", "- \n", "- \n", "- ).\n", "\n", "Matrices are ubiquitous in physics and\n", "mathematics, and their manipulation are the core of a large proportion of\n", "computational and numerical problems. Efficient matrix manipulation and\n", "characteristic computing must be very performant and scalable. That is why,\n", "normally, you should always use a specialized library to solve typical problems\n", "like $Ax=b$, or to compute eigen values and eigen vectors, determinants and so\n", "on, but , outside the realm of linear algebra matrices are heavily used so you\n", "need to know how to represent and manipulate them.\n", "\n", "Normally, a matrix will be indexed as $A_{ij}$, with the first index\n", "representing the row and the second one the column. You must be aware that\n", "sometimes the indexes start at 1 and sometimes at 0. In c++ your indexes start\n", "at 0. Even though the matrix is a 2D array, it will be represented as an array\n", "of arrays in memory given that the ram memory is linear, as shown in\n", " and\n", " .\n", "For 2D arrays, you have basically two options:\n", "\n", "**row-major**\n", "
\n", " \"Image\n", "
From: \"https://eli.thegreenplace.net/images/2015/row-major-2D.png\"
\n", "
\n", "\n", "**column-major**\n", "
\n", " \"Image\n", "
From: \"https://eli.thegreenplace.net/images/2015/column-major-2D.png\"
\n", "
" ] }, { "cell_type": "markdown", "id": "06dd3f61-193a-49da-8a13-ce834efa8afd", "metadata": {}, "source": [ "For instance, the following code shows the row-major ordering in c/c++, using\n", "old/primitive arrays:\n", "\n", "\n", "```c++\n", " // REF: https://scicomp.stackexchange.com/questions/37264/how-much-space-store-a-matrix-of-numbers\n", " #include \n", " \n", " int main(int argc, char **argv) {\n", " const int N = 4;\n", " int m[N][N] = {0}; //4 by4 matrix with all zeros (default value for ints) allocated on the stack\n", " for (unsigned int i = 0; i < N; ++i) {\n", " for (unsigned int j = 0; j < N; ++j) {\n", " std::cout << &m[i][j] << \"\\t\";\n", " }\n", " std::cout << \"\\n\";\n", " }\n", " return 0;\n", " }\n", "```\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
0x16b2cadb80x16b2cadbc0x16b2cadc00x16b2cadc4
0x16b2cadc80x16b2cadcc0x16b2cadd00x16b2cadd4
0x16b2cadd80x16b2caddc0x16b2cade00x16b2cade4
0x16b2cade80x16b2cadec0x16b2cadf00x16b2cadf4
\n", "\n", "Run this on your system. Do you get the same addresses? are they continuous?" ] }, { "cell_type": "markdown", "id": "1d5008ab-be6c-4bd7-a608-270cbc865af4", "metadata": {}, "source": [ "## Defining matrices from primitives in C++\n", "\n", "In C++, matrices there could be defined in several ways :\n", "\n", "- Primitive array in the stack (contiguous but limited in size)\n", " \n", "\n", " ```c++\n", " double A[M][N];\n", " ...\n", " A[2][0] = 59.323;\n", " ```\n", "\n", "- Primitive dynamic array (not guaranteed to be contiguous in memory):\n", " \n", "\n", " ```c++\n", " double **A;\n", " A = new double *[M];\n", " for (int ii = 0; ii < M; ++ii) {\n", " A[ii] = new double [N];\n", " }\n", " //...\n", " for (int ii = 0; ii < M; ++ii) {\n", " delete A[ii];\n", " }\n", " delete [] A;\n", " A = nullptr;\n", " ```\n", " \n", "- A array of arrays : In the stack, limited size\n", " \n", "\n", " ```c++\n", " std::array, M> A;\n", " ```\n", " \n", "- A vector of vectors : not guaranteed to be contiguous in memory\n", " \n", " \n", " ```c++\n", " std::vector> A;\n", " A.resize(M);\n", " for (int ii = 0; ii < M; ++ii) {\n", " A[ii].resize(N);\n", " }\n", " ```\n", "\n", "All the previous are fine if you don't care about memory size or contiguous\n", "memory. But in scientific computing, applied to physics, we need to create\n", "performant data structures. For that reason we need to just ask for block of memory\n", "of size `M*N` and then view it as a 2D array using index tricks." ] }, { "cell_type": "markdown", "id": "ef67ec7e-0ad0-4ad4-83aa-98f0082ae2d5", "metadata": {}, "source": [ "\n", "\n", "## 1D and 2D view\n", "\n", "Being efficient when accessing memory is key to get a good performance. How are\n", "2D arrays stored in memory? Assume that you have the following 2D array (data\n", "inside each cell correspond to indexes)\n", "\n", "
\n", " \"Image\n", "
From: \"https://eli.thegreenplace.net/images/2015/row-major-2D.png\"
\n", "
\n", "\n", "As you can see, in the 1d array in memory, the first element (1d index 0) has the coordinates (0, 0), the second element (1d index 1) corresponds to (0, 1), and so on. The sixth 1d element, with 1d index 5, has 2D coordinates (1, 2). Can you devise a rule to relate the 1d and the 2d indexes? \n", "\n", "For a general matrix of size `MxN`, you have the following mapping from 2D to 1D\n", "\\begin{equation}\n", "id = iN + j,\n", "\\end{equation}\n", "while , from 1D to 2D\n", "\\begin{align}\n", "i &= id/N,\\\\\n", "j &= id%N.\n", "\\end{align}\n", "\n", "
\n", " \"Image\n", "
From: \"1d-2d-mapping.png\"
\n", "
" ] }, { "cell_type": "markdown", "id": "810da876-67a5-47c3-901c-4fcaef7dbcd3", "metadata": {}, "source": [ "This is a bi-dimensional version for filling a given matrix and transposing it.\n", "```c++\n", "#include \n", "#include \n", "#include \n", "\n", "void fill_matrix(std::vector & data, int m, int n);\n", "void print_matrix(const std::vector & data, int m, int n);\n", "void transpose_matrix(const std::vector & indata, int m, int n,\n", " std::vector & outdata);\n", "\n", "int main(int argc, char **argv)\n", "{\n", " const int M = std::stoi(argv[1]);\n", " const int N = std::stoi(argv[2]);\n", "\n", " std::vector array2d(M*N, 0.0);\n", "\n", " fill_matrix(array2d, M, N);\n", " print_matrix(array2d, M, N);\n", "\n", " std::vector array2d_T(M*N, 0.0);\n", " transpose_matrix(array2d, M, N, array2d_T);\n", " print_matrix(array2d_T, N, M);\n", "\n", " return 0;\n", "}\n", "\n", "void fill_matrix(std::vector & data, int m, int n)\n", "{\n", " for (int ii = 0; ii < m; ++ii) {\n", " for (int jj = 0; jj < n; ++jj) {\n", " data[ii*n + jj] = ii*n+jj; // A_(i, j) = i*n + j = id\n", " }\n", " }\n", "}\n", "\n", "void print_matrix(const std::vector & data, int m, int n)\n", "{\n", " for (int ii = 0; ii < m; ++ii) {\n", " for (int jj = 0; jj < n; ++jj) {\n", " std::cout << data[ii*n + jj] << \" \";\n", " }\n", " std::cout << \"\\n\";\n", " }\n", " std::cout << \"\\n\";\n", "}\n", "\n", "void transpose_matrix(const std::vector & datain, int m, int n,\n", " std::vector & dataout)\n", "{\n", " for (int ii = 0; ii < m; ++ii) {\n", " for (int jj = 0; jj < n; ++jj) {\n", " dataout[ii*m + jj] = datain[ii*n + jj]; // Error here\n", " }\n", " }\n", "}\n", "```\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
01234
56789
1011121314
1516171819
     
0123 
5678 
10111213 
15161718 
19000 
     
\n", "\n", "\n", "" ] }, { "cell_type": "markdown", "id": "e89c3ee4-3850-4cc9-afcd-ac4941e52b6b", "metadata": {}, "source": [ "## Exercises\n", "\n", "For the following exercises, and unless stated otherwise, the size of the matrices\n", "must be read from the command line and also the matrices must be filled randomly\n", "using a uniform random distribution in [-1, 1]. Also, the random seed must be\n", "read from the command line. Furthermore, the solutions must be modularized as\n", "follows:\n", "\n", "- **`matrix_utils.h`:** All the functions declarations and includes must come here. Each\n", " exercise implies one or more new declarations in this file.\n", "- **`matrix_utils.cpp`:** All the functions implementations come here.\n", "- **`main_XXXX.cpp`:** Each exercise will be solved in a main file where the\n", " necessary functions are called. The filename will change (the `XXXX` part).\n", "\n", "Then you will need to compile as\n", "\n", " g++ -std=c++17 -fsanitize=undefined,address matrix_utils.cpp matrix_XXXX.cpp\n", "\n", "Here is a basic example is show for the three files. Remember: the files\n", "`matrix_utils.cpp/h` will grow with new functions implemented there. The file\n", "`main_XXXX.cpp` will have a different name for each exercise and call different\n", "functions.\n", "\n", "- **`matrix_utils.h`:** Declarations.\n", "```c++ \n", "#pragma once\n", "#include \n", "#include \n", "#include \n", "#include \n", "void print_matrix(const std::vector & M, int nrows, int ncols);\n", "void fill_matrix_random(std::vector & M, const int nrows, const int ncols, const int seed);\n", "```\n", "- **`matrix_utils.cpp`:** Implementations\n", "```c++ \n", "#include \"matrix.h\"\n", "void print_matrix(const std::vector & M, int nrows, int ncols){\n", " std::cout.setf(std::ios::scientific);\n", " std::cout.precision(15);\n", "\n", " for (int ii = 0; ii < nrows; ++ii) {\n", " for (int jj = 0; jj < ncols; ++jj) {\n", " std::cout << M[ii*ncols + jj] << \" \"; // (ii, jj) = ii*NCOLS + jj\n", " }\n", " std::cout << \"\\n\";\n", " }\n", " std::cout << \"\\n\";\n", "}\n", "void fill_matrix_random(std::vector & M, const int nrows, const int ncols, const int seed){\n", " std::mt19937 gen(seed);\n", " std::uniform_real_distribution<> dis(-1, 1);\n", " for (int ii = 0; ii < nrows; ii++){\n", " for (int jj = 0; jj < ncols; jj++){\n", " M[ii*ncols + jj] = dis(gen);\n", " }\n", " }\n", "\n", "}\n", "```\n", "- **`main_example.cpp`:** Simple example\n", "```c++ \n", "#include \"matrix.h\"\n", "\n", "int main(int argc, char **argv) {\n", " // read size of matrix\n", " if (argc != 3) {\n", " std::cerr << \"Error. Usage:\\n\"\n", " << argv[0] << \" M N \\n\"\n", " << \"M : Rows\\n\"\n", " << \"N : Columns\\n\";\n", " return 1;\n", " }\n", " const int M = std::stoi(argv[1]);\n", " const int N = std::stoi(argv[2]);\n", "\n", " // create matrix A\n", " std::vector A(M*N);\n", "\n", " // fill matrix randomly\n", " fill_matrix_random(A, M, N, 0); // 0 == seed\n", "\n", " // print matrix\n", " print_matrix(A, M, N);\n", "\n", " return 0;\n", "}\n", "```\n" ] }, { "cell_type": "markdown", "id": "5ce31940-1138-43aa-8ce2-0d58b3958960", "metadata": {}, "source": [ "### Trace of a matrix \n", "Compute the trace of a random matrix using 1D indexing." ] }, { "cell_type": "markdown", "id": "b7502a04-a96b-44c9-997d-14a58bcdb832", "metadata": {}, "source": [ "### Hilbert matrix \n", "Create a function to compute the trace of that matrix. Fill the matrix as the Hilbert Matrix, $A_{ij} = 1/(i +j +1)$" ] }, { "cell_type": "markdown", "id": "d9b78d45-e2f0-4b79-8774-15bc4706899c", "metadata": {}, "source": [ "### Implement matrix-matrix multiplication \n", "A function must received the two input matrices and an output matrix and write there the result." ] }, { "cell_type": "markdown", "id": "2388f8b0-3aba-4eac-a628-366a75956394", "metadata": {}, "source": [ "### $A \\times A^t$\n", " Implement a function to compute the matrix multiplication of matrix with its transpose." ] }, { "cell_type": "markdown", "id": "6589276d-35cc-40b7-9431-b4c24a3adef4", "metadata": {}, "source": [ "### $A^p$\n", " Implement a function to compute the power of a matrix performing successive multiplications." ] }, { "cell_type": "markdown", "id": "3d820e86-04ea-4c0f-921f-be38c19ac5bf", "metadata": {}, "source": [ "### Idempotent matrix \n", " Implement a function that receives a matrix and a power and checks if it is\n", " idempotent under some precision, that is, if $|A^p - A| < \\epsilon$, where the\n", " norm is evaluated checking the inequality element-by-element." ] }, { "cell_type": "markdown", "id": "f3979927-2d27-4d2e-a203-01f9eeec293c", "metadata": {}, "source": [ "### Checking inverse \n", " Implement a function that receives two matrices and check if they are inverse\n", " of each other, that is , if $|AB - I| < \\epsilon$ and $|BA - I| < \\epsilon$, where the check is done on each element." ] }, { "cell_type": "markdown", "id": "30cb0114-d95c-44b8-968a-2d1a6586e50b", "metadata": {}, "source": [ "### Matrix commute \n", " Write a function that receives two matrices and checks if they commute under\n", " a given precision, that is , $|AB - BA| < \\epsilon$." ] }, { "cell_type": "markdown", "id": "3e9ae05c-9b7b-4696-949f-e1e199b054fd", "metadata": {}, "source": [ "### Orthogonal matrix \n", " Implement a function that receives a matrix and checks if it is orthogonal,\n", " that is, if $|A A^T - I| < \\epsilon$, where the norm is evaluated checking the\n", " inequality element-by-element." ] }, { "cell_type": "markdown", "id": "e28ac5a3-bb5a-4114-bb70-e40acd7611ba", "metadata": {}, "source": [ "### Matrix of complex numbers \n", " Implement a function that receives a matrix OF COMPLEX NUMBERS and checks if it is hermitian,\n", " that is, if $|A A^\\dagger - I| < \\epsilon$, where the norm is evaluated checking the\n", " inequality element-by-element. Use the `` header. $A^\\dagger$ means the\n", " transpose conjugate." ] }, { "cell_type": "markdown", "id": "b12487c3-bfe0-4350-bfc0-c8300ce00d42", "metadata": {}, "source": [ "### Pauli matrix \n", " Write a function that receives three components of a vectors and saves the\n", " corresponding Pauli-vector-matrix:\n", " " ] }, { "cell_type": "markdown", "id": "9d2d13fe-9db6-4eff-a0c7-1177f0f63bd9", "metadata": {}, "source": [ "### Complex matmul and Pauli matrices \n", " Extend your multiplication routine to use complex numbers and check the\n", " Pauli matrix commutation relations:\n", " " ] }, { "cell_type": "markdown", "id": "ee41a317-83c7-4b1e-ba46-b7f38812e149", "metadata": {}, "source": [ "### Vandermonde matrix\n", " Write a routine to fill a matrix as a Vandermonde matrix,\n", " , and then compute and\n", " print its trace Read the number of rows from the command line" ] }, { "cell_type": "markdown", "id": "2a306d8e-7528-4863-8f02-2d95bd679005", "metadata": {}, "source": [ "### Matrix polynomial \n", " Implement a function that receives a matrix and an array representing the\n", " coefficients of a polynomial, $a_0, a_1, \\ldots, a_n$, and evaluates that\n", " polynomial on the given matrix, $a_0 I + a_1 A + a_2 A^2 + \\ldots + a_n A^n$,\n", " and stores the result in an output matrix.\n", " " ] }, { "cell_type": "markdown", "id": "0881d0a4-c130-4bd1-ab27-f1ecea4f24dc", "metadata": {}, "source": [ "### Matmul time performance \n", " Matrix multiplication time:\n", " Explore how the total multiplication time grows with the matrix size. The\n", " following code is a template. You should read it and understand it carefully.\n", " Implement what is needed.\n", "\n", "```c++ \n", " #include \n", " #include \n", " #include \n", " #include \n", " #include \n", " #include \n", " #include \n", " \n", " void multiply(const std::vector & m1, const std::vector & m2, std::vector & m3);\n", " \n", " int main(int argc, char **argv) {\n", " // read parameters\n", " const int N = std::atoi(argv[1]);\n", " const int SEED = std::atoi(argv[2]);\n", " \n", " // data structs\n", " std::vector A(N*N, 0.0), B(N*N, 0.0), C(N*N, 0.0);\n", " \n", " // fill matrices A and B, using random numbers betwwen 0 and 1\n", " std::mt19937 gen(SEED);\n", " std::uniform_real_distribution<> dist(0.,1.);\n", " // lambda function: a local function that captures [] something, receives something (), and return something {}\n", " // See: https://www.learncpp.com/cpp-tutorial/introduction-to-lambdas-anonymous-functions/\n", " std::generate(A.begin(), A.end(), [&gen, &dist](){ return dist(gen); }); // uses a lambda function\n", " std::generate(B.begin(), B.end(), [&gen, &dist](){ return dist(gen); }); // uses a lambda function\n", " \n", " // multiply the matrices A and B and save the result into C. Measure time\n", " auto start = std::chrono::high_resolution_clock::now();\n", " multiply(A, B, C);\n", " auto stop = std::chrono::high_resolution_clock::now();\n", " \n", " // use the matrix to avoid the compiler removing it\n", " std::cout << C[N/2] << std::endl;\n", " \n", " // print time\n", " auto elapsed = std::chrono::duration(stop - start);\n", " std::cout << elapsed.count() << \"\\n\";\n", " \n", " return 0;\n", " }\n", " \n", " // implementations\n", " void multiply(const std::vector & m1, const std::vector & m2, std::vector & m3)\n", " {\n", " const int N = std::sqrt(m1.size()); // assumes square matrices\n", " delete this line and implement the matrix multiplication here\n", " }\n", "```\n", " \n", "When run as\n", "\n", "```bash\n", " ./a.out 64 10\n", "```\n", "\n", "you should get something like\n", "\n", "| | Value |\n", "|---|-------|\n", "| 1 | 14.6676 |\n", "| 2 | 0.002486 |\n", "\n", " Plot, in log-log scale, the time as a function of the matriz size $N$ in the\n", " range $N \\in \\{4, 8, 16, 32, 64, 128, 256, 512, 1024\\}$. Normalize all times by\n", " the time at $N = 4$." ] }, { "cell_type": "code", "execution_count": null, "id": "3c1bc069-b1f3-47d7-b8c6-554f43cbf932", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" } }, "nbformat": 4, "nbformat_minor": 5 }